text:
# Density functional study of Ni bulk, surfaces and the adsorbate systems Ni(111) (√3×√3)R30∘-Cl, and Ni(111)(2×2)-K

## I Introduction

Adsorption on surfaces is one of the main focuses in surface science due to the enormous importance of topics such as corrosion or catalysis. A huge number of experiments has been performed to understand these processes. In addition, theoretical studies have become capable of modelling the adsorption process.

The adsorption of alkali metals or haloges on metallic surfaces may be considered as two prototypes of adsorption processes. They are fairly simple and experimentally well studied. However, there have been quite a few surprises. For example, the occupancy of top sites by alkali metals as adsorbates was unexpected (for reviews, see references  and more general, reference ). Halogens are another type of adsorbate which is experimentally well studied. On close-packed (111) surfaces, they usually occupy the face centered cubic (fcc) hollow site, and in a few systems a partial occupancy (around 15-25%) of hcp (hexagonal close-packed) sites has also been observed . In addition, there is now a growing set of simulations for halogens on metallic surfaces (e.g. chlorine , bromine  or iodine ). Both for alkali metals and for halogens as adsorbates, there are general questions such as the geometry, binding energies, diffusion barriers, charge transfer and the binding mechanism.

Magnetic surfaces pose additional difficulties and questions such as how the magnetic moment changes after adsorption. It is thus an interesting question which insights simulation can give for these systems. Nickel surfaces and adsorption thereon have been well studied by experimental groups and thus are a prototype system for a computer simulation. For example, hydrogen on Ni(111) , benzene on Ni(111) , CH3 dehydrogenation on Ni(111)  or CO on Ni(110) adsorption have been simulated recently. The present article will deal with the simulation of alkali and halogen adsorption, with Cl and K as examples.

The intention of this article is to give a broad overview over nickel and the adsorbate systems. The computational technique is based on a local basis set of Gaussian type functions. As nowadays most calculations, especially for metals, rely on plane-waves, the following section II was devoted to the fairly simple systems of Ni bulk and the low index surfaces. This way, it will be demonstrated that results from a code based on Gaussian type orbitals reproduce the ones from plane-wave codes. In sections III and IV, the results for chlorine and potassium adsorption are presented, respectively. Finally, the results are summarized and discussed. Additional details on the calculations are given in the appendix.

## II Ni bulk and low index surfaces

### A Computational Parameters

A local basis set formalism was employed where the basis functions are Gaussian type orbitals centered at the atoms as implemented in the code CRYSTAL . The calculations were done with spin-polarization, using the gradient corrected functional of Perdew, Burke and Ernzerhof (PBE) , and in some additional calculations the gradient corrected functional by Perdew and Wang (PW) . These functionals perform actually nearly identical (see, e.g. references , or with the CRYSTAL code for Ag ), and also in the present article the results with both functionals are very similar. The PBE functional may be viewed as slightly favorable because it was designed to have a simpler form and derivation, and already revised versions exist (see, e.g., reference ).

The nickel basis set from reference  was employed, with the innermost $`[4s3p1d]`$ exponents being unchanged. In addition, 2 $`sp`$ functions with exponent 0.63 and 0.13 and a $`d`$-function with exponent 0.38 as optimized in calculations for the bulk were used, so that the basis set as a whole was of the size $`[6s5p2d]`$. The chlorine and potassium basis sets were the ones used in previous studies on Cl/Cu(111) and K/Cu(111), respectively. For the calculations on the free atoms, enhanced basis sets with additional diffuse functions were used . For the numerical integration, the potential was fitted with auxiliary basis sets (for Cl the same as in reference , for K as in reference , and for Ni the same as for Cu in reference ). A $`\stackrel{}{k}`$-point net with 16 $`\times `$ 16 points in the surface Brillouin zone and 16 $`\times `$ 16 $`\times `$ 16 points for the bulk was used. A smearing temperature of 0.001 $`E_h`$ was applied to the Fermi function ($`1E_h`$=27.2114 eV=315773 K). This value for the temperature was chosen relatively low to ensure that the magnetic moment is not artificially modified by a too high value. For more details and tests on the computational parameters and basis sets, see also the appendix.

The calculation of the charge and spin population is based on the Mulliken population, which is actually due to McWeeny (see also reference ). The Mulliken charge is defined as $`_{\mu A}_\nu P_{\mu \nu }S_{\nu \mu }`$ (with density matrix $`P`$, overlap matrix $`S`$ and $`A`$ the set of basis functions for which the population is computed). In the case of the spin population the density matrix for the corresponding spin polarization has to be used. This type of population analysis thus depends on the basis sets. It will be questionable if basis sets of different quality are used for various atoms: if one basis set is very small and the other very large, and especially if the larger basis set also has more diffuse exponents (i.e. smaller exponents, which results in a Gaussian with a relatively large spatial extension), the larger basis set may describe a part of the charge of the atom with the smaller basis set as well, and thus the Mulliken charge will be artificially large for the atom with the larger basis set. It is thus important that balanced basis sets are used. Concerning the Mulliken populations computed in this article, the population can be expected to be fairly reliable in the case when identical basis sets are used for the various atoms in a slab, so that the magnetic moments computed for the various layers of a nickel slab should not have a too large error due to the choice of the basis set. Also, comparing results for various adsorption sites should not be too unreasonable, because the basis sets for the systems are identical. This means, that, for example, a relative comparison of the Mulliken charge for chlorine on the fcc site with the charge on the top site should not be too unrealistic, although the absolute value of the Mulliken charge may have a larger error. A more difficult issue is however, to compare the chlorine charge for Cl/Ni(111) with the one for Cl/Ag(111) or Cl/Cu(111) because different basis sets for the metals are used and because of the different geometries.

The adsorbate systems Ni(111) $`(\sqrt{3}\times \sqrt{3})`$R30-Cl and Ni(111)$`(2\times 2)`$-K were modeled with a slab consisting of 3 nickel layers and the adsorbate layer, with the adsorbate layer on one side of the slab. The model is truly two-dimensional and not periodically repeated in the third dimension. This is different from the more common approach where periodicity in three dimensions is used, and a number of vacuum layers is required to separate the periodically repeated slabs. In the latter case, this also introduces an artificial electrical field due to the periodic boundary conditions in the third dimension, which is usually compensated for by an additional dipole field . The approach used by the CRYSTAL code for slabs avoids this artificial field as there are no periodic boundary conditions in the third dimension, and thus there is no correction required, and there is no parameter for the number of vacuum layers either. The simulation was performed for the fcc, hcp, bridge and top site (see figures 1 and 2). To ensure the stability of the results with respect to the number of layers, additional calculations with thicker slabs were performed for the systems Cl/Ni(111) (which is easier to converge than K/Ni(111)).

### B Properties of bulk Ni and clean Ni surfaces

Results for cohesive properties of bulk Ni are displayed in Table I. The computed values in this work reproduce essentially the ones obtained with the gradient corrected functional of Perdew and Wang, with the plane-wave pseudopotential method . The binding energy is overestimated, compared to the experiment. Also, the energy splitting between $`d^8s^2`$ and $`d^9s^1`$ is not well reproduced (the $`d^9s^1`$ state should be lower by -0.03 eV, after averaging the experimental energies over spin-orbit components). This is however, not a serious problem, as the orbital occupancy in the solid is different from the atom and thus the error cancellation does not work too well, and a density functional description may be too crude to describe these effects. The splitting between $`d^8s^2`$ and $`d^9s^1`$ state is a problem well known in quantum chemistry. The relativistic effect to the energy splitting between these states was computed to be of the order of 0.3 eV  (i.e. in a non-relativistic calculation as the present one, the splitting should be -0.33 eV instead of -0.03 eV) and thus is not negligible. With highly accurate configuration interaction schemes , the calculated splitting is -0.18 eV; i.e. there is still an error of 0.15 eV, compared to the experimental splitting, after relativistic effects were subtracted from the experimental value. The value computed here (-0.035 $`E_h`$ $`\widehat{=}`$ -0.95 eV); is in reasonable agreement with other results from literature, e.g. in the range of -1.1 to -1.2 eV, both at the level of the local density approximation (LDA) and the level of the generalized gradient approximation (GGA), with the PW functional . This difficulty in obtaining the correct energy splitting between $`d^{n1}s^1`$ and $`d^{n2}s^2`$ was also identified as a general problem for the transition metals , with the LDA, where a bias towards d$`{}_{}{}^{n1}s_{}^{1}`$ of $``$ 1 eV was found. With a different gradient corrected functional (B-LYP), an energy splitting of -0.82 eV was computed. There is, however, no reason why other properties should be affected by this problem.

The band structure is displayed in figure 3. It agrees well with the early tight binding  or augmented plane wave results, and thus supports the validity of the applied approach, relying on Gaussian type orbitals.

Surface formation energies (in $`E_h`$ per surface atom) of the clean, unrelaxed surfaces are displayed in Table II. The surfaces were modeled with a finite number of Ni layers (typically from 1 to 9). The PBE energies obtained here are within literature values. For the (111) surface, the surface energy was additionally computed with the Perdew-Wang functional, and is identical to the one with the PBE functional.

Finally, in Table III, the magnetic moments of the layers of the unrelaxed slabs are given. In all cases, they were computed using a 9 layer slab. They are in good agreement with literature values. For the (111) surface it is found that the second layer in the slab has the highest magnetic moment and not the top layer, in agreement with two of the published articles. For comparison, for the (111) surface, relaxation of the top layer was additionally taken into account. This resulted in a small inwards contraction by 0.03 Å, which is consistent with experimental  and theoretical  findings. The magnetic moments were virtually unchanged, and thus relaxation does not seem to explain why in one publication it was observed that the outermost layer should have the highest moment . For comparison, also the Perdew-Wang functional was used to compute the moments of the unrelaxed (111) surface, and the results of the PBE functional were essentially reproduced. As explained earlier, the calculation of the magnetic moments is based on the Mulliken population analysis. The error should, however, not be too large as identical basis sets are used for all the nickel atoms of the slabs.

## III Cl on the Ni(111) surface

As a first adsorbate system, the adsorption of Cl on the Ni(111) surface was investigated. This system was studied experimentally with LEED (low energy electron diffraction), Auger electron spectroscopy, work function measurements and flash desorption , surface extended X-ray absorption fine structure (SEXAFS) and soft X-ray standing wave (XSW) , and angle-resolved photoemission fine-structure measurements (ARPEFS). At low coverages there were indications of a $`(2\times 2)`$ pattern, followed by a $`(\sqrt{3}\times \sqrt{3})`$R30 structure, and a complex $`\left(\genfrac{}{}{0pt}{}{\mathrm{2\; 1}}{\mathrm{4\; 7}}\right)`$ pattern. The $`(\sqrt{3}\times \sqrt{3})`$R30 structure is the one which has been extensively studied experimentally and thus also the calculations will focus on this pattern. There was agreement that the fcc hollow was the preferred site (with an fcc:hcp population ratio of around 85:15, and the Cl-Ni bond length was measured to be 2.40 $`\pm `$ 0.04 Å, 2.33 $`\pm `$ 0.02 Å or 2.332 $`\pm `$ 0.006 Å. This corresponds to an interlayer distance of 1.92 Å, 1.83 Å or 1.837 . The first nickel layer was found to be practically unrelaxed in two experiments, whereas in the other article, a 5% contraction with respect to the bulk value was deduced .

It is thus an interesting issue to compare with results from the simulations. These results are presented in Table IV. A structural optimization was performed with 3 and 4 nickel layers. An additional single point calculation with 5 nickel layers was performed for fcc and hcp site, where the geometry from the thinner slabs was used, with one additional nickel layer at the distance corresponding to the bulk. The vertical positions of the chlorine adatoms and of the first nickel layer were optimized. Lateral relaxations are not possible for fcc and hcp site. As bridge and top site are much higher in energy, lateral relaxations were not taken into account for these sites either. Comparing the data for three and four layers, we see that the geometry is well converged, and the difference in the energy splitting between the sites is within what is expected as the numerical noise (e.g. fcc vs. hcp site, the energy splitting is 0.7 $`mE_h`$ with three nickel layers, 0.2 $`mE_h`$ with four layers, and 0.6 $`mE_h`$ with five layers). Firstly, we note that the fcc hollow is indeed found to be lowest in energy, with the hcp hollow being nearly degenerate at an energy difference of only 0.2-0.7 $`mE_h`$ ($``$ 60-220 K). Although the numerical noise is not negligible when energy differences on this small scale are computed (see appendix), this energy difference would be compatible with the simple thermodynamic estimate given in reference , with copper as substrate. The bridge site is higher in energy by $``$ 3 $`mE_h`$ (or more than 900 K), the top site by 18 $`mE_h`$, and thus thermodynamically fcc and hcp sites are clearly favored. The energy splitting between fcc (or hcp site), and the bridge site, may also be a crude estimate for the diffusion barrier from one threefold-hollow site to another one. These findings are actually very similar to those  for Cu(111)$`(\sqrt{3}\times \sqrt{3})`$R30-Cl and Ag(111)$`(\sqrt{3}\times \sqrt{3})`$R30-Cl, and also for the system Ag(111)$`(1\times 1)`$-O : on the (111) surface, the highly coordinated sites are clearly lowest in energy for adsorbates with a small radius. The bond lengths from chlorine to the nearest nickel $`d_{\mathrm{Cl}\mathrm{Ni}\mathrm{nn}}`$ also increase with increasing coordination number and the binding energy increases with increasing coordination number, as is to be expected from Pauling’s argument : one strong Cl-Ni bond in the case of the top site will be shorter than the three weaker bonds in the case of fcc and hcp, but the three bonds as a whole are energetically more favorable. The computed bond length agrees well with the values from the literature .

There are two different experimental results for the interlayer distance between first and second nickel layer. The computed value from this article agrees very well with the experimental from references , whereas the large contraction of 0.10 Å  with respect to the bulk appears to be unlikely. A smaller contraction would also be in line with simulations with Cu or Ag as substrate .

An effective radius for Cl on Ni(111) could be computed by subtracting an effective nickel radius of $`\frac{3.53}{2\sqrt{2}}`$Å from the bond length $`d_{\mathrm{Cl}\mathrm{Ni}\mathrm{nn}}`$ of 2.34 Å. This would result in a radius of 1.09 Å, in agreement with the computed values  for Cl on Cu(111) (1.12 Å) or Cl on Ag(111) (1.17 Å). Thus the radius for Cl is intermediate between the atomic (0.99 Å) and ionic (Cl-; 1.81 Å) . The bond length $`d_{\mathrm{Cl}\mathrm{Ni}\mathrm{nn}}`$ is also in agreement with experimental values for systems with similar atomic radii for the substrate (e.g., bond lengths of 2.39 Å for the substrate Cu(111), 2.38 Å for Rh(111), 2.39 Å for Pd(111)).

In Table V, the Cl population is displayed. The data is extracted from the calculations where three nickel layers were used to model the substrate, but it is virtually identical to the data obtained with four nickel layers. The negative charge is largest for the top site, smaller for the bridge site and smallest for hcp and fcc site. This is consistent with the position of the 3$`s`$ core eigenvalue relative to the Fermi level, which is lowest for hcp and fcc, increasing for bridge and largest for the top site (because the negative charge of the chlorine adsorbate is increasing and thus the core levels are destabilized). Note that the core level for majority spin is also slightly lower (by $``$ 0.002 $`E_h`$) than for minority spin. The density of states is visualized in figure 4, projected on the chlorine basis functions, and projected on all basis functions (i.e. the total density of states). The different occupancy of the two spin states becomes visible. Note that the chlorine spin is mainly due to the different population of the second and third $`p`$-shell, and not of the fourth (if we order the shells according to the size of the exponents, starting with the highest). This means that the outermost and thus most diffuse exponent does not carry a huge spin and thus the chlorine spin is not an artifact of the Mulliken population (if the spin was on the most diffuse function only, one might argue that this basis function also describes the spin on the nickel atoms and thus the Mulliken population analysis would be questionable).

The computed Mulliken charge is also in agreement with the experimental estimates of 0.1 to 0.2 $`|e|`$ derived from work function measurements. The overlap population $`_{\mu A}_{\nu B}P_{\mu \nu }S_{\nu \mu }`$ (with density matrix $`P`$, overlap matrix $`S`$ and $`A,B`$ the set of basis functions of atoms $`A`$ and $`B`$ for which the overlap population is computed) is 0.09 for Cl and the nearest nickel atoms. This may indicate that there is no evidence for a strong covalent contribution to the binding mechanism, similar to non-magnetic substrates , and now also in the case of a magnetic substrate. For comparison, the overlap population between two nearest nickel atoms in the top nickel layer is 0.10. However, this data should not be over-interpreted, and the overlap population should rather be viewed as a rough qualitative measure and not a quantitative measure for covalency. Another indication that the covalency is not strong might be that in the case of the bridge site, where $`p_x`$ and $`p_y`$ orbitals do not have the same occupancy, the orbital with the larger overlap with the neighboring nickel atoms ($`p_x`$ in our choice of geometry) has a slightly smaller population than $`p_y`$ (the total Cl $`p`$ population is 3.73 for $`p_x`$, versus 3.82 for $`p_y`$).

Charge density, charge density difference and spin density are visualized in figure 5. As usual, the charge density difference is more informative as the charge density. The transfer from the top nickel layer to the chlorine atom becomes obvious. In contrast to the charge density, the spin density also offers more information: it is very interesting that the spin density varies relatively strong still for the atoms in the third layer. This is confirmed by the data from Mulliken population analysis for the spin (table VI). The magnetic moments depend relatively strong on the number of layers, and the analysis has been performed by extracting the data from a single-point calculation for a slab with 5 nickel layers, with an identical geometry as for the 3-layer slab, i.e. the additional layers were added at a distance as in the bulk.

The data is presented in such a way that one number is given, if all the atoms in a layer have the same distance to the chlorine atom, and two numbers, if there are two different distances for the atoms in one layer possible. For example, in the case of the fcc site, looking at the atoms in the first or second layer, all three nickel atoms within a supercell have the same distance to the adsorbed chlorine. However, in the third layer there will be one atom which is closer than the other two. The population of this closer atom (0.59) is displayed in the left column of the results for the $`3^{rd}`$ layer, the data for the other two atoms (0.64) in the right column.

It is interesting to note that those nickel atoms in the first layer which are closer to the chlorine atom have a reduced magnetic moment. This is strongest in the case of the top site where the atom vertically under the chlorine has its population reduced to 0.55, versus a moment of 0.76 for the other two atoms. It is also apparent for the bridge site (0.66 vs 0.73), and in the case of fcc and hcp site the data is identical because the nickel atoms in the first layer are identical by symmetry.

In the second layer, the spin population is higher for those atoms which are closer to the adsorbate (this is possible for the hcp and bridge structure). In the third layer, again those atoms which have a shorter distance to the adsorbate, have there moment reduced. These oscillations are still visible in the fourth and fifth layer.

Although this data depends relatively strong on the number of layers, it becomes apparent that there are oscillations in the spin density. It is also interesting to compare with the total charge: the latter quantity is virtually constant within the layers, from the second layer on, and has a value identical to that of a clean slab from the third layer on (also table VI). Only in the first layer, the atom(s) closer to the chlorine adsorbate have a slightly higher charge (i.e. more electronic charge) than the atom(s) which are further apart (visible for bridge and top site).

Finally, the chlorine adsorbate has a magnetic moment of $``$ 0.1 for all sites (the smallest value for the top site). The moment is always parallel to the nickel spin.

It is interesting to note that in the simulations for CO on Ni(110) , it was also found that the moment of the atom(s) closest to the CO was reduced.

## IV K on the Ni(111) surface

In addition, the system Ni(111)$`(2\times 2)`$-K was studied with this approach. Experimentally, in a LEED study at various coverages, this was found to be the only commensurate structure . In a further LEED study, the top site was found to be occupied, with a K-Ni distance of 2.82$`\pm `$ 0.04 Å, a vertical rumpling of the first Ni layer of 0.12 Å, and a lateral displacement of 0.06 Å of the three nickel atoms in the top layer which are not vertically under the K adsorbate. The adsorption site was later confirmed in SEXAFS measurements, and a larger bond length of 2.92$`\pm `$ 0.02 Å was obtained. With the same technique, also the bond length for the system K/Cu(111) was computed and deduced to be by 0.06 $`\pm `$ 0.04 Å larger. From ARPEFS, a top site adsorption with K-Ni bond-length of 3.02 $`\pm `$ 0.01 Å was obtained. In addition, the first to second layer spacing was reduced by 0.13 Å to 1.90 $`\pm `$ 0.04 Å, but no substrate rumpling in the top layer was observed. Lateral displacements were not found either (0.00 $`\pm `$ 0.09 Å) in this experiment. Finally, photoelectron diffraction (PhD)  resulted in a bond length of 2.87 $`\pm `$ 0.06 Å, and a reduced first to second layer distance of 1.86 $`\pm `$ 0.06 Å was observed, with virtually no rumpling in the first layer (0.01 $`\pm `$ 0.09 Å).

The experiments thus consistently support the top site as the adsorption site, but the geometrical parameters are not identical. Moreover, a comparison with the system K/Cu(111) is very interesting, so that the influence of the partially occupied $`d`$-bands can be investigated.

The results of the simulations are summarized in Table VII, with geometrical parameters as defined in figure 6. Initially, the potassium position and the position of the top nickel layer was optimized (i.e. the uniform relaxation of the top nickel layer). In a second step, a vertical rumpling of the nickel atoms in the top layer was allowed, as this is expected to be important for the top site occupancy. A further refinement would have been the additional possibility of lateral relaxations. These lateral relaxations were, however, found to be of minor importance for the related system K/Ag(111), where the energy was lowered by the order of typically 0.1 - 0.2 $`mE_h`$, when these relaxations were allowed. Thus, lateral relaxations were taken into account additionally only for the top site, which was found to be lowest in energy, and in order to compare with the experimental data for the lateral displacement.

Firstly, the top site is found to be the most favorable site. This is in agreement with the experimental findings . By disabling the possibility of substrate rumpling, it is also demonstrated that substrate rumpling is crucial for top site occupancy, as it lowers the energy by 1.6 $`mE_h`$, but only 0.4 $`mE_h`$ for the other sites. These findings are actually similar to the system K/Cu(111) . The top site thus becomes the most stable site because the nickel atom underneath is pushed into the substrate relative to its neighbors and the potassium with its relatively large radius thus also overlaps more with the second nearest nickel neighbors.

The computed bond length (2.79 Å) is also very similar to K/Cu(111) (2.83 Å)  so that these calculations do not indicate a major difference compared to Cu. Experimentally, using the same technique for both adsorbate systems, it was measured to be slightly shorter and it was speculated that this was due to greater adatom-substrate interactions involving the partially filled Ni $`d`$ bands. The possibility of a lateral displacement of the nickel atoms as observed experimentally (0.06 Å) was investigated. However, this could not be confirmed and the energy was found to be lowest for the case where this lateral displacement is exactly zero.

Bridge, fcc and hcp site are found to be energetically degenerate. The fact that all the sites are energetically very close (e.g. much closer than in the case of chlorine) is consistent with earlier findings  and was explained due to the large radius of the adsorbate which makes it experience only a small substrate electron density corrugation . This small energy splitting between the various sites is also consistent with the experiment because of various observations: firstly, the SEXAFS data indicated enhanced in-plane motion resulting from a wide and shallow potential well. Secondly, no SEXAFS signal was observed for low coverages below 0.13 monolayers. It was thus argued that lateral motion was essentially hindered by the presence of other adatoms. Finally, with increasing temperature (already in the range of 120-145 K) the SEXAFS amplitude was found to decrease. This was explained due to the occupancy of other sites, which leads to a destructive interference.

The bond length $`d_{\mathrm{K}\mathrm{Ni}\mathrm{nn}}`$ increases with the number of nearest nickel neighbors in the order top, bridge, fcc and hcp. The effective radius of the potassium adsorbate is 2.79-$`\frac{3.53}{2\sqrt{2}}`$ Å=1.54 Å for the top site. This number is in good agreement with experimental values for potassium when occupying the top site , or the computed value for K/Cu(111) (1.55 Å).

Looking at the populations in Table V, we see that they ($``$ 0.3 $`|e|`$) are very similar for all considered sites, only for the top site the charge is slightly less positive. In agreement with this, the core eigenvalues of the $`3s`$ and $`3p`$ orbital are virtually independent of the site. The core levels are identical for both spin polarizations which agrees with the finding that the magnetic moment of the potassium adatoms is negligible (see below). In the calculations for K/Ag(111) , it was recently shown that the adsorbate charge and thus also the core eigenvalue (and other properties such as the bond length) change with the coverage. Such a change in the core eigenvalues was observed experimentally for alkali metals on metallic surfaces, e.g. for the systems Na/Cu(111) and Na/Ni(111) , or for alkali metals on Ru(100)  or W(110) .

The overlap population is $``$ 0.03 for K and the nearest nickel atom, so that there is also in this case no evidence for a strong covalent contribution to the binding mechanism.

The density of states is displayed in figure 7. The integrated density of states, when projected on the potassium basis functions, differs so little for majority and minority spin that there is nearly no magnetic moment. For the top site, the spin is even antiparallel to the nickel sites, in contrast to the other sites (and in contrast to chlorine).

Again, the charge density is not very informative, and it is more interesting to look at the charge density difference (figure 8). Electrons have flown from the adsorbate mainly to the top nickel layer (and a relatively large part to the atom vertically under the potassium adatom). Thus, the charge of the potassium overlayer has decreased. It is apparent that the spin density of the nickel atom vertically under the top atom is slightly lower than the spin density of its neighbors in the same layer, as the contours are closer to the nucleus. This is confirmed in the Mulliken spin population (table VIII). Again, those atoms in the top layer which are closest to the potassium adsorbate, have their moment reduced. No noteworthy variation is found in the second and third nickel layer, and the values for the spin population are very close to the ones for the clean Ni slab. The moment of the potassium adsorbate is negligible, as its magnitude is $``$ 0.003. Also, the charge only varies in the first layer, and already the second and third layer have charges which are identical to that of a clean slab.

There are thus no oscillations of the spin population beyond the first nickel layer, probably because the potassium magnetic moment is much smaller than that of chlorine.

## V Summary

We have studied Ni bulk, the clean low index surfaces and the adsorption of chlorine and potassium on the Ni(111) surface.

The results for the adsorption geometry confirm the experimental findings. It turns out that there is no huge difference to copper as substrate; the preferred adsorption site and energies are very similar.

The reason for the preference of the top site for K is substrate rumpling which helps to increase the overlap of the K adsorbate with the second nearest nickel neighbors. This rumpling is thus crucial for the site preference. In the case of chlorine, the threefold hollow sites are very close in energy. In general, the adsorbate atoms prefer sites with the highest possible coordination number, or they try to increase the overlap to neighboring atoms by substrate rumpling.

There appear to be no strong covalent contributions to the binding so that the mechanism appears to be mainly an ionic bond (although the charge transfer is rather small), and in the case of potassium, additionally a metallic binding mechanism. However, the quantitative analysis of the binding mechanism is certainly difficult. The computed charges seem to be fairly reliable, but the overlap population is probably only a qualitative tool to estimate the degree of covalency.

Charges and core levels were found to be consistent, i.e. with increasing negative charge of the adsorbate, the core levels of the adsorbate are slightly destabilized. It is also demonstrated that the magnetic moments computed with the help of the Mulliken population analysis agree well with data from other schemes, e.g. pseudopotential plane wave. The calculation of work functions, however, seems to be difficult with a local basis set and the results are strongly basis-set dependent.

Finally, the magnetic moment of those nickel atoms closest to the adsorbate are always reduced. Oscillations of the spin population within one nickel layer are notable for the first three layers for Cl/Ni(111), but only for the first nickel layer for K/Ni(111), probably because of the very small magnetic moment of the potassium adsorbate, compared to chlorine.

## VI acknowledgments

All the calculations were performed at the computer centre of the TU Braunschweig (Compaq ES 45).

## A Test of computational parameters

The calculations require various sets of parameters, with the most important being: a finite number of $`\stackrel{}{k}`$-points, a finite temperature to smoothen the integrand and thus to facilitate the numerical integration of the exchange-correlation functional, basis set parameters, and truncation schemes to reduce the number of matrix elements (i.e. integrals of operators such as kinetic energy, nuclear attraction, coulomb and exchange, with the basis functions). In this first part of the appendix, the stability of the results with respect to the variation of parameters concerning the grid and the parameters concerning the selection of the integrals is tested. For further tests, see also the extensive tests performed for lithium , or tests for copper and silver ($`\stackrel{}{k}`$-points, smearing temperature). In tables IX and X, the binding energies are computed with the default grid and truncation parameters, and with an improved grid and two sets of stricter and thus more accurate truncation parameters.

The improved grid has roughly six times more sampling points than the default grid for the various structures. As is displayed in the tables, the grid has virtually no impact on the energy splitting. For Cl/Ni(111), the binding energy changes by at most 0.2 $`mE_h`$ for fcc, hcp and bridge site, and only for the top site an energy change of 0.9 $`mE_h`$ is observed, with respect to the default grid. This has, however, no impact as the top site is much higher in energy anyway. For K/Ni(111) the absolute value of the binding energy is shifted by 0.7 $`mE_h`$, and the relative energies of the various sites are identical.

There are two truncation parameters in calculations when no Fock exchange is involved: a general one (”ITOL 1”) which leads to a neglect of all integrals with an overlap below a certain threshold (default: 10-6), and a second parameter (”ITOL 2”) which leads to the evaluation of the coulomb integrals at a lower level of accuracy by means of a multipolar expansion, if the overlap is below another threshold (default: 10-6). As the results turned out to be more sensitive to these parameters, two sets of higher thresholds were used (10-7 for ITOL 1 and 2, and 10-8 for ITOL 1 and 2). In the case of chlorine, the largest error is found for the splitting between fcc and hcp site: it ranges from 0.7 to 2.1 $`mE_h`$, for the various parameters, with the fcc site always being more favorable. The bridge site is always clearly higher than the fcc site by 3.3 to 4.2 $`mE_h`$, and the top site is much higher in energy. For K/Ni(111), the impact of these parameters is smaller, and for example, the splitting from the top site to the site closest in energy ranges from 0.5 $`mE_h`$ to 0.8 $`mE_h`$.

As a whole, the impact of these computational parameters is certainly not negligible, but the conclusions are not affected. The fcc site is favored for Cl/Ni(111), with the hcp site being nearly degenerate. In the case of K/Ni(111), the top site is favored, with the other sites being very close in energy.

## B Basis set dependence of the results

The choice of the basis set is of high importance for the quality of the results. In this part of the appendix, the results of various tests are displayed, to investigate the dependence of the results on the basis set. These tests focus on nickel bulk, the clean surfaces, and Cl/Ni(111).

Firstly, one might consider reoptimizing the exponents for the clean surface, for example keeping the exponents of the inner layer like in the bulk and reoptimizing the exponents of the outermost layers. This would result (with a slab with three layers) in exponents of 0.12 ($`sp`$, instead of 0.13 for the bulk) and 0.38 ($`d`$, as in bulk nickel) for the outer layers. This is thus a tiny change only and the results are virtually identical, for example the surface energy per atom changes by 0.6 $`mE_h`$ which is negligible. If an outermost $`sp`$-exponent of 0.12 instead of 0.13 is used for the bulk, exactly the same results as in table I are obtained.

There is however a property which depends extremely strong on the basis set: the computed work function changes drastically when the outermost exponent is changed by a tiny value. For example, with the $`sp`$-exponent used for all the calculations (0.13), the work function of the Ni(100) surface is 0.142 $`E_h`$ (exp.: 0.192 $`E_h`$), of the Ni(110) surface 0.129 $`E_h`$ (exp.: 0.185 $`E_h`$), and of the Ni(111) surface 0.155 $`E_h`$ (exp.: 0.197 $`E_h`$). This data is practically stable with respect to the number of layers (i.e. 3 layers are sufficient to obtain a stable number). When, however, for example, the outermost $`sp`$-exponent of the outer layers is changed from 0.13 to 0.12, the work function of the Ni(111) surface changes by 0.01 $`E_h`$ to 0.165 $`E_h`$, and it further increases with smaller exponents (e.g. 0.183 $`E_h`$, if the outermost $`sp`$-exponent of the outer layers is 0.10 instead of 0.13). Also, increasing this exponent leads to a lower work function. The same problem shows also up with other metals, such as for example copper or silver (and also when the exponents of all the layers and not just the outermost layers are changed). It appears thus that the work function is a quantity which is not well described with a local basis set. The work functions may be qualitatively reasonable, but a quantitative comparison seems unadvised. Thus, also the work function data computed for Ag and the adsorption thereon  has to be considered with caution. A problem with computing work functions with a local basis set was already discussed earlier .

The results thus indicate that the only property which appears to be strongly basis set dependent is the work function (total energy, surface energy, bulk modulus and lattice constant are essentially stable with respect to small variations in the basis set). Mulliken charge and magnetic moment are stable if an identical basis is used for all the nickel atoms in a slab. Of course, if basis sets of different type are used for the inner layer of a slab and the outermost layers, the Mulliken population also gets less reliable.

Further tests were performed for the adsorbate system Cl/Ni(111). One test was done by choosing a more tight outermost $`sp`$-exponent of 0.15 instead of 0.09 for chlorine. When a structural optimization is performed, virtually the same structural data as in table IV is obtained (the maximum deviation is 0.01 Å). The binding energy is lower by $`3mE_h`$ because of the poorer chlorine basis set. However, what is important is that the energy splitting for the various sites remains practically identical (the binding energy is -0.1326 $`E_h`$ for the fcc site, -0.1320 $`E_h`$ for the hcp site, -0.1292 $`E_h`$ for the bridge site, and -0.1151 $`E_h`$ for the top site). Also, properties such as the population, magnetic moment and the relative position of the core levels with respect to the Fermi energy are essentially stable with respect to this change of the basis set (but again the work function changes strongly).

As a whole, the basis set needs of course to be carefully chosen. When this is done, the results for the geometry, energies, core eigenvalues, populations, and magnetic moments are reliable and not too strongly basis set dependent. The accurate determination of the work function, however, appears to be difficult with a local basis set.


format: text

text:
# Absence of magnetic order in Yb3Ga5O12: relation between phase transition and entropy in geometrically frustrated materials

## Abstract

From muon spin relaxation spectroscopy experiments, we show that the sharp peak ($`\lambda `$ type anomaly) detected by specific heat measurements at 54 mK for the ytterbium gallium garnet compound, Yb3Ga5O12, does not correspond to the onset of a magnetic phase transition, but to a pronounced building up of dynamical magnetic pair-correlations. Beside the $`\lambda `$ anomaly, a broad hump is observed at higher temperature in the specific heat of this garnet and other geometrically frustrated compounds. Comparing with other frustrated magnetic systems we infer that a ground state with long-range order is reached only when at least $`1/41/3`$ of the magnetic entropy is released at the $`\lambda `$ transition.

A magnetic phase transition to a long-range order occurs at low temperature for most crystallographically ordered compounds containing a three dimensional lattice of magnetic ions. Its classical signature is a sharp anomaly in the magnetic specific heat (SH) at temperature $`T_\lambda `$ which corresponds to the phase transition temperature; see e.g. Ref. Morrish (1965). In this letter, for convenience, we label such anomaly a $`\lambda `$-anomaly although it may not have all the properties attributed to such an anomaly Stanley (1971).

However, for some particular lattice structures, the long-range magnetic ordering of the magnetic ions may be impeded by their geometric arrangement which gives rise to frustration of their magnetic interactions Villain (1979); Diep (1994). Most of the experimental studies focus nowadays on pyrochlore and gallium garnet compounds, R2T2O7 and R3Ga5O12 respectively. R denotes a rare earth atom and T a transition element. The R ions are arranged on a motif of corner sharing tetrahedra for the pyrochlore structure. In the garnet case, the R atoms form two interpenetrating, non-coplanar, corner sharing triangular sublattices. For an experimental review, see e.g. Ref. Ramirez (2001).

The results from SH measurements are particularly intriguing. Although $`|\theta _{\mathrm{CW}}|/T_\lambda 1`$, where $`\theta _{\mathrm{CW}}`$ is the Curie-Weiss temperature, long-range order may not be present below $`T_\lambda `$.

For example, let us review data on pyrochlore compounds. According to SH, Gd2Ti2O7 has two magnetic phase transitions at 1 K and $``$ 0.75 K and no defined structure in the SH data is observed above 1 K Ramirez et al. (2002). The detection of magnetic Bragg reflections from neutron scattering at 50 mK shows that a magnetic structure is established Champion et al. (2001). Similarly, Gd2Sn2O7 has a large $`\lambda `$-anomaly in SH and Mössbauer spectroscopy Bonville et al.  and muon spin relaxation ($`\mu `$SR) Yaouanc et al.  are consistent with the presence of long-range ordering. A $`\lambda `$-anomaly is detected at $``$ 1.2 K for Er2Ti2O7 Blöte et al. (1969) and long-range order is observed at low temperature Champion et al. (2003). For Yb2Ti2O7 a $`\lambda `$-anomaly is found at $`T_\lambda 0.21`$ K Blöte et al. (1969) but there are no long-range magnetic correlations below $`T_\lambda `$ but rather dynamical hysteretic short-range correlations Hodges et al. (2002). We note that beside the $`\lambda `$ peak, a broad peak centered near 2 K is present in SH. In the popular spin-ice systems Ho2Ti2O7 and Dy2Ti2O7, only a broad anomaly is present in SH Ramirez et al. (1999) and at low temperature the magnetic moment are frozen with no long-range order Harris et al. (1997); Bramwell and Gingras (2001).

Concerning the garnets, only a broad SH peak is detected at low temperature for Gd3Ga5O12 Schiffer et al. (1994). This compound does not display any long-range order Petrenko et al. (1998). Dy3Ga5O12 displays an SH anomaly at $`T_\lambda 0.37`$ K, overlapping with a broad peak Filippi et al. (1980a). Neutron scattering results show the presence of a long-range magnetic order below $`T_\lambda `$ Filippi (1981). Finally, while a $`\lambda `$-anomaly is found for Yb3Ga5O12 at $`T_\lambda `$ = 54 mK Filippi et al. (1980b), no information is available in the literature about the existence of a short or long-range magnetic order. We note that a relatively large SH hump is also present, centered at $``$ 0.2 K.

Since no $`\lambda `$-anomaly is present in the spin-ice systems and in Gd3Ga5O12 we shall not consider these systems any longer. To gain further understanding of the relationship between entropy and magnetic correlations in Yb3Ga5O12, we performed $`\mu `$SR measurements on this system. Here we assume, that the measured SH corresponds to magnetic dipole moment degrees of freedom and not to another exotic order parameter Paixão et al. (2002).

In the garnet lattice (space group $`Ia\overline{3}d`$), the crystal field acts on the $`{}_{}{}^{2}\mathrm{F}_{7/2}^{}`$ state of a Yb3+ ion to leave a well isolated ground state Kramers doublet and three closely grouped excited doublets having an average energy corresponding to $`\mathrm{\Delta }_{\mathrm{ave}}/k_\mathrm{B}`$ 850 K Buchanan et al. (1967). In a good approximation, the ground state of Yb3+ ions diluted in Y3Ga5O12 is described by an effective spin $`\mathrm{S}^{}`$ = 1/2 with an isotropic $`g`$-factor, $`g=3.43`$ Carson and White (1960); Wolf et al. (1962). The same description of Yb3+ in Yb3Ga5O12 is expected to be valid. This is supported by the electronic entropy variation below 1 K which is very close to $`R\mathrm{ln}2`$ Filippi et al. (1980b), as expected for a doublet. $`R`$ is the gas constant. From the measured $`g`$-factor and bulk magnetization Filippi et al. (1980b), the mean value of the Yb3+ magnetic moment at low temperature is found to be $`1.7\mu _\mathrm{B}`$. The interionic interaction is essentially antiferromagnetic as shown from the susceptibility above $`T_\lambda `$ which is lower than that extrapolated from high temperature for a Curie-Weiss behaviour Filippi et al. (1980b). From the position of the maximum of the bump of the SH an estimate for this interaction strength corresponds to at least 0.2 K. The antiferromagnetic nature of the interactions together with the crystallographic structure lead to frustration.

A sample was prepared by heating the constituent oxides to 1100 C four times with intermediate grindings. Its quality was checked by x-ray diffraction and susceptibility measurements. Zero-field $`\mu `$SR measurements were performed at the ISIS and PSI muon facilities covering the temperature range from 21 mK to 290 K. Additional spectra were recorded with a longitudinal field.

170Yb Mössbauer spectroscopy experiments Hodges et al. (2003) performed between 36 mK and 4.2 K show no resolved hyperfine structure at any temperature including below $`T_\lambda `$ as depicted in the inset of Fig. 1. This means that the moments fluctuate at a frequency larger than $``$ 300 $`\mu `$s-1 down to 36 mK. This also indicates the absence of magnetic correlations either short- or long-range.

We now present zero-field $`\mu `$SR results. The technique consists of implanting polarized (along direction $`Z`$) muons into a specimen and monitoring $`P_Z^{\mathrm{exp}}(t)`$ which is the evolution of the muon ensemble polarization projected onto direction $`Z`$ Dalmas de Réotier and Yaouanc (1997). The quantity actually measured is the so-called asymmetry corresponding to $`a_0P_Z^{\mathrm{exp}}(t)`$, where $`a_0`$ $``$ 0.24

For a paramagnet, $`P_Z^{\mathrm{exp}}(t)`$ tracks the dynamics of the magnetic field at the muon site, $`𝐁_{\mathrm{loc}}`$, reflecting the dynamics of the electronic moments. In the fast fluctuation or motional narrowing limit, $`P_Z^{\mathrm{exp}}(t)`$ takes an exponential form characterized by a relaxation rate, $`\lambda _Z`$. A stochastic approach to relaxation in zero-field leads to the relation $`\lambda _Z=2\gamma _\mu ^2\mathrm{\Delta }_{\mathrm{ZF}}^2\tau _c`$, assuming a single correlation time, $`\tau _c`$, for the magnetic moment dynamics. $`\mathrm{\Delta }_{\mathrm{ZF}}`$ is the standard deviation of the components of $`𝐁_{\mathrm{loc}}`$ and $`\gamma _\mu `$ the muon gyromagnetic ratio ($`\gamma _\mu `$ = 851.615 Mrad s-1 T-1). The motional narrowing limit is valid if $`\gamma _\mu \mathrm{\Delta }_{\mathrm{ZF}}\tau _c1`$.

Two examples of $`\mu `$SR spectra are shown in Fig. 1. All the spectra were satisfactorily fitted to an exponential function. The temperature dependence of $`\lambda _Z`$ is shown in Fig. 2. Since there is no qualitative change in $`P_Z^{\mathrm{exp}}(t)`$ and since $`\lambda _Z`$ increases monotonously as the temperature is lowered, there is neither long-range nor short-range magnetic order in Yb3Ga5O12. This means that the specific heat anomaly at $`T_\lambda `$ does not correspond to the onset of a conventional phase transition, in agreement with the 170Yb Mössbauer results. This is our first result.

The fact that we observe an exponential relaxation function below $`T_\lambda `$ is another important point. It implies that we are in the fast fluctuation limit (we shall check this quantitatively below) and therefore the Yb3+ moments continue to fluctuate rapidly down to the lowest temperature investigated. This behaviour is very different from that in Yb2Ti2O7 where the moments abruptly slow down below $`T_\lambda `$ where they are quasi-static Hodges et al. (2002).

We now comment and interpret the shape of $`\lambda _Z(T)`$.

As for nuclear magnetic resonance Moriya (1962), $`\lambda _Z`$ can be expressed in terms of the static wavevector dependent susceptibility, $`\chi (𝐐)`$, and linewidth of the quasielastic peak of the imaginary component of the dynamical susceptibility, $`\mathrm{\Gamma }(𝐐)`$. These two functions are assumed to be scalar, consistent with the quasi-isotropic nature of the Kramers ground-state doublet. Following Ref. Yaouanc et al. (1993) we write

$`\lambda _Z={\displaystyle \frac{\mu _0\gamma _\mu ^2}{16\pi ^2}}k_\mathrm{B}T{\displaystyle 𝒞(𝐐)\frac{\chi (𝐐)}{\mathrm{\Gamma }(𝐐)}\frac{\mathrm{d}^3𝐐}{(2\pi )^3}}.`$ (1)
$`\mu _0`$ is the permeability of free space and $`k_\mathrm{B}`$ the Boltzmann constant. The term $`𝒞(𝐐)`$ accounts for the interaction of the muon spin with the Yb3+ effective spins. The integration is over the first Brillouin zone.

For a Heisenberg magnet at high temperature, i.e. if the thermal energy is larger than the exchange energy, the $`𝐐`$ dependence of $`\chi (𝐐)`$ becomes negligible (see e.g. Ref. Lovesey (1984)) and it is simply given by the Curie law. Since $`\mathrm{\Gamma }(𝐐)`$ is temperature independent in the same limit Lovesey (1984), we deduce that $`\lambda _Z`$ should be temperature independent. This is effectively observed for $`0.4T80`$ K.

However, above 80 K, $`\lambda _Z(T)`$ decreases steadily as the sample is heated. Such a behaviour is sometimes associated with muon diffusion, but it would be quite unusual in an insulating oxide. In fact the decrease of $`\lambda _Z`$ is due to the relaxation of the Yb3+ magnetic moments resulting from an Orbach process, i.e. a two-phonon real process with an excited crystal-field level as intermediate Orbach (1961); Orb . A global fit is achieved with the formula $`\lambda _Z^1=A+B_{\mathrm{me}}\mathrm{exp}\left[\mathrm{\Delta }_e/\left(k_\mathrm{B}T\right)\right]`$. $`\mathrm{\Delta }_e`$ is the energy of the excited crystal-field level involved. $`A\lambda _Z^1`$ when $`\lambda _Z`$ saturates at low temperature (but still $`T>0.4`$ K). The constant $`B_{\mathrm{me}}`$ models the magneto-elastic coupling of the Yb3+ spin with the phonon bath. The result of the fit for $`0.4T290`$ K is shown in Fig. 2. Taking $`\mathrm{\Delta }_e/k_\mathrm{B}`$ = $`\mathrm{\Delta }_{\mathrm{ave}}/k_\mathrm{B}`$ = 850 K, we get $`A`$ = 0.69 (2) $`\mu `$s and $`B_{\mathrm{me}}`$ = 250 (70) $`\mu `$s. Combining the value of $`\lambda _Z`$ in the range 0.4 $`<T<`$ 80 K and the value $`\tau _c`$ 38 ps obtained from perturbed angular correlation and Mössbauer data in the same temperature interval Hodges et al. (2003), we compute $`\mathrm{\Delta }_{\mathrm{ZF}}`$ 0.16 T which is of the expected magnitude. For comparison $`\mathrm{\Delta }_{\mathrm{ZF}}`$ 0.08 T was found for Yb2Ti2O7 above $`T_\lambda `$ Yaouanc et al. (2003).

Information relevant to the influence of the frustrated nature of the magnetic interactions in Yb3Ga5O12 is obtained from $`\lambda _Z(T)`$ at low temperatures. As shown in Fig. 2, $`\lambda _Z`$ starts to deviate at $``$ 0.4 K from the behaviour expected for the Heisenberg Hamiltonian in the high-temperature limit. As we cool down the sample, we first observe a mild linear increase of $`\lambda _Z`$, with slope $`2.6(3)`$ $`\mu `$s-1K-1, followed by a sharp increase below $`T_\lambda `$ with slope $`21(7)`$ $`\mu `$s-1K-1.

In the temperature range $`T_\lambda <T<0.4`$ K, we note that the broad SH peak centered at $``$ 0.2 K covers the temperature régime where $`\lambda _Z(T)`$ varies slowly, suggesting that both features reflect the same physics. Since there are no excited crystal-field energy levels at low energy, the hump can not be of a Schottky-type. This extra SH reflects short-range correlations among groups of spins, see e.g. Ref. Morrish (1965). This means that the spin dynamics is wavevector dependent below $``$ 0.4 K. This dependence is expected to generate a mild increase of $`\lambda _Z`$ as shown by the explicit computation of Paja and coworkers Paja et al. (2002) for a simple model.

Now we turn our attention to the region $`T<T_\lambda `$. Being in an interstitial site, the muon spin can be strongly influenced by magnetic pair-correlations. In contrast, since a 170Yb Mössbauer nucleus is embedded in a Yb atom which is magnetic, it is mainly sensitive to the self correlations. Taking into account that the dynamics measured by Mössbauer spectroscopy is not wildly different above and below $`T_\lambda `$ Hodges et al. (2003), we infer that the sharp increase of $`\lambda _Z`$ occuring right below $`T_\lambda `$ is the signature of the building up of magnetic pair-correlations. Down to the lowest temperature, the $`\mu `$SR spectra have been recorded in the motional narrowing limit, since we compute $`\gamma _\mu \mathrm{\Delta }_{\mathrm{ZF}}\tau _c0.021`$ at low temperature, taking as an estimate for $`\tau _c`$ the hyperfine field correlation time from Mössbauer ($`0.3`$ ns) and computing $`\mathrm{\Delta }_{\mathrm{ZF}}`$ $``$ 0.08 T from the relation $`\lambda _Z=2\gamma _\mu ^2\mathrm{\Delta }_{\mathrm{ZF}}^2\tau _c`$. Therefore, the exponential form of the relaxation is fully justified.

An experimental estimate for $`\tau _c`$ can in principle be obtained from the analysis of the field dependence of $`\lambda _Z`$ assuming that the properties of the system are not modified by an external field. In the case of Yb3Ga5O12 we have found that a field as low as 0.3 T has a strong influence, increasing rather than quenching $`\lambda _Z`$. We conclude that the field influences the system: this is not completely astonishing recalling that a field of 0.7 T induces a phase transition for the isomorphous compound Gd3Ga5O12 Schiffer et al. (1994).

Yb3Ga5O12 is therefore a second system with frustrated magnetic interactions, after Yb2Ti2O7, where no long-range order is found below $`T_\lambda `$. In the pyrochlore a sharp first-order transition appears in the fluctuation rate of the magnetic moments and below $`T_\lambda `$ they continue to fluctuate slowly (in the megahertz range) at a temperature independent frequency. In the garnet we have a new scenario: the fluctuations are still rapid below $`T_\lambda `$ and no abrupt change in their frequency is observed. Only dynamical short-range correlations build up. The continuous rise of $`\lambda _Z`$ is magnetic in origin, and therefore excludes a dimer phase to be formed below $`T_\lambda `$.

Among the different frustrated pyrochlores and garnets with a $`\lambda `$ peak in SH that we have considered at the beginning of this letter, i.e. Gd2Ti2O7, Gd2Sn2O7, Er2Ti2O7, Yb2Ti2O7, Dy3Ga5O12 and Yb3Ga5O12, some of them order and others do not. One can tentatively find a condition for a magnetic order to be present. For all these compounds the variation of magnetic entropy at low temperature, say below $``$ 10 K, is close to the expected value deduced from the number of electronic degrees of freedom. This is $`R\mathrm{ln}8`$ for the Gd based compounds and $`R\mathrm{ln}2`$ for the others. Inspecting Refs. Ramirez et al. (2002); Bonville et al. ; Blöte et al. (1969); Hodges et al. (2002); Filippi et al. (1980b, a); Filippi (1981), the fraction of the total entropy frozen at $`T_\lambda `$ is roughly 60 % for Er2Ti2O7 and Dy3Ga5O12, 40 % for Gd2Sn2O7 and 35 % for Gd2Ti2O7 whereas it is only $``$ 20 % for Yb2Ti2O7 and Yb3Ga5O12. Therefore we estimate that a long-range magnetic order is only present when at least $`1/41/3`$ of the magnetic entropy is released at $`T_\lambda `$.

While the entropy change at the $`\lambda `$-anomaly temperature is about the same for both Yb3Ga5O12 and Yb2Ti2O7, the Kramers ground-state doublet is approximately isotropic for the former and strongly anisotropic for the latter. This suggests that a high magnetic anisotropy impedes the system from exploring different configurations, leading to an abrupt slowing down of the correlations. $`\mathrm{\Delta }_{\mathrm{ZF}}`$ is reduced by a factor $``$ 14 Hodges et al. (2002); Yaouanc et al. (2003) when crossing $`T_\lambda `$ from above for Yb2Ti2O7 whereas this reduction is much less ($``$ 2) for Yb3Ga5O12. This is again consistent with the anisotropy difference between both systems: we expect $`\mathrm{\Delta }_{\mathrm{ZF}}`$ to be small if the magnetic moments are confined to specific orientations.

In conclusion, according to wisdom, a $`\lambda `$-type anomaly is indicative of a second order phase transition at $`T_\lambda `$ for a three-dimensional system. However, due to the presence of frustrated interactions, only dynamical short-range correlations build up below $`T_\lambda `$ in Yb3Ga5O12. By comparison with results from other three-dimensional frustrated systems, we infer that the long-range order of the order parameter does not occur if a too large fraction of the entropy is released at a broad bump located above $`T_\lambda `$. The anisotropy of the system favors a freezing of the magnetic correlations.

We are grateful to J.A. Hodges for initiating this project and for his constant interest. We thank A. Forget for preparing the sample, P. Bonville and J.A. Hodges for communication of their Mössbauer spectroscopy results prior to publication and useful discussions and M.E. Zhitomirsky for a careful reading of the manuscript. The authors from the Netherlands are grateful to the Dutch Scientific Organisation (NWO) for its financial support for the use of ISIS. Support from the European Community through its Access to Research Infrastructure action of the Improving Human Potential Program is acknowledged.


text:
# On the Bottom Magnetic Fields of the Millisecond Pulsars

## References

Brown, F.B., & Bildsten, L. 1998, ApJ, 496, 915

Cheng, K.S.,& Zhang, C.M. 1998, A&A, 337, 441; A&A, 361, 1001

Payne, R., & Melatos, A. 2003, in this Proceedings.

Romani, G.W. 1990, Nature, 347, 741

Shapiro, S.L., & Teukolsky, S.A. 1983, Black Holes, White Dwarfs and Neutron Stars (Wiley, New York)

Shibazaki, N., & Murakami, T., & Shaham, J., & Nomoto, K. 1989, Nature, 342, 656

Taam, R.E., & Van den Heuvel, E.P.J. 1986, ApJ, 305, 235

Urpin, V., & Geppert, U. 1995, MNRAS, 275, 1117

van den Heuvel, E.P.J., & Bizaraki, O. 1995, A&A, 297, L41

White, N., & Zhang, W. 1997, ApJ, 490, L87


text:
# Atmosphere-Corrected Phase-Referencing

## 1. Introduction

Phase-referencing has become a standard technique in Very Long Baseline Interferometry (VLBI) to measure accurate relative positions of radio sources. The concept of phase-referencing is based on the assumption that the phase errors of two sources with a small angular separation on the sky are similar. One observes a target source between two scans on a calibrator. Then one can interpolate the calibration between the two scans and apply this interpolated phase corrections to the phase of the target source. The phase-referenced difference phase contains information only about the target source structure and its position relative to the calibrator and contains noise from interpolation errors.

One major problem of phase-referencing VLBI observations is the unknown tropospheric zenith delay at each antenna. The VLBA correlator model uses a seasonally averaged and latitude-dependent atmospheric model (Niell 1996), which can misestimate the zenith delay by a few centimeters. The excess path length caused by the dry troposphere to a radio signal from a source with zenith angle $`Z`$ is given by

$`l0.228{\displaystyle \frac{\mathrm{cm}}{\mathrm{mb}}}P_0\mathrm{sec}Z(10.0013\mathrm{tan}^2Z)`$ (1)
where $`P_0`$ is the total pressure at the surface (see equation 13.41 in Thompson, Moran, & Swenson 2001). This excess path length will be different for the calibrator and the target source since both sources usually have different zenith angles. Hence, a residual error will remain after phase-referencing. These residual errors degrade the quality of the phase-referenced image and limit the achievable astrometric accuracy (e.g., Beasley & Conway 1995). If the correlator model misestimates the zenith delay $`l_0`$ by 3 cm ($`0.1`$ ns), one would expect after phase-referencing delay errors of 0.2 cm for zenith angles of $`60^{}`$ and source separations of $`1^{}`$. However, these errors can be corrected if the true zenith delay at each antenna is known.

## 2. Atmosphere-Corrected Phase-Referencing

### 2.1. Phase fitting

The fringe-phase of the target source after phase-referencing is mainly a sum of the position offset and the effects of the zenith delay error if the source structure is negligible. The phase induced by the position offset shows a simple 24 hour cosine behavior, while the phase errors from the zenith delay error follow a more complex dependence on the zenith angle according to equation 1. Because of the different behavior of the two contributions, it is possible to separate both effects and to estimate the position offset as well as the zenith delay error.

One can compute a model phase with a position offset (two parameters) and a zenith delay error at each antenna (one parameter for each antenna) and perform a least-squares fit to the phase data. This method has been applied successfully by Reid et al. (1999) to measure the proper motion of Sgr A. In Brunthaler et al. (2002) we applied this method to phase-referencing observations of extragalactic H2O masers in IC 10 with respect to the background quasar J0027+5958 and achieved an astrometric accuracy of $`10\mu `$as.

Second epoch observations of the H2O masers in IC 10 were made in January 2002 and two observations were made within five days to check the accuracy and repeatability of our results. The data were calibrated with standard techniques using the AIPS software package. A priori amplitude calibration was applied using system temperature measurements and standard gain curves. A fringe fit was performed on J0027+5958 and the solutions were applied to IC 10.

Then we modeled the phase data of IC 10 as described above. The position offset for the two observations are given in Table 1 (A). The positions obtained from the two observations are consistent and have an rms of $`10\mu `$as. The zenith delay offsets are typically of the order of a few centimeters. The values of the first observation on 2002 January 12 are given in the second column of Table 2.

### 2.2. Phase correction

The fitted zenith delay offsets can be used to correct the uv-data itself. This correction was done using a new version of the AIPS task CLCOR with OPCODE=’ATMO’. Fig. 1 shows the phase-referenced images of a strong H2O maser component in IC 10 with (b) and without (a) the corrections. The noise in the image is reduced while the peak flux increased in the corrected image. The peak-to-noise ratio in the corrected image is 101 compared to 77 in the uncorrected image. The positions of the maser feature in the corrected maps for the two observations are given in Table 1 (B) and agree within 10 $`\mu `$as with the positions derived from the fit to the phases. The positions obtained from the fit to the phases and from the corrected maps are not independent and the agreement is not surprising. However, it is a consistency check for the method.

### 2.3. Geodetic-like observations

A second, independent method to estimate the zenith delay offset is to use geodetic-like observations. Here, a number of bright quasars with positions known to better than 1 mas were observed for 45 minutes before and after the phase-referencing observation. The quasars were observed at different elevations with a frequency setup of 8 IFs spread over 450 MHz. A fringe-fit on the quasars yielded then a multi-band delay and a rate for each baseline and source. These multi-band delays and rates were then fit with a model that consisted of a zenith delay offset and drift at all antennas as well as a clock offset and drift at all antennas except the reference antenna. The parameters of the best model-fit are given in Table 2. The differences between the observed multi-band delay and the model are less than 0.1 ns.

The zenith delay offsets obtained from the fit to the phases and from the geodetic-like observations are in general comparable although they can be different for some antennas. The differences between the two methods are not surprising, since the model in the geodetic-like observations includes not only a zenith delay offset but also a drift, a clock offset and a clock drift.

The fit values of the zenith delay and clock offsets can also be used to correct the uv-data. A corrected image is seen in Figure 1 (c) and the peak-to-noise ratio is 100. The positions of the maser feature in the corrected maps are given in Table 1 (C). They are also consistent with the positions of the fit to the phases.

## 3. Conclusions

We have presented a technique to correct phase-referencing observations for inaccurate zenith delays in the VLBA correlator model. The quality of the atmosphere-corrected phase-referenced images can be improved by this technique. The peak-to-noise ratio in the corrected images increases by $`30\%`$. By using this technique, an astrometric accuracy of 10 $`\mu `$as can be achieved.

Two independent methods were presented to estimate the zenith delay offsets at each antenna. The images obtained from the correction with the zenith delay offsets from the fit to the phases (Fig.1 (b)) and from the geodetic-like observations (Fig.1 (c)) are very similar and the positions measured from these images are consistent. This is a proof of concept for our technique.

The astrometric accuracy of this technique enables us to detect significant proper motions of galaxies in the Local Group out to a distance of 800 kpc within a few years. One needs at least three epochs of observations to make reliable statements about measured proper motions. Further observations are underway and the results will be presented elsewhere.

## References

Beasley, A., & Conway, J. E. 1995, VLBI Phase-Referencing, in ASP Conf. Ser. 82: Very Long Baseline Interferometry and the VLBA, eds. J. A. Zensus, P. J. Diamond & P. J. Napier (San Francisco: ASP), 327

Brunthaler, A., Reid, M. J., Falcke, H., Greenhill, L. J., & Henkel, C. 2002 in Prodeedings of the 6th European VLBI Network Symposium, E. Ros, R. W. Porcas, A. P. Lobanov, and J. A. Zensus (eds.), MPIfR Bonn, Germany, p. 189.

Niell, A. 1996, J. Geophys. Res., 101, 3227

Reid, M. J., Readhead, A. C. S., Vermeulen, R. C., & Treuhaft, R. N. 1999, ApJ, 524, 816

Thompson, A. R., Moran, J. M.,& Swenson, G. W. 2001, Interferometry and Synthesis in Radio Astronomy (Second ed.). Wiley & Sons, New York


text:
# Electronic properties of structural twin and antiphase boundaries in materials with strong electron-lattice couplings

## Abstract

Using a symmetry-based atomic scale theory of lattice distortions, we demonstrate that elastic textures, such as twin and antiphase boundaries, can generate intricate electronic inhomogeneities in materials with strong electron-lattice coupling, as observed in perovskite manganites and other functional materials.

Recent advances in imaging techniques have revealed the presence of rich elastic textures in functional materials such as colossal magnetoresistive (CMR) manganites,Uehara99 ; Fukumoto99  ferroelectrics, Abakumov00  ferroelastics,Cao01  and shape memory alloys.Ren02  In particular, recent experiments on certain perovskite manganites compoundsFukumoto99  have shown the correlation between electronic transport properties and the presence of meandering antiphase boundaries (APBs) within insulating charge ordered domains, interpreted as the existence of metallic regions forming around APBs. It is also reported that strains near grain boundaries in thin film can modify electronic properties in manganites.Soh00  The interplay between elastic texture and electronic heterogeneity is thus central to understanding multiphase coexistence and the resultant electronic properties in CMR and other functional materials.

In this Report, we illustrate the importance of elastic inhomogeneities in modifying electronic properties in materials with strong electron-lattice coupling. In particular, we study the electronic properties of APBs and twin boundaries (TBs) on a two-dimensional (2D) lattice. We first show how a recently-developed symmetry-based atomic-scale theory of lattice distortions can be used to find atomic configurations of twin and antiphase boundaries. Within our framework, we illustrate the differences and similarities between TBs and APBs from the point of view of localization of long \[short\] wavelength modes inside APBs \[TBs\], evolution upon energy relaixtion and roughness TB.v.APB  \[smoothness\] of APBs \[TBs\]. We subsequently perform a tight binding calculation with a Su-Schrieffer-Heeger (SSH) type model of electron-lattice coupling to predict the distribution of electronic density of states, which can be related to the results of scanning tunneling microscope (STM) measurements. Our work thus forms the basis for predicting electronic properties from predesigned materials microstructures.

Twin boundaries separate domains related by the rotation of crystalline axes, whereas APBs represent boundaries at which the sequence of alternating distortions, such as alternating rotational directions of oxygen octahedra in perovskite oxides, change their phase (i.e., broken translational symmetry). Although our method can be applied to 2D or 3D lattices with mono- or multiatomic bases, we illustrate our ideas with a square lattice in 2D space with a monatomic basis, for which the appropriate atomic scale distortion variables are the modes shown in Fig. 1 (Ref. Ahn03, ). These modes have important advantages over displacement variables – they reflect the symmetries of the lattice and can serve as order parameters (OP) in structural phase transitions. Therefore, energy expressions with desired ground states can be written in a simpler way in these variables than in usual displacement variables. Moreover, since the lattice distortions are decomposed into the modes at $`\stackrel{}{k}=(0,0)`$ (long-wavelength or intercell modes) and $`\stackrel{}{k}=(\pi ,\pi )`$ (short-wavelength or intracell modes), the approach using these modes shows the differences between long and short wavelength lattice distortions in a natural way.

We consider APBs, such as the one shown in Fig. 2 where open circles represent the distorted atomic positions, for $`s_x`$ or $`s_y`$ modes. The simplest energy expression yielding a ground state with either pure $`s_x`$ or $`s_y`$ mode lattice distortion is

$`E_{\text{sxsy}}`$ $`=`$ $`{\displaystyle \underset{\stackrel{}{i}}{}}{\displaystyle \frac{D}{2}}[s_x(\stackrel{}{i})^2+s_y(\stackrel{}{i})^2]+{\displaystyle \frac{G_1}{4}}[s_x(\stackrel{}{i})^4+s_y(\stackrel{}{i})^4]`$ (1)
$`+`$ $`{\displaystyle \frac{G_2}{2}}s_x(\stackrel{}{i})^2s_y(\stackrel{}{i})^2+{\displaystyle \underset{\stackrel{}{i},n=1,2,3}{}}{\displaystyle \frac{C_n}{2}}e_n(\stackrel{}{i})^2,`$
where the coefficients $`D`$, $`G_1`$, $`G_2`$ (with $`G_1<G_2`$) can be obtained from interatomic forces and $`C_n`$ are the associated elastic moduli.

Since the ground state has $`\stackrel{}{k}=(\pi ,\pi )`$ component lattice distortions, we define variables with tilde by multiplying $`(1)^{i_x+i_y}`$ with all intercell/intracell/displacement variables \[e.g., $`\stackrel{~}{s}_x(\stackrel{}{i})`$=$`s_x(\stackrel{}{i})(1)^{i_x+i_y}`$, $`\stackrel{~}{d}_x(\stackrel{}{i})`$=$`d_x(\stackrel{}{i})(1)^{i_x+i_y}`$\]. In wavevector space, this corresponds to interchanging the $`(0,0)`$ and $`(\pi ,\pi )`$ points. Therefore, $`\stackrel{~}{s}_x`$ and $`\stackrel{~}{s}_y`$ are the modes near $`\stackrel{}{q}=(0,0)`$, and $`\stackrel{~}{e}_1`$, $`\stackrel{~}{e}_2`$, and $`\stackrel{~}{e}_3`$ are the modes near $`\stackrel{}{q}=(\pi ,\pi )`$, where $`\stackrel{}{q}`$ represents the wavevector for the variables with tilde. The constraint equations Ahn03  in the continuum limit are $`\stackrel{~}{e}_{1,3}(\stackrel{}{r})=(_y\stackrel{~}{s}_x\pm _x\stackrel{~}{s}_y)/(2\sqrt{2})`$ and $`\stackrel{~}{e}_2(\stackrel{}{r})=(_x\stackrel{~}{s}_x+_y\stackrel{~}{s}_y)/(2\sqrt{2})`$. These relations show that spatial variations of intracell modes ($`\stackrel{~}{s}_x`$, $`\stackrel{~}{s}_y`$) always generate intercell modes ($`\stackrel{~}{e}_1`$, $`\stackrel{~}{e}_2`$, $`\stackrel{~}{e}_3`$), which becomes microscopic origin of the phenomenological gradient energy term in continuum elasticity theory. Therefore, intercell modes are present inside an APB at which an intracell mode with tilde changes its sign. Similar arguments apply to a TB Ahn03  – the spatial variations of intercell modes always generate intracell modes inside TB through their constraints – and lead to a similar solitary profile for TB and APB. Barsch84 ; Cao01  However, there exists a fundamental difference between intercell and intracell modes: In the $`k0`$ limit, $`e_1`$, $`e_2`$, and $`e_3`$ are given as a first derivative of $`d_x`$ and $`d_y`$, whereas in the $`q0`$ limit, $`\stackrel{~}{s}_x`$ and $`\stackrel{~}{s}_y`$ are $`2\stackrel{~}{d}_x`$ and $`2\stackrel{~}{d}_y`$. Therefore, $`e_1`$, $`e_2`$, and $`e_3`$ are related by constraint equations, whereas $`\stackrel{~}{s}_x`$ and $`\stackrel{~}{s}_y`$ are not constrained by each other.

We consider the APB between the two ground states with $`\stackrel{~}{s}_x(\stackrel{}{i})=\pm \sqrt{D/G_1}`$ (Fig. 2). On expanding Eq. (1) around the ground state and retaining leading order terms, we write the energy $`E_{\text{APB}}=E_{\text{APB,OP}}+E_{\text{APB,NOP}}`$, where

$`E_{\text{APB,NOP}}`$ $`=`$ $`{\displaystyle \underset{\stackrel{}{i},n=1,2,3}{}}{\displaystyle \frac{C_n}{2}}\stackrel{~}{e}_n(\stackrel{}{i})^2+{\displaystyle \underset{\stackrel{}{i}}{}}{\displaystyle \frac{D^{}}{2}}\stackrel{~}{s}_y(\stackrel{}{i})^2`$ (2)
$`E_{\text{APB,OP}}`$ $`=`$ $`{\displaystyle \underset{\stackrel{}{i}}{}}{\displaystyle \frac{D}{2}}\stackrel{~}{s}_x(\stackrel{}{i})^2+{\displaystyle \frac{G_1}{4}}\stackrel{~}{s}_x(\stackrel{}{i})^4,`$ (3)
and $`D^{}=D(G_2/G_11)`$. The energy $`E_{\text{APB}}`$ has the identical form to that for a TB with $`e_3`$ as OP Ahn03  – a double-well potential for OP and harmonic potentials for non-OP. For both TB and APB problems, the non-OP energy terms mediate an anisotropic interaction between OPs through the constraints. The minimization of $`E_{\text{APB,NOP}}`$ using Lagrange multipliers leads to $`E_{\text{APB,NOP}}^{\text{min}}=_\stackrel{}{q}\frac{1}{2}\stackrel{~}{s}_x(\stackrel{}{q})\stackrel{~}{U}(\stackrel{}{q})\stackrel{~}{s}_x(\stackrel{}{q})`$, where $`\stackrel{~}{U}(\stackrel{}{q})`$ has a $`q^2`$ leading order term with an anisotropic coefficient $`\stackrel{~}{U}_2(\theta _q)=[(C_1+C_2+C_3)+(C_2C_1C_3)\mathrm{cos}2\theta _q]/16`$. When transformed into real space, the $`q^2`$ leading order term gives rise to a short range $`R^4`$ interaction between OP $`\stackrel{~}{s}_x`$, where $`R`$ is the distance between two sites, unlike the long-range $`R^2`$ interaction between OP $`e_3`$ for the TB case.Shenoy99 ; LongRange  Such different range of interaction is consistent with well-known understandingTB.v.APB  of smoothness of TB and roughness of APB, which our simulation reproduces below.

The physical origin of the long-range interaction between intercell modes and short-range interaction between intracell modes lies in the difference in the symmetry operations relating the two domains separated by a TB or APB. For a TB, the sign change in OP corresponds to a change in orientation, which has no intrinsic length scale and thus gives rise to a long range interaction. For an APB, the sign change in OP signifies translation of the configuration by one atomic spacing, implying the presence of an intrinsic atomic length scale that is responsible for the fast decay of the interaction between OPs. These considerations apply to any 2D/3D lattice with mono- or multiatomic bases. We can also expect it to be valid for the rough and fluctuating stripes suggested in high-$`T_C`$ cupurates, which are examples of magnetic APBs. Arguments similar to those in Ref. Ma.SM, , related to the range of the interactions, explain the characteristic smoothness and roughness of TBs and APBs.

The full expression of the kernel $`\stackrel{~}{U}(\stackrel{}{q})`$ in $`E_{\text{APB,NOP}}^{\text{min}}`$ is given by

$`\stackrel{~}{U}(\stackrel{}{q})`$ $`=`$ $`{\displaystyle \frac{C_1+C_3}{2}}T_y(\stackrel{}{q})^2+{\displaystyle \frac{C_2}{2}}T_x(\stackrel{}{q})^2`$
$``$ $`{\displaystyle \frac{(C_1+C_2C_3)^2T_x(\stackrel{}{q})^2T_y(\stackrel{}{q})^2}{4D^{}+2(C_1+C_3)T_x(\stackrel{}{q})^2+2C_2T_y(\stackrel{}{q})^2}},`$
where $`T_x(\stackrel{}{q})=\mathrm{tan}(q_x/2)`$ and $`T_y(\stackrel{}{q})=\mathrm{tan}(q_y/2)`$. We choose similar parameter values for both APB and TB cases and use the Euler method Shenoy99  to relax the lattice starting from random initial conditions to see difference in evolution between TB and APB. The results, $`e_3`$ for TB and $`\stackrel{~}{s}_x`$ for APB, on a $`64\times 64`$ lattice with periodic boundary conditions are shown in Figs. 3 and 4, respectively, where the main panels show real space distributions and insets the $`k`$ or $`q`$-space distributions. Figures 3(a) and 4(a) correspond to early stages of the relaxation, which show the characteristics of the OP distribution at high temperature($`T`$), whereas Figs. 3(b) and 4(b) reflect the OP distributions for late stages at low $`T`$. Even at the early stage, the presence of the long range correlation between intercell OP along $`45^o`$ and $`135^o`$ can be identified in the main panel in Fig. 3(a), reminiscent of tweed structures in martensitic materials. Shenoy99  The $`k`$-space distribution plotted in the inset of Fig. 3(a) shows strong preference of $`e_3(\stackrel{}{k})`$ along $`45^o`$ and $`135^o`$ orientations. Such long-range correlation and anisotropy are absent for the phase of the intracell mode distortion, as shown in Fig. 4(a). The late stage of the relaxation depicted in Figs. 3(b) and 4(b) show characteristic smooth TBs Ren02  and rough APBs Fukumoto99 . The TB is metastable and cannot be removed by further relaxation, unless large noise is applied. In contrast, the APB in Fig. 4(b), which is a ring due to the periodic boundary condition, shrinks and disappears upon further relaxation. Although the solitary-wave profile of the smooth APB along a certain direction is a metastable state,Cao01  the absence of a long range interaction between the phase of intracell mode distortions prevents the relaxation of random initial configuration from reaching such metastable states. This indicates that lattice defects or boundary conditions are necessary to reach the metastable solitary-wave APB configuration in materials, and influence the geometry of APBs.

To study functional electronic aspects associated with these TB and APB microstructures, we consider the modulation of electronic properties based on the following SSH electron-lattice coupling Hamiltonian,

$$H_{\text{SSH}}=\underset{\stackrel{}{i},a=x,y}{}t_0[1\alpha (d_{\stackrel{}{i}+\widehat{a}}^ad_\stackrel{}{i}^a)](c_\stackrel{}{i}^{}c_{\stackrel{}{i}+\widehat{a}}+c_{\stackrel{}{i}+\widehat{a}}^{}c_\stackrel{}{i}),$$
(5)
where $`c_\stackrel{}{i}^{}`$ is the creation operator for an electron at $`\stackrel{}{i}`$. We use $`d^x(\stackrel{}{i})`$ and $`d^y(\stackrel{}{i})`$ obtained from our atomistic model as inputs to the SSH Hamiltonian. For $`t_0=1`$, $`\alpha =1`$, and the TB and APB results shown in Figs. 3(b) and 4(b), we find all energy levels and eigenstates numerically, and calculate the local density of state (DOS) at each site and the distributions of local DOS and charge density for chosen Fermi energies ($`E_F`$), as shown in Figs. 56 and 7.

The blue curve in the inset of Fig. 5(a) represents the local DOS within the domain of Fig. 3(b), close to the bulk DOS for the homogeneous phase. The energy difference between the two peaks is proportional to $`\alpha |e_3|`$. Within the TB, $`|e_3|`$ is small and therefore the local DOS spectral weight moves toward $`E=0`$, as shown with a red curve. This local shift of the DOS weight can be measured using direct local probes, such as STM. Uehara99 ; Davis01  For $`E_F`$ shown in this inset, the local DOS is smaller inside the TB than inside the domain, and the main panel of Fig. 5(a) shows the real space distribution of local DOS at this $`E_F`$. The oscillation of local DOS within the domain is related to Friedel oscillations. The feedback from electron to lattice neglected here may generate a similar Friedel oscillation of the lattice. The charge density, which is the area in local DOS below $`E_F`$, is a constant (0.5) if $`E_F=0`$. If $`E_F<0`$, more electrons are depleted from the TB than the domain, as shown in Fig. 5(b). The charge density also has a Friedel oscillation, though not clearly visible in Fig. 5(b).

The inset of Fig. 6(a) shows the local DOS for the APB case. Within the domain, each site has a neighbor along $`x`$ axis that is closer in one direction and farther in the other (Fig. 2). Such a distortion pattern results in a V-shaped local DOS with zero DOS at $`E=0`$, plotted as a blue curve in the inset of Fig. 6(a) (the small finite value of local DOS at $`E=0`$ is due to the finite energy bin size). Within the APB, the average of the two bond lengths along $`x`$-axis has a $`(\pi ,\pi )`$ component. For example, along sites at the center of the $`90^o`$ APB shown in Fig. 2, the neighboring sites in the $`x`$-direction are alternately closer or farther apart. The red line in the inset of Fig. 6(a) represents the local DOS for the sites with two farther neighbors within the center of APB. This has a peak near $`E=0`$. The green line represents the local DOS at sites with two closer neighbors, right next along $`y`$ direction to the sites for the red curve. The difference between the red and green curves emphasizes the difference that occurs in density of states due to atomic-scale changes in microstructure as along an APB. Such atomistically sharp changes at structural interfaces have been seen in CMR manganites with STM (Renner et al. Uehara99 ). The local DOS at $`E_F=0`$ is zero within the domain and finite only around the APB with a $`(\pi ,\pi )`$ component modulation, as depicted in the main panel of Fig. 6(a). This means that the electronic states created by the APB dominate the low energy properties of the whole system, e.g., conductivity or specific heat. The charge density for $`E_F`$ slightly higher than zero is plotted in Fig. 6(b). Checkerboard-type modulations of electronic properties around APBs, reflecting the electronic states coupled to the $`(\pi ,\pi )`$ type lattice distortions, demonstrate that STM studies of TB and APB can be used to reveal underlying electronic properties of materials (usually hidden in a homogeneous phase), similar to the way single impurity atoms have been used to reveal superconducting properties in high TC cupurates. Results for different values of $`E_F`$ are displayed in Fig. 7. Local DOS modulations, particularly the checkerboard type modulations, tend to reach farther into the domain than the elastic texture itself. We also note that the checkerboard type modulations with a wave vector $`(\pi ,\pi )`$ are present at any $`E_F`$, in contrast to Friedel oscillations for TBs which change wavevectors with $`E_F`$.

Although the orbital states and type of electron-lattice coupling in CMR manganites are quite different from the simple model presented here, our results suggest that simple APBs from double-well type potentials, such as Eq.(3), may give rise to electronic heterogeneities, but not undistorted metallic regions with uniform charge densities near APBs. An energy landscape with a local minimum at undistorted state, such as the one considered in Ref. Ahn04, , would be more likely to nucleate metallic domains at APBs, and create percolating conducting paths in CMR manganites.

In summary, we have shown that in functional materials with strong electron lattice coupling, the electronic properties are modified near elastic textures such as TBs and APBs, which can be directly measured by STM. The results also show that the heterogeneities of electron local DOS are not just confined within TBs and APBs, but can propagate into domains in the form of Friedel oscillations for TBs and with the wave vector of short wave length lattice distortions for APBs.

This work was supported by the US DOE.


text:
# An Infrared Renormalization Group Limit Cycle in QCD

The successes of the renormalization group (RG) range from the explanation for universality in critical phenomena to the nonperturbative formulation of quantum field theories. The RG can be reduced to a set of differential equations that define a flow in the space of dimensionless coupling constants. Asymptotic scale invariance at long distances, as in critical phenomena, can be explained by RG flow to an infrared fixed point. Asymptotic scale invariance at short distances in quantum field theories can be explained by RG flow to an ultraviolet fixed point. A fixed point is the simplest topological feature that can be exhibited by a RG flow. As pointed out by Wilson in 1971, the next simplest possibility is a limit cycle, a 1-dimensional orbit that is closed under the RG flow . The RG flow carries the system around a complete orbit every time the ultraviolet cutoff $`\mathrm{\Lambda }`$ increases by some factor $`\lambda _0`$. One of the signatures of an RG limit cycle is discrete scale invariance: symmetry with respect to scaling by a factor $`\lambda _0^n`$ only for integer values of $`n`$. Simple discrete Hamiltonian systems with RG limit cycles have been constructed . An ultraviolet limit cycle has been discovered in a familiar 2-dimensional field theory . Ultraviolet limit cycles have also been found in effective field theories (EFT’s) for bosons with large scattering length  and for low-energy nucleons . In this talk, we point out that Quantum Chromodynamics (QCD) may have an infrared limit cycle in the 3-nucleon sector for critical values of the up and down quark masses.

In 1970, Efimov discovered that if the scattering length $`a`$ of nonrelativistic bosons is much larger than the range of the interaction, there are shallow 3-body bound states whose number increases logarithmically with $`|a|`$ . In the resonant limit $`a\pm \mathrm{}`$, there are infinitely-many arbitrarily shallow Efimov states. The ratio of successive binding energies rapidly approaches a universal constant $`\lambda _0^2=e^{2\pi /s_0}515`$, where $`s_01.00624`$ is a transcendental number. The Efimov effect was revisited by Bedaque, Hammer, and van Kolck (BHvK) within the EFT framework . Bosons with large scattering length $`a`$ can be described by a nonrelativistic field theory with a 2-body contact interaction $`(g_2/4)(\psi ^{}\psi )^2`$. In the 2-body sector, renormalization can be implemented by adjusting the coupling constant $`g_2(\mathrm{\Lambda })`$ as a function of the ultraviolet momentum cutoff $`\mathrm{\Lambda }`$ so that the scattering length is $`a`$. As $`\mathrm{\Lambda }\mathrm{}`$, the dimensionless coupling constant $`\mathrm{\Lambda }g_2(\mathrm{\Lambda })`$ flows to an ultraviolet fixed point: $`\mathrm{\Lambda }g_2(\mathrm{\Lambda })4\pi ^2`$. The renormalization of $`g_2`$ eliminates ultraviolet divergences from 3-body observables, but they still depend on the ultraviolet cutoff as periodic functions of $`\mathrm{ln}(\mathrm{\Lambda })`$ with period $`\pi /s_0`$. BHvK showed that the 3-body sector could be fully renormalized by including a 3-body contact interaction $`(g_3/36)(\psi ^{}\psi )^3`$ . As $`\mathrm{\Lambda }\mathrm{}`$, the dimensionless coupling constant $`\mathrm{\Lambda }^4g_3(\mathrm{\Lambda })`$ flows to an ultraviolet limit cycle: $`\mathrm{\Lambda }^4g_3(\mathrm{\Lambda })144\pi ^4H(\mathrm{\Lambda })`$ , where

$`H(\mathrm{\Lambda })`$ $`=`$ $`{\displaystyle \frac{\mathrm{cos}[s_0\mathrm{ln}(\mathrm{\Lambda }/\mathrm{\Lambda }_{})+\mathrm{arctan}(s_0)]}{\mathrm{cos}[s_0\mathrm{ln}(\mathrm{\Lambda }/\mathrm{\Lambda }_{})\mathrm{arctan}(s_0)]}}`$ (1)
and $`\mathrm{\Lambda }_{}`$ is a scaling-violation parameter introduced by dimensional transmutation. In the resonant limit $`a\pm \mathrm{}`$, there is an infrared limit cycle. The asymptotic discrete scaling symmetry with scaling factor $`\lambda _022.7`$ is manifest in the spectrum of Efimov states.

The spin-singlet and spin-triplet scattering lengths $`a_s`$ and $`a_t`$ for nucleons are both significantly larger than the range of the nuclear force. Efimov used this observation as the basis for a qualitative approach to the few-nucleon problem in which the nuclear forces are approximated by zero-range potentials with depths adjusted to reproduce the scattering lengths $`a_s`$ and $`a_t`$ . This approach gives a reasonably accurate prediction for the deuteron binding energy. In the 3-nucleon sector, the Efimov effect makes it necessary to impose a boundary condition on the wavefunction at short distances that can be fixed by using the triton binding energy $`B_t`$ as input. The resulting prediction for the neutron-deuteron scattering length is surprisingly accurate. Efimov’s program has been revisited by BHvK within the EFT framework . The EFT involves an isospin doublet $`N`$ of Pauli fields with two independent 2-body contact interactions: $`N^{}\sigma _iN^cN^c\sigma _iN`$ and $`N^{}\tau _kN^cN^c\tau _kN`$, where $`N^c=\sigma _2\tau _2N^{}`$. Renormalization in the 2-body sector requires the two coupling constants to be adjusted as functions of $`\mathrm{\Lambda }`$ to obtain the correct values of $`a_s`$ and $`a_t`$. Renormalization in the 3-body sector requires a 3-body contact interaction $`N^{}\sigma _iN^cN^c\sigma _jNN^{}\sigma _i\sigma _jN`$ with a coupling constant proportional to (1). Thus the renormalization involves an ultraviolet limit cycle. The scaling-violation parameter $`\mathrm{\Lambda }_{}`$ can be determined by using the triton binding energy $`B_t`$ as input.

The low-energy few-nucleon problem can also be described by an EFT that includes explicit pion fields. Such an EFT has been used to extrapolate nuclear forces to the chiral limit of QCD in which the pion is exactly massless . The extrapolation to larger values of $`m_\pi `$ predicts that $`a_t`$ diverges and the deuteron becomes unbound at a critical value in the range 170 MeV $`<m_\pi <`$ 210 MeV. It is also predicted that $`a_s`$ is likely to diverge and the spin-singlet deuteron become bound at some critical value of $`m_\pi `$ not much larger than 150 MeV. Both critical values are close to the physical value $`m_\pi =135`$ MeV.

Chiral extrapolations can also be calculated using the EFT without pions . The inputs required are the chiral extrapolations $`a_s(m_\pi )`$, $`a_t(m_\pi )`$, and $`B_t(m_\pi )`$, which can be calculated using an EFT with pions. As an illustration, we take the central values of the error bands for the inverse scattering lengths $`1/a_s(m_\pi )`$ and $`1/a_t(m_\pi )`$ from the chiral extrapolation in Ref. . Since the chiral extrapolation of the triton binding energy $`B_t(m_\pi )`$ has not yet been calculated and since $`\mathrm{\Lambda }_{}`$ should vary smoothly with $`m_\pi `$, we approximate it by its physical value $`\mathrm{\Lambda }_{}=189`$ MeV for $`m_\pi =138`$ MeV. In Fig. 1, we show the resulting 3-body spectrum in the triton channel as a function of $`m_\pi `$. Near $`m_\pi 175`$ MeV where the decay threshold comes closest to $`\kappa =0`$, an excited state of the triton appears. This excited state is a hint that the system is very close to an infrared limit cycle. In the case illustrated by Fig. 1, the value of $`m_\pi `$ at which $`a_t`$ diverges is larger than that at which $`a_s`$ diverges. If they both diverged at the same value of $`m_\pi `$, there would be an exact infrared limit cycle.

We conjecture that QCD can be tuned to this infrared limit cycle by adjusting the up and down quark masses $`m_u`$ and $`m_d`$. As illustrated in Fig. 1, the tuning of $`m_\pi `$, which corresponds to $`m_u+m_d`$, is likely to bring the system close enough to the infrared limit cycle for the triton to have one excited state. We conjecture that by adjusting the two parameters $`m_u`$ and $`m_d`$ to critical values, one can make $`a_t`$ and $`a_s`$ diverge simultaneously. At this critical point, the deuteron and spin-singlet deuteron would both have zero binding energy and the triton would have infinitely-many increasingly-shallow excited states. The ratio of the binding energies of successively shallower states would rapidly approach a constant $`\lambda _0^2`$ close to 515. It may be possible to use a combination of lattice gauge theory and EFT to demonstrate the existence of this infrared RG limit cycle in QCD.