text:
# Metallic mean-field stripes, incommensurability and chemical potential in cuprates

## Abstract

We perform a systematic slave-boson mean-field analysis of the three-band model for cuprates with first-principle parameters. Contrary to widespread believe based on earlier mean-field computations low doping stripes have a linear density close to 1/2 added hole per lattice constant. We find a dimensional crossover from 1D to 2D at doping $`0.1`$ followed by a breaking of particle-hole symmetry around doping 1/8 as doping increases. Our results explain in a simple way the behavior of the chemical potential, the magnetic incommensurability, and transport experiments as a function of doping. Bond centered and site-centered stripes become degenerate for small overdoping.

It is now a well established fact that doped holes in cuprates tend to self-organize in antiferromagnetic (AF) domain wallsTranquada et al. (1995, 1997); Yamada et al. (1998); Arai et al. (1999, 2000). These quasi one-dimensional (1D) structures called stripes where predicted by mean-field theoriesZaanen and Gunnarsson (1989); Poilblanc and Rice (1989) inspired by the problem of solitons in conducting polymersSu et al. (1979).

In the superconducting systems observed stripes are parallel to the Cu-O bond (hereafter called vertical stripes) and for low doping have a linear filling fraction $`\nu 1/2`$ of added holes per lattice constant ($`a1`$). Unless a one-dimensional (1D) instability opens a gap at the Fermi levelZaanen and Oleś (1996), half-filled ($`\nu =1/2`$) stripes are expected to be metallic. On the contrary mean-field theories predicted $`\nu =1`$ insulating stripes which led to an early rejection of stripesCheong et al. (1991).

Here we perform a systematic analysis of slave-boson (SB) mean-fieldKotliar and Ruckenstein (1986) solutions for a three-band Hubbard model as appropriate for the cuprates. Contrary to the early mean-field analysis we find that for low doping $`x`$ in a material like La2-xSrxCuO4 (LSCO) the most favorable stripes have indeed $`\nu 1/2`$ (and are vertical) reconciling the mean-field picture with experiment. Small doping introduces more stripes at fixed $`\nu `$ reducing the stripe distance $`d`$ and explaining the observed behavior of the magnetic incommensurability $`ϵx`$ with $`ϵ=1/(2d)`$Tranquada et al. (1995, 1997); Yamada et al. (1998); Arai et al. (1999, 2000).

We find low doping stripes to be centered on oxygen, bridging two Cu sites, hereafter named bond centered (BC) stripes after the terminology of one-band modelsWhite and Scalapino (1989a); Fleck et al. (2000) (see Fig. 1).

Significant interstripe overlap sets in at $`x0.1`$ triggering a crossover from 1D behavior at low doping to 2D behavior at large doping. This does not affect the stability of stripes up to $`x=1/8`$ where a change of regime occurs. Further doping occurs at a fixed stripe distance $`d=4`$ and $`\nu `$ starts to grow. This explains the change of behavior observed in both the incommensurabilityTranquada et al. (1997); Yamada et al. (1998); Arai et al. (1999, 2000) and the chemical potentialIno et al. (1997); Harima et al. (2001) around doping $`x=1/8`$. Furthermore for $`x<1/8`$ the chemical potential is at the center of an approximately particle-hole symmetric band (Fig. 2) whereas for doping $`x>1/8`$ particle-hole symmetry is broken in good agreement with the picture deduced from transport experimentsNoda et al. (1999); Wang and Ong (2001). Finally at $`x0.21`$ BC stripes and site-centered (SC) stripes become degenerate suggesting a regime of soft lateral stripe fluctuations possibly relevant for superconductivity. This quasidegeneracy was found in more sophisticated treatments of (less realistic) one-band modelsWhite and Scalapino (1989a); Fleck et al. (2000) however the ordering of the states was reversedFleck et al. (2000).

One-band models can provide qualitative understanding but may miss subtle but important details specific for the cuprates. Therefore we use a three-band Hubbard model with a parameter set deduced from a first-principle LDA computation for LSCOMcMahan et al. (1990). We view our computation as the last stage of a first principle computation and thus without free parameters. The largest repulsion $`U_d`$McMahan et al. (1990) is treated within a SB approachKotliar and Ruckenstein (1986) fully equivalent to the Gutzwiller approximation and the other interactions are treated within the Hartree-Fock (HF) approximation (an equivalent formalism has been used in Ref. Sadori and Grilli (2000) for electron-doped cuprates). We solve the mean-field equations both in real and in momentum space. The real space computations are completely unrestrictedSeibold et al. (1998) except for the spins which are assumed to be along the $`z`$ axis.

Typically we minimize the energy in large clusters with $`N_s`$ stripes ($`N_s=2,4`$) of a maximum length $`L=30`$. Depending on boundary conditions, size and filling, solutions with and without charge order in the stripe direction are obtained. In general due to the 1D character of stripes one expects that a 1D instability will develop in large systems at low doping rendering the system insulating even for $`\nu <1`$. Indeed if superconductivity is suppressed by pulsed magnetic fields the system is insulating below a critical doping and a crossover temperature $`T^{}`$Boebinger (1998). In this work for simplicity we concentrate on the metallic state and restrict to solutions without charge order in the stripe direction. Since $`T^{}`$ is quite small compared with the typical electronic scales we can set the temperature to zero in our computations.

Fig. 1 depicts the typical charge and spin distribution for a BC vertical stripe. Added charge is accumulated in the antiphase boundary between the AF regions and mainly resides on oxygen sites.

In Fig. 2 we show the band structure for the electrons in an extended zone scheme. Two bands appear well inside the charge transfer gap. The upper one at $`0.9`$ eV is more Cu like and quite flat (even in the stripe direction) whereas the lower one has substantial dispersion, is mainly O like, and crosses the Fermi level. Counting both ingap bands the $`\nu =1/2`$ system corresponds to $`3/4`$ hole ($`1/4`$ electron) band filling. The gap between the bands is related to the magnetic character of the stripe. (It is absent in SC solutions with a non-magnetic core.) The electronic structure close to the chemical potential is well represented by a half-filled cosine like band, hereafter referred to as the active band.

The dispersion of the active band perpendicular to the stripe is quite flat indicating quasi 1D behavior. Small oscillations appear in this direction for $`d5`$ which can be identified as a crossover interstripe distance from 1D to 2D behavior and corresponds to a crossover doping $`x=0.1`$ for $`\nu =0.5`$.

The band structure in the stripe direction (right panel) is quite similar in all BC stripes solutions regardless of $`d`$ which instead determines the periodicity in the perpendicular direction (left panel). Comparison with photoemission data requires matrix elements considerations and will be presented elsewhere. Here we notice that the Fermi surface crossing at vertical momentum $`k_y\pi /4`$ has been recently observed for $`d=4`$ stripesZhou et al. (1999) and the flat portion of the active band in the direction perpendicular to the stripe correlates well with the flat bands observed in many cuprates around $`𝐤(\pi ,0)`$M. King et al. (1994).

To evaluate the stability we compute the energy per added hole for $`N_h`$ holes added to the system with $`N`$ CuO2 units, $`e_h[E(N+N_h)E(N)]/N_h`$. Here $`E(N+N_h)`$ \[$`E(N)`$\] is the total energy of the doped \[undoped AF\] solution. In Fig. 3 (a) we show $`e_h`$ as a function of the filling fraction $`\nu N_h/(N_sL)`$ for BC vertical stripes. Each curve corresponds to a fixed $`d`$. The curves have a sharp minimum at $`\nu =0.5`$. This is in contrast to early one-bandPoilblanc and Rice (1989); Zaanen and Oleś (1996) and three-band mean-field computationsZaanen and Gunnarsson (1989) for which $`\nu =1`$ is the most favorable hole filling. We can reproduce this behavior by setting artificially $`t_{pp}=0`$ in our parameter set (crosses). The difference in behavior can be understood by analyzing the band structure. In the case $`t_{pp}=0`$ the active band becomes as narrow as the 0.9eV band (Fig. 2) and a gap develops between the active band and the band immediately bellow instead of overlapping as in Fig. 2. This gap produces the cusp singularity at $`\nu =1`$ \[Fig. 3 (a)\] via the discontinuity in $`E/N`$ which makes $`\nu =1`$ very stable. As $`t_{pp}`$ is increased the gap closes, the cusp singularity disappears, and a minimum develops in the $`\nu <1`$ part of the curve. A simple model explains the formation of this minimum. Consider a flat density of states (DOS) of width $`W`$ to estimate the kinetic energy of the active band and neglect interactions. In this approximation $`e_h=e_0/\nu +W\nu /2`$ where $`e_0`$ is the energy cost per unit length to create a domain wall in the undoped system ($`e_0J/2`$ with $`J`$ the superexchange constant). The optimum filling is given by $`v_{\mathrm{opt}}=\mathrm{min}(\sqrt{2e_0/W},1)`$. Increasing $`t_{pp}`$ increases $`W`$ without substantially affecting $`e_0`$ and produces a local minimum for $`\nu <1`$. The stability close to $`\nu =1/2`$ is further enhanced because the DOS of a 1D band is not flat but has a minimum at $`\nu =1/2`$. For the same reason a gap produces a cusp in the energy, a minimum of the DOS enhances the local curvature of $`e_h`$ vs. $`\nu `$ and tends to shift the optimum filling closer to $`\nu =1/2`$.

We have checked that for the most common parameter sets which contain $`t_{pp}`$McMahan et al. (1990) results are similar. However, as soon as $`t_{pp}`$ is set to zero, as in Ref. Zaanen and Gunnarsson (1989), half-filled stripes become unstable in favor of filled ones.

The $`t_{pp}=0`$ result is in accord with HF computations in the simplest one-band model. Including a second-neighbor hopping ($`t^{}`$) we have found that partially filled stripes can be stabilized also in one-band models by a similar mechanism, which will be reported elsewhere.

Partially filled stripes where also found in the one-band Hubbard model with $`t^{}=0`$ in a much stronger correlation regime and going beyond a static mean-fieldFleck et al. (2000). In our case, because the active band is broad, correlations are less important and we expect that a static mean-field provides a good starting point.

Vertical stripes should be compared with other possible ground states. We find that within the present SB approach they are lower in energy than polaron solutionsLorenzana and Yu (1993) \[see Fig. 3 (b)\] and diagonal stripes solutions (not shown). Very low nonsuperconducting dopings ($`x<0.05`$) were not explored in detail since we believe that a careful consideration of other effects is required in this case. Especially long-range order perpendicular to the planes favors loop configurations of stripesWhite and Scalapino (1989a) and long-range Coulomb effectsLöw et al. (1994); Lorenzana et al. (2001a) are expected to become important.

If one uses the HF approximation instead of SB the minimum also occurs at $`\nu 1/2`$ for BC vertical stripes. However a polaron lattice is the ground state in HF and diagonal stripes are lower in energy than vertical ones. Thus the overall success of the present mean-field computation with respect to the earlier onesZaanen and Gunnarsson (1989) is due to both a more accurate mean-field approximation SB instead of HF) and a more accurate parameter set.

Further stabilization of the half-filled stripes can occur if due to many-body effects a gap or pseudogap tends to open at the commensurate filling $`\nu =1/2`$Zaanen and Oleś (1996). Although we do not need this effect to explain $`\nu 1/2`$ it is quite possible that this produces a fine tuning for $`T<T^{}`$. In this regard it is interesting to remark that the “V” shape form of the curves in Fig. 3 (a) is due to the gap produced by the discretization of the levels in a finite system. For larger systems \[see $`d=4`$ curve in Fig. 3 (a)\] the cusp becomes rounded, however the minimum is still close to $`\nu =1/2`$. One should be aware that the curves with the cusp are for already quite long stripes ($`L=12`$) and it is not clear whether stripes in real materials will be much longer. Thus finite size data could turn out to be more realistic than infinite size one.

For $`d>4`$ all curves coincide close to $`\nu =0.5`$ in Fig. 3 (a) whereas for $`d=3,4`$ the curves are shifted up. This feature is due to the width of the domain wall of $`45`$ lattice sites (see Fig. 1) which forces stripes to overlap leading to an increase of $`e_h`$.

In Fig. 3 (b) we show $`e_h`$ for various values of the stripe separation $`d`$ and as a function of doping $`x`$ assuming each dopant introduces one hole. Since $`x=\nu /d`$ curves with larger $`d`$ appear at lower concentration. The locus of the minimum of $`e_h`$ as a function of doping in Fig. 3(b) is expected to form a continuous curve in the thermodynamic limit by combining different $`d`$ solutions in new solutions with larger periodicity. Up to $`x1/8`$ the stripe filling is fixed at $`\nu 0.5`$ and consequently the density of stripes increases with doping. This explains the behavior of the incommensurability $`ϵ=x/(2\nu )x`$ as seen in neutron scattering experiments in this doping range (inset of Fig. 3)Tranquada et al. (1995, 1997); Yamada et al. (1998); Arai et al. (1999, 2000). For $`x>1/8`$ the right branch of the $`d=4`$ solution is more stable than the $`\nu 0.5`$ and $`d=3`$ solution due to the stripe overlap effect. Therefore the incommensurability remains locked at $`ϵ1/8`$ in good agreement with the change of behavior in $`ϵ`$ observed around $`x1/8`$Tranquada et al. (1997); Yamada et al. (1998); Arai et al. (1999, 2000).

As doping increases further, BC stripes become degenerate with SC ones at $`x0.21`$. This suggests that as doping increases lateral fluctuations of the stripe will become soft and possibly mediate pairing between holes.

For doping $`x>0.225`$ we find the $`d=3`$ stripe ($`ϵ=1/60.17`$) to become the lowest energy solution with an initial filling fraction $`\nu =0.675`$. Experimentally the situation is not clear. $`ϵ=1/6`$ has been reported for YBa2Cu3O6+δArai et al. (2000) (YBCO) but not for LSCO where $`ϵ`$ remains in the $`ϵ=1/8`$ line up to $`x=0.25`$ (inset of Fig. 3). In La2-x-yNdySrxCuO4 (LNSCO), where stripes are pinned by the low-temperature tetragonal lattice distortion, the incommensurability is substantially increased beyond $`ϵ=1/8`$ at doping $`x=0.2`$ but without reaching $`ϵ=1/6`$Tranquada et al. (1997).

It is possible that lateral stripe fluctuations become so strong at this doping range that an effectively isotropic state is reached Kivelson et al. (1998). Another possibility is that the $`d=3`$ phase is skipped due to phase separation among the $`d=4`$ phase and the overdoped Fermi liquid which for our parameters becomes the lowest energy solution close to $`x=0.4`$ (see Ref. Uemura (2001) for a related scenario). Further theoretical and experimental work should be done to clarify this point. Especially being the system charged a careful analysis of phase separation is neededLorenzana et al. (2001a).

It has been recently emphasizedHarima et al. (2001) that there is a close connection between the doping dependent incommensurability as discussed above and the chemical potential in cuprates. In fact the chemical potential for the electrons can be related to $`e_h`$ via: $`\mu =(e_h+x\frac{e_h}{x})`$. From Fig. 3(b) one can then deduce that $`\mu `$ is approximately constant for $`x0.1`$ and decreases for $`x0.1`$ in qualitative agreement with the observed behaviorIno et al. (1997); Harima et al. (2001). The rate of change of $`\mu `$ with doping, being a high derivative of the energy, is very sensitive to finite size effects and, moreover, few experimental points are available in this doping range in order to allow for a precise comparison. A rough estimate indicates that the theoretical rate of change of $`\mu `$ with doping for $`x>1/8`$ is approximately a factor of 2 larger than the experimental oneIno et al. (1997); Harima et al. (2001). This may be due to phase separation among the $`d=4`$ stripe solution and the paramagnetic overdoped Fermi liquid as mention above.

We finally turn to the discussion of our calculation in light of recent measurements of the Hall coefficient RH in LNSCONoda et al. (1999) and YBCOWang and Ong (2001). In the former compound RH displays an abrupt decrease below the charge-stripe ordering temperature T0 and for concentrations $`x1/8`$. It has convincingly been argued that quasi 1D transport is not enough to explain this anomaly and instead reflects a remarkable cancellation due to particle-hole symmetry in the stripe stateWang and Ong (2001); Prelovšek et al. (2001). A partial suppression of RH (and simultaneously the thermopower S) below some temperature Tmax has also been observed in YBCO up to oxygen contents corresponding to doping $`x1/8`$ Wang and Ong (2001). For distant stripes we indeed observe that the chemical potential crosses an approximately particle-hole symmetric band (see Fig. 2). Interstripe hopping becomes significant for $`x>0.1`$ however it does not break particle-hole symmetry. As discussed above the chemical potential shifts from the center of the band for $`x>1/8`$, thus breaking particle-hole symmetry and consequently both $`R_H`$ and $`S`$ start to grow in modulus. Unfortunately the magnitude and even the sign of the transport coefficients cannot be predicted by a knowledge of the band structure alone. For example $`R_H`$ for our Fermi surface (which turns out to be open) is determined by the unknown anisotropy of the scattering path lengthOng (1991).

To conclude we have shown that, contrary to widespread believe, an appropriate mean-field theory of stripe phases in cuprates is not in contradiction with experiments. Our theory explains in a simple way the change of behavior of several experimental quantities (incommensurability, chemical potential, transport coefficients) around the “magic” doping $`x=1/8`$. Obviously the present stripes are static whereas those in real materials are dynamic unless pinned by a structural distortion. Many properties however will be largely independent of the dynamical or static character of the stripes. The energy for example is determined by short range correlations and is insensitive to the long-distant behavior. This is analogous to the problem of spin waves in a 1D antiferromagnet where the spin wave approximation gives a reasonable value for the energy although the system does not have true long-range orderAnderson (1952). Since the properties we have discussed are mainly determined by the energetics we believe that the present mean-field theory is a good starting point for an understanding of the electronic structure of hole doped cuprates.


format: text

text:
# Non-equilibrium chemistry in the atmospheres of brown dwarfs

## 1. Introduction

The discovery of strong methane bands in the spectrum of Gl 229B indicated at once that it has a very low effective temperature ($`1000`$K) and it was hailed as the first clearly identified brown dwarf (Nakajima et al. 1995; Oppenheimer et al. 1995). This follows from the chemistry of carbon in stellar atmospheres. In cool M dwarfs, carbon is predominantly in the form of carbon monoxide (CO), while in Jovian planets methane (CH4) is the dominant form. The transition occurs in the brown dwarf regime, at $`T1100`$K for a pressure of 1 bar (Fegley & Lodders 1996).

The presence of methane bands is so striking that it prompted the creation of the new spectral class of T dwarfs (Kirkpatrick et al. 1999). Gl 229B is a T6 dwarf. In this context, the discovery of a strong CO band in the 4.5 – 5$`\mu `$m spectrum of Gl 229B by Noll, Geballe & Marley (1997) is quite remarkable. Analysis of the spectrum reveals an abundance of CO that is over 3 orders of magnitude larger than expected from chemical equilibrium calculations (Noll et al. 1997; Griffith & Yelle 1999; Saumon et al. 2000), as anticipated by Fegley & Lodders (1996). A similar situation is observed in the atmosphere of Jupiter where CO is detected while chemical equilibrium overwhelmingly favors CH4 as the reservoir of carbon (Prinn & Barshay 1977).

The excess of CO in Gl 229B can be understood if the chemistry of carbon is driven away from equilibrium by vertical mixing on a time scale shorter than the rate of conversion of CO into CH4 (Fegley & Lodders 1996; Noll et al 1997; Griffith & Yelle 1999).

The two reservoirs of nitrogen in brown dwarfs are molecular nitrogen (N2) and ammonia (NH3). Ammonia is the dominant molecule at low temperatures. While N2 is invisible in the near infrared, Saumon et al. (2000) detected NH3 features in Gl 229B in the $`K`$ band and inferred its presence in the $`H`$ band. Each detection corresponds to a different level in the atmosphere. NH3 abundances determined separately for each band reveal a value in agreement with chemical equilibrium in the $`H`$ band and a depletion by at least a factor of 4 in the $`K`$ band. These results can be interpreted consistently with the CO observation by invoking vertical mixing in the atmosphere on a time scale shorter than the rate of conversion of N2 into NH3.

The modeling of non-equilibrium chemistry of CO and other molecules due to mixing in planetary atmospheres is a mature field (e.g. Barshay & Lewis 1978; Fegley & Prinn 1985; Fegley & Lodders 1994; Lodders & Fegley 1994). Application to brown dwarfs have so far been limited to modeling the excess of CO in Gl 229B (Griffith & Yelle 1999) and a detailed study of the kinetics of CNO chemistry (Lodders & Fegley 2002). Here we present the first systematic exploration of the effects of mixing on the chemistry of carbon and nitrogen in brown dwarfs atmospheres. We first summarize the main features of non-equilibrium chemistry of carbon and nitrogen caused by vertical mixing. Using a simple parametrization for the rate of mixing, we then compute the resulting non-equilibrium abundances of CO, CH4, N2, NH3 and H2O as a function of effective temperature and mixing efficiency. The new abundance profiles lead to significant changes in the emergent spectra of cool brown dwarfs.

## 2. Carbon and Nitrogen Chemistry

The relative abundances of CO and CH4 can be driven out of equilibrium in the cool atmospheres of brown dwarfs and giant planets because of the strong asymmetry in the reaction rates of the (net) reaction

$$\mathrm{CO}+3\mathrm{H}_2\mathrm{CH}_4+\mathrm{H}_2\mathrm{O}.$$
$`(1)`$
Because of the large binding energy of CO, the reaction proceeds much more slowly to the right than to the left. Vertical mixing dredges up hot gas that is relatively rich in CO to the cooler part of the atmosphere, where CO is converted very slowly into CH4, while the CH4 that is brought down is quickly converted into CO. The net result is an overabundance of CO, with the CO/CH4 ratio in the upper atmosphere fixed to its value at the level where the mixing time scale and the chemical time scale are equal. Nearly all of the elemental oxygen is partitioned between condensates, carbon monoxide and water. It follows that an overabundance of CO implies an underabundance of H2O.

The conversion of N2 (favored at higher temperatures) into NH3 is entirely analogous to that of CO/CH4. The net reaction is

$$\mathrm{N}_2+3\mathrm{H}_22\mathrm{N}\mathrm{H}_3.$$
$`(2)`$
The large binding energy of N2 leads to a slow reaction rate to the right, while the reaction to the left proceeds rapidly. Thermochemical and kinetic data indicate that all other molecular and atomic species identified in the spectra of brown dwarfs react on very short time scales and are expected to be in chemical equilibrium. In brown dwarfs, non-equilibrium chemistry induced by vertical mixing is thus limited to the five species CO, CH4, H2O, NH3 and N2.

The correct chemical pathway and reaction time scale for the conversion of CO into CH4 remains somewhat uncertain (Griffth & Yelle 1999; Lodders & Fegley 2002). For the purposes of this exploratory calculation, however, this is not a critical consideration and we adopt the chemical time scales for the net reactions (1) and (2) given in Lodders & Fegley (2002).

## 3. Vertical Transport

In radiative atmospheres, vertical mixing can occur by eddy turbulence. This process is analogous to diffusion and occurs over a characteristic time scale

$$\tau _{\mathrm{mix}}\frac{H^2}{K},$$
$`(3)`$
where $`H`$ is the pressure scale height and $`K`$ is the coefficient of diffusion. While $`K`$ is a free parameter in our calculation, it ranges from $`10^2`$ to $`10^5`$cm2/s in planetary stratospheres.

Atmospheres of brown dwarfs all become convective below a certain depth and convection is a very effective mixing mechanism. The time scale for convective mixing is

$$\tau _{\mathrm{mix}}=\tau _{\mathrm{conv}}\frac{H_c}{v_c}=\frac{H^2}{K_c},$$
$`(4)`$
where $`H_c`$ and $`v_c`$ are the convective mixing length and velocity, respectively. For convenience, we choose $`H_c=H`$. An effective diffusion coefficient, $`K_c`$, can be defined for convective mixing (Eq. 4) and evaluated using the mixing length formalism. In brown dwarf atmospheres, we find $`K_c10^810^9`$cm2/s. Thus, the transition from a convective zone to a radiative zone in the atmosphere is accompanied by a sudden decrease in the efficiency of mixing.

## 4. Non-Equilibrium Chemistry in Brown Dwarf Atmospheres

We model the effects of vertical transport on the chemistry of carbon and nitrogen in cloudless atmospheres with $`700T_{\mathrm{eff}}2000`$K, $`g=10^5`$cm/s2, and solar abundances. The eddy diffusion coefficient in the radiative zone is varied from $`K=10^2`$ to $`10^{10}`$cm2/s while in the convective zone the mixing time scale is calculated from the mixing length formalism.

An example of the effect of vertical mixing on the chemical abundance profiles in a model atmosphere is shown in Figure 1. Depth in the atmosphere is indicated by the local temperature. The equilibrium abundances, shown with dashed lines, clearly show the transition in carbon chemistry from CO deep in the atmosphere to CH4 near the surface. The oxygen freed by the conversion of CO into CH4 toward the top of the atmosphere is used to form H2O. The H2O abundance thus increases in concert with that of CH4. In this relatively hot model, nearly all nitrogen is in the form of N2 and the NH3 abundance remains very low throughout the atmosphere.

The chemical time scales for the conversion of CO and N2 are shown by the two heavy solid lines. These vary by more than 25 orders of magnitude throughout the atmosphere! Mixing occurs rapidly in the convection zone but is considerably less efficient in the radiative zone for our choice of $`K=100`$cm2/s. In regions where $`\tau _{\mathrm{chem}}<\tau _{\mathrm{mix}}`$, chemical equilibrium prevails. Where $`\tau _{\mathrm{chem}}>\tau _{\mathrm{mix}}`$, however, the mole fractions are determined by their values at the crossing point. As can be seen in Fig. 1, there can be more than one level where $`\tau _{\mathrm{chem}}=\tau _{\mathrm{mix}}`$, resulting in alternating zones in and out of chemical equilibrium. For simplicity, we consider only the uppermost (lowest T) crossing. It turns out that in all cases, this is an excellent approximation to the more complex, full solution because the species affected have either 1) very low abundances or 2) have nearly flat abundance profiles where the deeper crossings occur and are unaffected by mixing.

The resulting non-equilibrium abundances are shown by solid lines in Fig. 1. The abundance of CO is considerably enhanced by vertical transport while both CH4 and H2O are reduced by more modest factors. N2 is “frozen” at a deeper level than CO. Its abundance is barely affected in this model, while the NH3 abundance decreases. The net effect of vertical transport in all models is an increase in the abundances of CO and N2, and a decrease in the abundances of CH4, H2O and NH3. Shortening the mixing time scale in the radiative zone (i.e. increasing $`K`$) produces larger changes in the abundances for fixed $`T_{\mathrm{eff}}`$ and gravity.

The effect of non-equilibrium chemistry on the spectra of brown dwarfs is shown in Figure 2. The change in the strength of the 4.7$`\mu `$m band of CO increases gradually toward lower $`T_{\mathrm{eff}}`$ up to a maximum at $`T_{\mathrm{eff}}900`$K. The corresponding depletion in CH4 results in weaker CH4 bands with a maximum effect at $`T_{\mathrm{eff}}1400`$K. The 3.3$`\mu `$m band of CH4 is the most strongly affected, followed by the 2.2$`\mu `$m and 1.6$`\mu `$m bands.

The strongest features of NH3 appear in the 10–11$`\mu `$m region and have not yet been observed in any brown dwarf. Ammonia forms only at relatively low temperatures ($`T_{\mathrm{eff}}<1200`$K) and the strength of the features increases monotonically with decreasing $`T_{\mathrm{eff}}`$. Vertical mixing systematically reduces the abundance of NH3 compared to its equilibrium value. Since NH3 is detectable at three different levels in the atmosphere (in the $`N`$, $`K`$, and $`H`$ bands; Saumon et al. 2000), the flat abundance profile predicted by this model (Fig. 1) can in principle be tested observationally.

Saumon et al. (2000) argued that the most reliable diagnostic of metallicity for T dwarfs at present is the forest of H2O lines in regions of the spectrum where water is the dominant absorber. A determination of the metallicity based on H2O lines would be affected by non-equilibrium chemistry, which can reduce its abundance by up to $`0.3`$ dex (increasing with $`K`$ and $`T_{\mathrm{eff}}`$).

## 5. Observational Evidence

Gl 229B is the only T dwarf with a published 4.5 – 5$`\mu `$m spectrum (Noll et al. 1997; Oppenheimer et al. 1998) but $`M^{}`$ band photometry is now available for several other brown dwarfs (Leggett et al. 2002; Golimowski et al. 2002). The $`M^{}`$ band is a very good probe of the strength of the CO band (Fig. 2). Non-equilibrium spectra are compared with the spectrum of Gl 229B in Figure 3. A fair agreement with the rather noisy 4.5 – 5$`\mu `$m spectrum is obtained for $`K10^2`$cm2/s.

Absolute $`M^{}`$ photometry is shown in Figure 4. The equilibrium models (dashed curve) overestimate the $`M^{}`$ flux at lower $`T_{\mathrm{eff}}`$. Non-equilibrium models can reproduce the observed fluxes quite well with $`K10^210^4`$cm2/s. Figure 4 suggests that increased CO absorption resulting from vertical transport is not a peculiarity of Gl 229B but a feature common to late L dwarfs and T dwarfs.

## 6. Conclusions

Motivated by observations in the T6 dwarf Gl 229B of a strong CO band where none was expected and of a depletion of NH3, we have explored the effects of vertical mixing in brown dwarf atmospheres on the most prominent molecules visible in their spectra. Vertical mixing, by convection or eddy turbulence, can occur fast enough to keep abundances of CO, CH4, H2O and NH3 from reaching their values at chemical equilibrium. Increased efficiency of mixing leads to an enrichment of CO and a depletion of CH4, H2O, and NH3. The model of vertical mixing can explain the CO and NH3 observations in Gl 229B self-consistently. Furthermore, observations of NH3 in three different bandpasses can test the model prediction of a flat abundance profile throughout the atmosphere (Saumon et al. 2002). In this study, the efficiency of mixing is treated as a free parameter but it can eventually be calibrated with data.

New observational evidence suggests that the overabundance of CO observed in Gl 229B may be a common feature of late L dwarfs and of T dwarfs. If this is indeed the case, it becomes essential to include non-equilibrium chemistry in modeling the atmospheres and spectra of cool brown dwarfs. Because it affects the relative abundances of major constituents of the atmosphere, vertical mixing further complicates the determination of the composition of brown dwarfs.

Finally, vertical mixing can result in a suppression of the $`M^{}`$ flux by up to one magnitude for $`T_{\mathrm{eff}}`$ between $`700`$ and 1100 K (Fig. 4). This must be taken into account in plans to image extrasolar giant planets by capitalizing on their predicted surperthermal $`M`$ band flux (Marley et al. 1996, Burrows et al. 1997).

### Acknowledgments.

We thank S. K. Leggett and D. A. Golimowski for kindly providing data in advance of publication. This work was supported in part by NSF grant AST-0086288 and NASA grant NAG5-9273 to M. S. M. and grant NAG5-4970 to R. S. F. Work of K. L. is supported by NASA grant NAG5-11958 and NSF grant AST-0086487.

## References

Barshay, S. S. & Lewis, J. S. 1978, Icarus, 33, 593

Burrows, A., Marley, M. S., Hubbard, W. B., Lunine, J. I., Guillot, T., Saumon, D., Freedman, R. S., Sudarsky, D., & Sharp, C. 1997, ApJ, 491, 856

Fegley, B. Jr., & Lodders, K. 1994, Icarus, 110, 117

Fegley, B. Jr., & Lodders, K. 1996, ApJ, 472, L37

Fegley, B. Jr., & Prinn, R. G. 1985, ApJ, 299, 1067

Golimowski, D. A., 2002, in preparation

Griffith, C. A. & Yelle, R. V. 1999, ApJ, 519, L85

Kirkpatrick, J. D., Reid, I. N., Liebert, J., Cutri, R. M., Nelson, B., Beichman, C. A., Dahn, C. C., Monet, D. G., Gizis, J. E. & Skrutskie, M. F. 1999, ApJ, 519, 834

Leggett, S. K. et al. 2002, ApJ, 564, 452

Leggett, S.K., Toomey, D.W., Geballe, T.R. & Brown, R.H. 1999, ApJ, 517, L139

Lodders, K. & Fegley, B. Jr. 1994, Icarus, 112, 368

Lodders, K. & Fegley, B. Jr. 2002, Icarus, 155, 393

Marley, M. S., Saumon, D., Guillot, T., Freedman, R. S., Hubbard, W. B., Burrows, A. & Lunine, J. I. 1996, Science, 272, 1919

Nakajima, T., Oppenheimer, B. R., Kulkarni, S. R., Golimowski, D. A., Matthews, K., & Durrance, S. T. 1995, Nature, 378, 463

Noll, K. S., Geballe, T. R. & Marley, M. S. 1997, ApJ, 489, L87

Oppenheimer, B. R., Kulkarni, S. R., Matthews, K. & Nakajima, T. 1995, Science, 270, 1478

Oppenheimer, B. R., Kulkarni, S. R., Matthews, K. & van Kerkwijk, M. H. 1998, ApJ, 502, 932

Prinn, R. G. & Barshay, S. S. 1977, Science, 198, 1031

Saumon, D., Geballe, T. R., Leggett, S. K., Marley, M. S., Freedman, R. S., Lodders, K., Fegley, B. Jr. & Sengupta, S. K. 2000, ApJ, 541, 374


text:
# Sole-Manrubia model of biological evolution: some new insights

## I Introduction

Several fundamental questions on the origination, evolution and extinction of species have been addressed in the recent years by theoretical treatments of dynamic network models using the concepts and techniques of non-equilibrium statistical mechanics drossel ; newman . Kauffman kauffman  pioneered a model which describes biological evolution as a “walk” in a “fitness landscape” and argued that the eco-system adaptively evolves so as to be delicately poised “at the edge of chaos”. On the other hand, Bak and Sneppen bak , suggested an alternative way of implementing the Darwinian principle of selection (the “survival of the fittest”); this dynamics naturally leads the eco-system to a “self-organized” critical statejensen . These models implicitly assume a primarily biotic cause of extinctions. However, models based on purely abiotic causes, which arise from “environmental stresses” also lead to similar pattern of extinctions of species newman ; newman2 .

A few years ago Sole and Manrubia (SM) sole  introduced one such model of “macro”-evolution that takes into account the inter-species interactions explicitly. They reported several statistical properties of the model as well as their biological implications. In this paper we report some other statistical properties of this model; some of these features are qualitatively similar to the corresponding results obtained in a more recent model chowstau . Moreover, we re-interpret some of the quantities in their approach to extract some relevant biological implications of the model. Furthermore, we introduce a new correlation function throwing light on the “volatility” of the fitness and the level of “coherence” in the changes of fitness of the species in this model.

## II The model

For the sake of completeness, we briefly present the SM model and its main published results.

The system consists of $`N`$ species, each labelled by an index $`i`$ ($`i=1,2,\mathrm{}N`$). The states of the $`i`$-th species is represented by a two-state variable $`S_i`$; $`S_i=0`$ or $`1`$ depending on whether it is extinct or alive, respectively. The inter-species interactions are captured by the interaction matrix $`𝐉`$; the element $`J_{ij}`$ denotes the influence of the species $`j`$ on the species $`i`$. If $`J_{ij}>0`$ while, simultaneously, $`J_{ji}<0`$ then $`i`$ is the predator and $`j`$ is the prey. On the other hand, if both $`J_{ij}`$ and $`J_{ji}`$ are positive (negative) they the two species cooperate (compete).

The dynamics of the system consists in updating the states of the system (i.e., to determine the state at the time step $`t+1`$ from a complete knowledge of the state at the time $`t`$) in the following three steps:

Step (i): One of the input connections $`J_{ij}`$ for each species $`i`$ is picked up randomly and assigned a new value drawn from the uniform distribution in the interval $`[1,1]`$, irrespective of its previous magnitude and sign.

Step (ii): The new state of each of the species is decided by the equation

$$S_i(t+1)=\mathrm{\Theta }\left(\underset{j=1}{\overset{N}{}}J_{ij}S_j(t)\theta _i\right)$$
(1)
where $`\theta _i`$ is a threshold parameter for the species $`i`$ and $`\mathrm{\Theta }(x)`$ is the standard step function, i.e., $`\mathrm{\Theta }(x)=1`$ if $`x>0`$ but zero otherwise. If $`S(t+1)`$ becomes zero for $`m`$ species, then an extinction of size $`m`$ is said to have taken place.

Step (iii): All the niches left vacant by the extinct species are refilled by copies of one of the randomly selected non-extinct species.

Sole and Manrubia recorded extinctions of sizes as large as $`500`$ and the distributions of the sizes of these extinctions could be fitted to a power law of the form $`N(m)m^\alpha `$ with an exponent $`\alpha 2.3`$. Moreover the periods of stasis $`t_s`$ were also found to obey a power law $`N(t_s)t_s^\gamma `$ with the exponent $`\gamma 3.0`$. However, surprisingly, in none of their papers sole , Sole and his collaborators reported the distributions of the lifetimes of species which, according to some claims (see, for example, refs.drossel ; newman  for references to the experimental literature and data analysis), also follows a power law.

## III Results and discussions

In most of the food webs recorded so far only the sign (but not the magnitude) of the inter-species interactions are noted. Moreover, all the qualitative features of the SM model remain unchanged if, instead of allowing the elements $`J_{ij}`$ to take any real value in the domain $`[1,1]`$ with equal probability, we allow $`J_{ij}`$ to be only $`+1`$ or $`1`$ with equal probability. Therefore, we work with the simpler version of the SM model where each of the elements of $`𝐉`$ are discrete binary variables that can be either $`+1`$ or $`1`$.

In our simulations we let the system “warm up” for $`t_w`$ time steps, starting from a random initial condition, and then run it for further $`t_m`$ time steps during which collect the data.

### III.1 Distribution of lifetimes

It has been claimed in the literature drossel ; newman  that the distributions of the lifetimes of the extinct species, estimated from the fossil data, can be fitted to a power law of the form

$$P_t(\tau )\tau ^x$$
(2)
where several possible values of $`x`$, ranging from $`1`$ to $`2`$ have been suggested. However, the log-log plot of our data in fig.1, obtained from the simulation of the SM model, demonstrates that an approximate power-law-like behaviour is seen only over nearly one decade of short lifetimes while strong deviations from power-law are seen for longer lifetimes. While the sharp drop of the curves in the tail region can be attributed to finite time effects, the deviations from the early power-law-like region at around $`\tau 10`$ is neither an artefact of finite time of simulation nor that of finite size of the system chow .

### III.2 Fitness

It was originally pointed out by SM sole  that $`f_i=_jJ_{ij}S_j\theta _i`$ may be taken as a measure of the “fitness” of the species $`i`$; the larger the $`f_i`$ the higher is the fitness of the corresponding species $`i`$. This definition of “fitness” is, of course, consistent with the expectation in the Darwinian scenario that, at any arbitrary stage of evolution, the species with positive fitness would not become extinct. Instead of studying the statistics of $`P(f_i)`$, or its time-average, SM investigated the relation between the species-averaged value of $`f_i`$, i.e., $`f_T=_if_i/N`$, and the size of the extinctions.

We have now computed the time-averaged probability distribution $`<P(f_i)>_t`$ (following SM, we also take $`\theta _i=0`$ for all $`i`$). From this distribution (see fig.2) it is clear that, on the average, there are very few species which do not possess adequate fitness required for survival and, therefore, become extinct. However, the distribution $`P(f_i)`$ fluctuates with time and exhibits significant negative fitnesses when extinction occurs. In these stochastic network models of ecosystems, the species become extinct because of “bad luck” rather than “bad genes” raup .

We define the incremental change of fitness of a species, say, $`i`$, per unit time step as

$$\mathrm{\Delta }f_i(t)=f_i(t)f_i(t1)$$
(3)
The distribution of $`\mathrm{\Delta }f_i`$ is shown in fig.3; it falls roughly exponentially in the tail region.

### III.3 Volatility and Coherence

Borrowing the concept of “volatility” in financial markets stanley ; bouchaud , we now define the “volatility” in the eco-system by

$$V(t)=C_1\underset{i}{}(\mathrm{\Delta }f_i)^2$$
(4)
and the time-dependent correlation function by

$$C(t)=C_2\underset{i<j}{}(\mathrm{\Delta }f_i)(\mathrm{\Delta }f_j)$$
(5)
where the normalization factors $`C_1`$ and $`C_2`$ can be chosen appropriately. For simplicity, we worked with $`C_1=1=C_2`$. We have also measured the sizes of the extinctions by the fractions $`m/N`$ where $`m`$ is the numbers of species that become extinct simultaneously at the same time step and $`N`$ is the total number of species in the system.

Large extinctions occur when volatility is also large (see fig.4). $`C`$ is a measure of “coherence” in the changes in the fitness of the species during the same period of time. Interestingly, this correlation function is very rarely negative. Even when it takes negative value, the corresponding magnitude is not large. So, most of the species do not profit from the extinction of others.

### III.4 Single common ancestor

Each node $`i`$ ($`i=1,2,\mathrm{},N`$) carries a label $`L(i)`$ that indicates the identity of its ancestor. In the initial state we assign $`L(i)=i`$. During the evolution of the system, whenever a node $`\beta `$ becomes extinct and is, therefore, re-populated by the species occupying the node $`\alpha `$ we re-label the node $`\beta `$ with $`L(\alpha )`$, i.e., the node $`\beta `$ also carries the label $`L(\alpha )`$. From our simulation we have always found that, after sufficiently long time, all the nodes carry the same label thereby indicating that all the living organisms at that stage share a single common ancestor from which all the wide variety of species have appeared as a result of massive diversification. This is perfectly consistent with theory harris  and the fossil record benton . The distributions of the sizes of phylogenic trees had already been studied by Sole and Manrubia sole .

## IV Summary

We have computed the distribution of the lifetimes of the species in SM model; the deviations from the pure power-law behavior observed for the SM model is similar to that in a recent model that describes both “macro”-evolution as well as “micro”-evolution chowstau . We have also interpreted $`f_i=_jJ_{ij}S_j`$ as the fitness of the species $`i`$ and investigated the distributions of fitness as well as the incremental changes of fitness per unit time. Finally, we have introduced the concepts of volatility and coherence; computing these for the SM model as functions of time, we have demonstrated their close relation with the sizes of extinctions.

## V Acknowledgement

This work is supported in part by the Alexander von Humboldt Foundation and by the Deutche Forschungsgeminschaft (DFG) through a joint Indo-German research project. We also thank the Supercomputer center, Jülich, for computer time on their CRAY-T3E.


text:
# Interplay of superexchange and orbital degeneracy in Cr-doped LaMnO3

## I Introduction

It was recently shown that in La1-xSrxMnO3 the realization of an insulating and ferromagnetic ground state for $`0.1<x<0.15`$ results from a superexchange (SE) driven rearrangement of the high-temperature orbital order (OO) established by the Jahn-Teller (JT) distortion.Endoh et al. (1999); Paraskevopoulos et al. (2000) Upon Sr-doping the JT distortion in the parent compound gets suppressed and the degeneracy of the Mn3+ eg electrons becomes almost restored. In this case the two-orbital model by Kugel and Khomskii becomes valid and ferromagnetic SE interactions between Mn3+ ions show up.Kugel (1982) These drive the system into an orbitally ordered and ferromagnetic state. The exact nature of the OO is still unclear, although some proposals have been made by Khomskii recently,Khomskii (2001) based on the idea of the so called orbital polaron by Kilian and Khaliullin. Kilian and Khaliullin (1999) The realization of a ferromagnetic and insulating state has been recently adressed by Khomskii and Sawatzky Khomskii and Sawatzky (1997) and also by Mizokawa et al. Mizokawa et al. (2000) For even higher Sr-doping levels ($`x0.17`$) the increasing double-exchange (DE) interactions become dominant and establish a ferromagnetic and metallic ground state.

Since the first studies of the series LaMn1-xCrxO3 in the fifties,Jonker (1956); Gilleo (1957); Bents (1957) it has been known that upon Cr-doping the antiferromagnetic Mott insulator LaMnO3 develops a ferromagnetic component. Moreover, Taguchi et al. reported (for T $`>`$ 300 K) that within the complete concentration range $`0x1`$ the samples show semiconducting behavior.Taguchi et al. (1999) Sun et al. and Zhang et al. came to similar conclusions for the temperature range T $`<`$ 300 K, recently.Sun et al. (2001); Zhang et al. (2000) Due to the ceramic nature of the samples a direct comparison of the resistivity values of these studies is vain, but both yield an increase in resistivity with increasing Cr concentration (at fixed temperatures) indicating the persistence of an insulating ground state throughout the whole concentration range. Hence, the LaMn1-xCrxO3 system is a promising candidate for studying the close relationship between JT distortion, SE interactions and orbital degeneracy, especially because in the simple ionic picture the Mn3+ (3d4) ions are partially substituted by isoelectronic Cr3+ ions (3d3), which have the same electronic configuration as Mn4+ (e.g. t$`{}_{}{}^{3}{}_{2g}{}^{}`$) in these compounds. Therefore no mobile charge carriers should be present in contrast to the Sr-doped case.

However, Sun et al. and Zhang et al. invoked the possibility of Mn3+\- O - Cr3+ double-exchange interactions in LaMn1-xCrxO3 in order to correlate their electronic transport data with the magnetic properties in their samples.Sun et al. (2001); Zhang et al. (2000) Magnetoresistance measurements (T $`>`$ 77 K) have been reported by Gundakaram et al. in the system $`Ln`$Mn1-xCrxO3 ($`Ln`$ = La,Pr,Nd,Gd) indicating the absence of DE interactions.Gundakaram et al. (1996)

The influence of Cr-doping in mixed-valence manganites has been investigated by Cabeza et al. in the system La0.7Ca0.3Mn1-xCrxO3.Cabeza et al. (1999) Their conclusion, namely that the Cr-ions do not contribute to the DE mechanism, is in accordance with the observations of Kimura et al. and Troyanchuk et al., who investigated the influence of Cr-doping in Nd0.5Ca0.5Mn1-xCrxO3 and Nd0.6Ca0.4Mn1-xCrxO3, respectively. Kimura et al. (2000, 2001); Troyanchuk (2002) In contrast, Sun et al. again found indications for DE in La0.67Ca0.33Mn1-xCrxO3.Sun et al. (2000)

In the present paper we try to shed some light on the controversial findings by comparing the influence of Cr-doping on the magnetic properties of LaMnO3 with the intensively studied DE systems La1-x(Ca,Sr)xMnO3+δ.

## II Experimental details

The polycrystalline specimens were prepared using conventional ceramic techniques. Ultra pure oxide powders (La2O3 (4N), Mn2O3 (4N) and Cr2O3 (4N), Alpha) were dried, mixed in the appropriate amounts and carefully ball-milled to ensure homogeneous samples. The LaMn1-xCrxO3 ($`0x0.5`$) samples have been pressed into pellets and were prepared by heating in flowing pure nitrogen at 1400 C for 120 hours and then slowly cooled to room temperature. This procedure was repeated four times. Powder diffraction patterns were collected employing Cu-Kα1 radiation at room temperature. The magnetic susceptibility and magnetization were measured using a dc superconducting quantum interference device (SQUID) magnetometer ($`050`$ kOe, 1.5 K $`T`$ 400 K).

ESR measurements were performed with a Bruker ELEXSYS E500 CW-spectrometer at X-band frequencies ($`\nu `$ 9.35 GHz) equipped with continuous gas-flow cryostats for He (Oxford Instruments) and N2 (Bruker) in the temperature range between 4.2 K and 680 K. The polycrystalline samples were powdered and filled into quartz tubes and fixed with either paraffin (at low temperatures 4 K $`T`$ 300 K) or NaCl (at 300 K $`T`$ 680 K).

## III Experimental Results

### III.1 Structural properties

In Figure 1 we show the orthorhombic lattice parameters (a, b, c) as a function of the Cr-concentration $`x`$. All parameters show a linear dependence on $`x`$ up to 0.3. A transition from the Jahn-Teller distorted O-phase to the pseudocubic ($`abc/\sqrt{2}`$) O-phase occurs between $`0.3<x<0.4`$. In the Jahn-Teller distorted phase the Mn-O bond lengths in the MnO6 octahedra are highly anisotropic and become isotropic in the pseudocubic O phase. Our findings are in agreement with the structural data of Gundakaram et al.Gundakaram et al. (1996) The authors also showed that the sample-preparation method described above ensures the absence of Mn4+ ions in these compounds confirmed by analysis of the samples via iodometric titration. Therefore we estimate the deviation from the ideal oxygen value to be less than 1%. The persistence of the orthorhombic structure up to $`x=0.5`$ has also been reported by Nakazono et al.Nakazono et al. (2001)

Sun et al. reported that a transition from rhombohedral to orthorhombic symmetry occurs around $`x=0.2`$.Sun et al. (2001) Since no detailed analysis of the lattice parameters and information about the atmosphere used in the preparation procedure was provided, it is difficult to argue where this discrepancy results from: A closer look on the X-ray-diffraction spectra presented by Sun et al. might reveal the persistence of the orthorhombic distortion or a mixed phase up to $`x=0.3`$. The samples in the study by Zhang et al. ($`x0.3`$) reportedly were all found to exhibit rhombohedral symmetry,Zhang et al. (2000) whereas Taguchi et al. found rhombohedral symmetry for $`x0.4`$ and orthorhombic for $`x0.4`$ and a high Mn4+ content due to oxygen non-stoichiometry. Taguchi et al. (1999)

It is known that if upon Sr-doping the percentage of Mn4+ in LaMnO3 becomes larger than 20%, the room-temperature structure is rhombohedral.Urushibara et al. (1995) The influence of oxygen content on the structural properties has been intensively studied in LaMnO3+δ by Prado et al.,Prado et al. (1999) who found a rhombohedral structure at room temperature for an oxygen excess of $`\delta >0.09`$, which corresponds to a Mn4+ content similar to a Sr concentration of $`x=0.18`$. Therefore, it is possible that oxidation (i.e. oxygen non-stoichiometry) of the samples holds for the differences in the structural, electronic and magnetic properties.

Two main reasons account for the different behavior of Cr3+ and Mn4+ ions (both are not JT active): Firstly, the ionic radius of the Cr3+ ion is larger than the corresponding one of the Mn4+ ion and hence the destabilization of the Mn3+ Jahn-Teller matrix is weaker than in the Sr-doped samples. Additionally, the mobile nature of the induced charge carriers (holes) in the case of Sr-doping allows for a dynamical change of the local electronic configuration and a homogenous distribution of the Mn4+ ions in the JT matrix, whereas the Cr3+ ions are fixed to the lattice sites at the time of sample preparation.

### III.2 Magnetization

The magnetization $`M`$ below the magnetic ordering temperatures increases upon increasing Cr content, as can be seen from the temperature dependence of the magnetization in Figure 2. This is in agreement with the results of Sun et al. and Nakazono et al.,Sun et al. (2001); Nakazono et al. (2001) who additionally observed a FC/ZFC splitting up to $`x=0.5`$ indicating a cluster-glass-like behavior. Since the magnetic transition temperatures cannot be unambiguously determined from the dc measurements, we restrict our discussion to the doping dependence of the Curie-Weiss (CW) temperature $`\mathrm{\Theta }`$. In the inset of Fig. 2 the inverse susceptibility $`vs.`$ temperature is plotted for $`x=0.1`$ and $`x=0.5`$.

The susceptibilities for all samples follow a CW law $`(T\mathrm{\Theta })^1`$ in the paramagnetic temperature regime. From this measurements we deduce the effective paramagnetic moments $`p_{\mathrm{eff}}`$ (Fig. 3a) and the paramagnetic CW temperature $`\mathrm{\Theta }`$ (Fig. 3b). From the doping dependence of $`p_{\mathrm{eff}}`$ it is noticeable that for $`x`$ values smaller than 0.1 the measured values are slightly enhanced with respect to the simple ionic spin-only values, while for $`x0.2`$ they lie well below the spin-only curve. Both parameters show a significantly different doping dependence in the O and O phases respectively. While in the Jahn-Teller distorted regime (O-phase) a linear dependence can be observed up to $`x=0.3`$, both $`p_{\mathrm{eff}}`$ and $`\mathrm{\Theta }`$ become almost constant in the O phase.

In Figure 4 hysteresis measurements up to 50 kOe at $`T=5`$ K are presented. The evolution of a ferromagnetic component upon increasing Cr-doping is clearly evident. Additionally, upon increasing Cr concentration the samples become magnetically softer, as indicated by the decreasing coercive fields in the hysteresis loops. As shown in the inset of Figure 4, the spontaneous magnetization $`M_\mathrm{s}`$ (obtained by extrapolation from the data at the highest magnetic fields, see lower panel of Fig. 4) increases linearly with the concentration of Cr3+ and reaches a maximum at $`x=0.3`$. In the O-phase $`M_\mathrm{s}`$ shows a constant behavior, similar to $`p_{\mathrm{eff}}`$ and $`\mathrm{\Theta }`$.

In principle, the behavior of the magnetic parameters discussed above is in accordance with the early studies by Jonker and Bents.Jonker (1956); Bents (1957)

### III.3 Electron Spin Resonance

Electron spin resonance detects the power $`P`$ absorbed by the sample from the transverse magnetic microwave field as a function of the static magnetic field $`H`$. The signal-to-noise ratio of the spectra is improved by recording the derivative $`dP/dH`$ with lock-in technique. ESR spectra, which are characteristic for the paramagnetic regime, are presented in Fig. 5, illustrating their evolution with Cr concentration $`x`$ (left column) and temperature $`T`$ (right column). Within the whole paramagnetic regime the spectrum consists of a broad, exchange narrowed resonance line, which is well fitted by a Lorentzian line shape as described previously.Ivanshin et al. (2000)

The integrated intensity $`I(T)`$ of the resonance line measures the spin susceptibility $`\chi _{\mathrm{ESR}}`$ of the ESR probe. For ferromagnetically coupled ions its temperature dependence usually follows a CW law, where $`\mathrm{\Theta }`$ is the CW temperature. Indeed, we observed a linear behavior of $`1/I(T)`$ in the paramagnetic regime within the whole concentration range. The CW temperatures as determined from the ESR experiments agree well with values obtained from susceptibility measurements (Fig. 3(b)). By comparison with the intensity of the reference compound Gd2BaCuO5,Goya et al. (1996) we obtained the effective moments, which are in good agreement with the magnetic susceptibility data except for $`x=0.5`$.

In the non-JT distorted orthorhombic O phase $`(x=0.4)`$ the resonance field $`H_{\mathrm{res}}`$ yields an effective g value $`g_{\mathrm{eff}}1.99`$ slightly below the free-electron value, which is characteristic for transition-metal ions with a less than half filled d-shell.Abragam and Bleaney (1970) In Fig. 6 we show $`\mathrm{\Delta }H_{\mathrm{res}}=H_{\mathrm{res}}(T)H_{\mathrm{res}}(670\mathrm{K})`$ for concentrations $`x0.15`$. At 670 K the resonance fields for these samples reach values characteristic for the O phase $`H_{\mathrm{res}}(670\mathrm{K})3.34`$ kOe. For $`x<0.15`$ the resonance shift is even larger, but the structural phase transition seems to take place beyond our accessible temperature range. However, the resonance shift at $`T=400`$ K for $`x0.15`$ (see inset of Fig. 6) clearly mirrors the structural phase transition as the shift decreases with increasing Cr concentration and drops to almost zero for $`x0.4`$.

Finally approaching the ordering temperature from above, the whole resonance becomes seriously distorted and is strongly shifted to lower fields due to the internal fields caused by the onset of magnetic order. A similar evolution of the resonance line has been reported by Sun et al. and interpreted as an indication for DE.Sun et al. (2001) However, we restrict our discussion to the paramagnetic phase ($`T>T_N`$), because in order to characterize the magnetic order accurately one needs single-crystals of defined shape instead of randomly oriented powder which gives rise to the observed distortion of the lineshape.

Figure 7 shows the characteristic temperature dependence of the ESR linewidth. All samples with $`x0.3`$ show a broad maximum in the JT distorted phase and a minimum above the critical broadening on approaching magnetic order. For the parent compound LaMnO3 the maximum cannot be detected within the investigated temperature range, but has been reported by Causa et al. Causa et al. (1999) For the highest temperatures under investigation the linewidth of all samples with $`x>0.05`$ seem to converge indicating the range of the high-temperature limit. For $`x>0.3`$ the maximum is not observable anymore and the linewidth shows a monotonous increase with increasing temperatures. The inset of Fig. 7 shows the dependence of the linewidth at 450 K on the Cr concentration. At this temperature the samples with $`x0.3`$ reveal the broad maximum, but have not yet reached the high-temperature limit. The obvious jump at the critical concentration reflects the transition from the O\- to the O-phase (Fig. 1) thus revealing a close correlation between spin relaxation and structural distortion.

## IV Discussion

### IV.1 ESR linewidth and resonance field

The temperature dependence of the linewidth and the resonance field of LaMn1-xCrxO3 for $`x0.3`$ is very similar to the ones we observed in La1-xSrxMnO3 single crystals for $`x0.15`$.Ivanshin et al. (2000) Especially for a Sr concentration of 5% the features of linewidth and resonance field can be compared to the present system. La0.95Sr0.05MnO3 is still an antiferromagnetic insulator ($`T_\mathrm{N}=140`$ K) dominated by Mn3+-O-Mn3+ SE interactions. The presence of the non-JT active Mn4+ ions, however, weakens the cooperative JT matrix resulting in a reduced JT transition temperature $`T_{\mathrm{JT}}=600`$ K. Below this temperature both $`\mathrm{\Delta }H`$ and $`H_{\mathrm{res}}`$ show an anisotropic behavior, which can be regarded as a characteristic feature of the O phase:Ivanshin et al. (2000); Alejandro et al. (2001)

If the external field is applied within the ferromagnetically coupled ab-plane, the linewidth reveals a broad maximum within the JT distorted O phase, whereas $`\mathrm{\Delta }H`$ is nearly constant, if the external field is parallel to the antiferromagnetically coupled c-axis. We also investigated a powdered La0.95Sr0.05MnO3 sample and found that the linewidth shows the same features as the single crystal for the external field applied within the $`ab`$-plane. Correspondingly, we observed a resonance shift to higher magnetic fields within the JT distorted O phase, when the external field is within the $`ab`$-plane, and a constant behavior along the c-axis. In our recent approach we presented a detailed analysis for the angular dependence of the linewidth and the resonance field in La0.95Sr0.05MnO3 in the O phase by taking into account the crystal field (CF) and the Dzyaloshinsky-Moriya (DM) interaction,Deisenhofer et al. (2002) which were shown to yield by far the largest contribution to the linewidth in manganites by Huber et al.Huber et al. (1999) Due to the similarities in the temperature dependence of $`\mathrm{\Delta }H`$ and $`H_{\mathrm{res}}`$ we assume that these two interactions also account for the ESR properties in Cr-doped LaMnO3.

Unfortunately, in the case of polycrystalline samples such a detailed analysis is not possible due to the random orientation of the powder. However, a similar behavior of $`\mathrm{\Delta }H`$ vs. $`T`$ has also been reported in polycrystalline LaMnO3+δ by Tovar et al.,Tovar et al. (1999) who found a broad linewidth maximum in the O phase up to $`\delta =0.04`$. The corresponding JT transition temperatures decrease from 700 K for $`\delta =0`$ to 400 K for $`\delta =0.05`$, and for $`\delta >0.05`$ the cooperative JT distortion vanishes concomitantly with the broad maximum. The authors describe the temperature dependence of the exchange-narrowed linewidth by following Huber et al. as

$$\mathrm{\Delta }H(T)=\frac{\chi _0(T)}{\chi (T)}\mathrm{\Delta }H_{\mathrm{}}$$
(1)
with the free Curie susceptibility $`\chi _0T^1`$ and the static susceptibility $`\chi (T)(T\mathrm{\Theta })^1`$.Huber et al. (1999) The high-temperature linewidth $`\mathrm{\Delta }H_{\mathrm{}}`$ depends on the contributions of both the CF and the DM interaction. Thus, we fitted the linewidth data for all samples under investigation by using eq. 1 (see Fig. 7), which can neither describe the critical broadening of the linewidth in the vicinity of a magnetic transition nor the structural transition at $`T_{\mathrm{JT}}`$, but fits very nicely in between these transitions. The obtained parameters $`\mathrm{\Theta }`$ and $`\mathrm{\Delta }H_{\mathrm{}}`$ are shown in Fig. 8. Though the obtained CW temperatures are, in comparison with the values in Fig. 3b, slightly lower for $`x0.3`$ (O phase) and enhanced for $`x>0.3`$ (O phase), the overall tendency is the same. Moreover $`\mathrm{\Delta }H_{\mathrm{}}`$ shows the same behavior as the linewidth values taken at 450 K (see inset of Fig. 7), confirming that the maximum linewidth is a good measure for the strength of the interactions.

The difference in $`\mathrm{\Delta }H_{\mathrm{}}`$ observed between the O and the O-phase is due to the fact that concomitantly with the disappearance of the cooperative JT effect the contribution of the single-ion anisotropy of the CF vanishes at $`T_{\mathrm{JT}}`$ whereas the DM contribution is due to the tilting of the octahedra, which for LaMnO3 has been shown to remain nearly unchanged through the transition.Huang et al. (1997); Zimmermann et al. (1998)

The obtained values in both phases are in agreement with those reported for La1-xCaxMnO3 and LaMnO3+δ of $`\mathrm{\Delta }H_{\mathrm{}}2.7`$ kOe for the O-phase and $`\mathrm{\Delta }H_{\mathrm{}}2.2`$ kOe for the O phase yielding values of the order of 1 K for both interaction.Huber et al. (1999); Tovar et al. (1999)

The influence of DE on the ESR linewidth in the paramagnetic regime is not fully understood and still subject of controversial discussion: Shengelaya et al. found that in La0.8Ca0.2MnO3, where the cooperative JT distortion is already suppressed, the temperature dependence of the linewidth correlates with the conductivity and both parameters can be described by a small-polaron hopping model $`1/T\mathrm{exp}(E_\mathrm{a}/k_\mathrm{B}T)`$ with similar activation energies $`E_\mathrm{a}`$.Shengelaya et al. (2000) In contrast, in the optimally doped region $`x=0.33`$ the linewidth has been described previously by eq. (1) only.Causa et al. (1998); Oseroff et al. (1996) Some authors argue that due to the different time scales the DE interaction, which manifests itself in the paramagnetic phase through the hopping of the JT polaron, cannot have any effect on the spin relaxation.Huber et al. (1999); Rivadulla et al. (2000) However, we want to emphasize that in our case the linewidth can be satisfactorily described without any additional relaxation process connected to polaron hopping.

### IV.2 Magnetic Properties and ESR intensity

Though all observed parameters show a close correlation to the structural changes induced upon Cr-doping, the origin of the observed ferromagnetic component of the magnetization cannot easily be explained. From the values of the high-field magnetization (5 T) we deduce that the magnetic phase is not purely ferromagnetic, since the deviation from the theoretical values of $`34\mu _\mathrm{B}/\mathrm{f}.\mathrm{u}.`$ for the full magnetic moment of all Mn3+ and Cr3+ spins cannot be attributed to temperature induced disorder, only. Therefore, we have to consider all exchange interactions between the Mn3+ and Cr3+ ions that could contribute to the magnetic state of the system:

In the parent compound LaMnO3 the JT distortion lifts the degeneracy of the $`e_g`$ orbitals and the Mn3+\- O - Mn3+ SE interaction yields antiferromagnetic coupling along the c-axis and ferromagnetic coupling in the ab-plane. In the case of degenerated $`e_g`$ orbitals, the Kanamori-Goodenough rules Goodenough (1963) do not apply and the two-orbital Kugel-Khomskii model yields a purely ferromagnetic Mn3+\- O - Mn3+ SE interaction.Kugel (1982) Due to the weakening of the JT distortion upon Cr-doping, the degeneracy of the $`e_g`$ orbitals can become almost restored (similar to the ferromagnetic insulating state of La1-xSrxMnO3) and favor ferromagnetism (FM). Taking into account the Kanamori-Goodenough rules Goodenough (1963) the Mn3+\- O - Cr3+ superexchange is ferromagnetic and therefore can also account for the increasing values of $`\mathrm{\Theta }`$ in the O-phase with increasing Cr-doping. Correspondingly, the antiferromagnetic Cr3+\- O - Cr3+ superexchange should contribute considerably with increased Cr content and could account for the stabilization of the magnetic state for $`x>0.3`$, where the observed parameters become almost constant. Finally, a DE interaction between Mn3+\- O - Cr3+ ions has been proposed.Sun et al. (2001); Zhang et al. (2000)

Regarding the existence of DE, we can compare our data to observations of Prado et al. Prado et al. (1999) in LaMnO3+δ, where one can clearly see the impact of Mn4+ \- O - Mn3+ double exchange on the magnetic properties: Though the hysteresis loops show a similar behavior upon increasing oxygen surplus as upon Cr-doping, the magnetic phase transition from AFM to FM at $`\delta 0.04`$ is characterized by the appearance of the full magnetic moment of $`34\mu _\mathrm{B}/\mathrm{f}.\mathrm{u}.`$ These features are also found upon doping with Sr.Paraskevopoulos et al. (2000) How difficult it is to distinguish the involved interactions from the dc magnetization data only can be seen in comparison with the behavior of $`M`$ vs. $`x`$ in electron doped Ca1-xLaxMnO3,Neumeier and Cohn (200) which looks very similar to our system. However, this behavior has been found to correlate with the electron mobility and attributed to the dynamical charge transfer due to DE interactions. As neither such correlations nor the appearance of the full magnetic moment has been observed in LaMn1-xCrxO3, we will consider the influence of DE in the following as negligible.

Due to the delicate interplay of the above SE interactions, Jonker favored the existence of a complex ferrimagnetism (FiM),Jonker (1956) and Bents found in his neutron diffraction study that not only FM, but also A-type AFM (for $`x0.15`$) and G-type AFM (for $`x>0.3`$) are present in the system.Bents (1957) A theoretical density-functional study by Yang et al. based on the data of these early publications supports the idea of FiM on account of positive Mn3+\- O - Mn3+ and negative Cr3+\- O - Cr3+ superexchange interactions. However, positive Mn3+\- O - Cr3+ interaction has not been found in the calculations.Yang et al. (2000) Very recently, Ono et al. suggested another scenario based on their soft X-ray magnetic circular dichroism experiments:Ono et al. (2000, 2001) They found that the Mn ions are aligned parallel to the direction of the magnetization, whereas the Cr magnetization seems to disappear. In order to simulate such a behavior by using the Monte Carlo method, they had to assume a negative Mn3+\- O - Cr3+ interaction in contrast to the Kanamori-Goodenough rules. Assuming antiferromagnetic Mn3+\- O - Cr3+ superexchange the latter authors explain the cluster-glass-like behavior of the system as a frustration of the spins due to the competing SE interaction.Nakazono et al. (2001) Although the reasons for the disappearance of the Cr magnetization are rather unclear, such a scenario is in accordance with the observed values of the spontanous magnetization $`M_\mathrm{S}`$ values (see inset of Fig. 4): The linear increase in the JT distorted O phase can be attributed to the partially restored degeneracy of the $`e_g`$ orbitals,Kugel (1982) weakening the antiferromagnetic Mn3+\- O - Mn3+ SE (along the c-axis) and finally yielding only ferromagnetic Mn3+\- O - Mn3+ interaction for $`x>0.3`$. Thus, the theoretical values of 1.6 -2 $`\mu _\mathrm{B}/\mathrm{f}.\mathrm{u}.`$ are in agreement with our data.

Astonishingly, we also observe that in contrast to Ca/Sr-doped LaMnO3, where the ESR signal has been attributed to both Mn3+ and Mn4+ ions, Causa et al. (1998); Ivanshin et al. (2000) the ESR signal in the paramagnetic regime seems to be originated by the Mn3+ ions only (see Fig. 3): The effective moments determined by ESR follow the theoretical curve for Mn3+ up to $`x=0.5`$ very nicely. The corresponding values from the magnetization measurements, however, are slightly enhanced in comparison to the ESR data, but the overall tendency is similar. In contrast to the SQUID measurements (see inset of Fig. 2), the ESR experiments were performed up to 670 K allowing a more reliable determination of the effective moments. Thus, the frustration of the Cr3+ ions even seems to influence the paramagnetic phase.

## V Conclusions

We presented a detailed magnetic analysis of the doping effect of Cr on LaMnO3. The parent compound is an A-type AFM with orbital order induced by a strong JT effect due to the double degeneracy of the $`e_g`$-orbitals of the Mn3+ ions. Upon doping with Cr3+ ions, ferromagnetic interactions show up. The magnetically ordered state can be understood in terms of a competition of positive Mn3+-O-Mn3+ and negative Mn3+-O-Cr3+ and Cr3+-O-Cr3+ SE interactions. Assuming the disappearance of the Cr magnetization the increasing ferromagnetic component can be ascribed to the Mn3+ ions the partial restoration of the double degeneracy of the $`e_g`$-orbitals due to the weakening of the JT matrix by non-JT active Cr3+-ions. Additionally, the Jahn-Teller driven orbital order becomes weaker and finally disappears for doping levels between 0.3 and 0.4. All parameters, $`i.e.`$ $`\mathrm{\Theta }`$, $`H_{res}`$ and $`\mathrm{\Delta }H`$, reflect this transition in their doping dependence. Thus, the system can be considered as a ferromagnetic insulator with respect to Mn3+ in the sense of the two-orbital Kugel-Khomskii model, as the competing negative SE interactions seem to produce a disappearance of the Cr magnetization. The ESR properties can be naturally explained by the structural disorder induced upon Cr-doping and the effects of the CF and DM interactions. Moreover, the Cr3+ ions seem not to contribute to the ESR signal. From our experimental data we do not see any need to invoke the possibility of DE interaction between Cr and Mn ions. The questions about the origin of the disappearance of the Cr magnetization and why the Cr3+ ions seem to be ESR-silent remain unsolved problems and challenge further experimental and theoretical investigations.

###### Acknowledgements.

We want to thank Ch. Hartinger and A. Spushk for fruitful discussions and the authors of Ref. 20 for providing us with a preprint of their work. This research was supported by the BMBF via the contract number VDI/EKM 13N6917 and partly by the Deutsche Forschungsgemeinschaft via SFB 484 (Augsburg).


text:
# Melting of sodium clusters

## I Introduction

The study of the thermal stability and melting transition of sodium clusters is nowadays a very active field of research. This interest is motivated, in part, by the fact that just recently it was possible to make direct comparisons between theoretical results, mainly extracted from computer simulations , and experimental data  obtained for the melting-like transition of sodium clusters of different sizes. In fact, during the last few years the caloric curve, heat capacity, and melting temperatures of NaN, $`N`$ = 55-200, clusters have been measured using the temperature dependence of the fragmentation mass spectra . The melting points of these clusters were found to be on average 33 % (120 K) lower than the bulk value, and more surprisingly, large variations were observed in the melting temperatures ($`\pm `$ 30 K) with changing cluster size .

On the theoretical side, molecular dynamics (MD) and Monte Carlo (MC) methods have been used to provide microscopic descriptions on the melting mechanisms of sodium clusters in the size range of 8-147 atoms. In these studies, the metallic bonding of the sodium clusters has been described using different levels of approximation, from empirical many-body model potentials  and tight-binding hamiltonians  to first-principles density functional theory .

Despite the large amount of information obtained from the above theoretical studies, several questions on, for example, the irregular variations of the melting temperatures with respect to the cluster size remain unsolved. One difficulty existing in the theoretical approaches to study the melting-like transition of sodium clusters is related with the interplay between the geometric and electronic structure effects in such systems . Although MD simulations that use density functional theory explicitly take into account the electronic degrees of freedom and their effects in the cluster geometry, the limitation of this approach is related to the relatively short simulation times (a few picoseconds) during which time-averages of the structural and dynamical properties are calculated. This problem, caused by the extremely large computational effort involved in the first-principles MD simulations, is especially critical in the study of the melting transition where large thermal fluctuations are present, and therefore much longer simulation times are required to calculate converged time averages.

On the other hand, MD simulations based on many-body model potentials allow extension of the simulation time up to the nanosecond regime, employing reasonable computational resources. However, in this case the description of the metallic bonding does not explicitly include the electronic degrees of freedom. In the present work, we adopt the latter approach to perform molecular dynamics simulation of the melting-like transition of NaN, $`N`$ = 13, 20, 55, 135, 142, and 147, clusters using a many-body Gupta potential , and simulation times of $``$ 50 nanoseconds. The objectives of this study are: (1) to test the phenomenological many-body Gupta potential in the description of the melting transition using adequate simulation times and (2) to compare the predictions of this simulation with those obtained from the same many-body potential but using the MC method  for the averaging procedure, and also with the results obtained from first-principles MD  using shorter simulation times. These tests and comparisons will provide additional insights on the melting mechanisms in sodium clusters and provide useful information on the performance of the different simulation methods.

In Section II we describe the model potential and provide details on the simulation method. Results on the caloric curves, heat capacities and thermal stability properties of different cluster sizes are given in Sec. III. A summary of this work is presented in Sec. IV.

## II Metallic potential and simulation procedure

The many-body model potential used in this work corresponds to the Gupta potential  that is based on the second moment approximation of a tight-binding hamiltonian . Its analytical expression is given by:

$$V=\underset{i=1}{\overset{N}{}}V_i$$
(1)
$$V_i=A\underset{ji}{}e^{p\left(\frac{r_{ij}}{r_0}1\right)}\xi \left(\underset{ji}{}e^{2q\left(\frac{r_{ij}}{r_0}1\right)}\right)^{\frac{1}{2}}$$
(2)
where $`r_0`$, $`A`$, $`\xi `$, $`p`$, and $`q`$ are adjustable parameters . For sodium clusters these parameters have been fitted to band structure calculations . Their values are: $`A`$=0.01595 eV, $`\xi `$=0.29113 eV, $`r_0`$=6.99 bohr, $`p`$=10.13, and $`q`$=1.30 . This type of potential has been extensively utilized in other metal cluster simulations , obtaining good agreement with results generated from first-principles methods.

In order to study the cluster melting-like transition, we use the constant-energy MD method to calculate the structural changes as a function of the cluster energy. Within this method, Newton’s equations of motion are solved for each atom within the cluster using the Verlet algorithm . Through this procedure, one obtains the atomic positions and momenta as a function of time, that are used to calculate time-averages of physical quantities characterizing the cluster structure and dynamics. A typical calculation consists in heating up a cluster from its zero temperature configuration until it transforms into a liquid-like cluster. To simulate this procedure the cluster total energy is increased in a step-like manner by scaling up the atomic velocities and therefore augmenting the kinetic energy of the cluster.

The simulation starts by slightly perturbing the coordinates corresponding to the lowest-energy configuration of the cluster and setting up the atomic momenta to a zero value, in order to eliminate the translation of the cluster center of mass and its rotational motion. The time step of the MD runs is 2.4 fs, which provides total energy conservation within 0.001 %. For each initial condition the cluster is equilibrated during 104 time steps and the time averages are calculated using 107 time steps in the solid and liquid regions, but 2 times longer simulation times are used in the melting-like transition region. This amount of averaging time is a necessary condition to obtain converged physical quantities characterizing the melting process in these finite systems.

To obtain the lowest-energy structure of each cluster size we combine simulated quenching techniques  and evolutive algorithms , which are able to perform global searches of the potential energy surface (PES) in a very efficient way despite the large number of degrees of freedom involved in the cluster structure optimization. These methods do not only provide the lowest-energy configuration but the distribution of isomers in the low energy range. This optimization procedure has been successfully used in other cluster structure optimizations .

The behavior of the structural and thermal properties during the heating process of the cluster is monitored by calculating the temperature, heat capacity, and root-mean-square (rms) bond length fluctuations as a function of the cluster total energy using the following expressions:

$$T=\frac{2<E_k>}{(3N6)k_B}$$
(3)
$$\frac{C}{Nk_B}=\left[NN(1\frac{2}{3N6})<E_k><E_k^1>\right]^1$$
(4)
$$\delta =\frac{2}{N(N1)}\underset{i<j}{\overset{N}{}}\frac{[<r_{ij}^2><r_{ij}>^2]^{1/2}}{<r_{ij}>}$$
(5)
where $`E_k`$ is the cluster kinetic energy, $`k_B`$ is the Boltzmann constant, $`<`$$`>`$ denotes a time average, and $`r_{ij}`$ corresponds to the distance between atoms $`i`$ and $`j`$. The above mathematical expressions were introduced in Refs.  in order to calculate structural and thermal properties of clusters from computer simulation studies.

## III Results and Discussion

### A Na13

The lowest-energy configuration of the Na13 cluster corresponds to an icosahedral structure shown in Fig. 1(a). This geometry was used to initiate the heating up procedure through the MD method. The cluster temperature (caloric curve) and specific heat as a function of the cluster total energy, calculated using Eqs. (3) and (4), are displayed in Figs. 2(a) and 3(a), respectively. The change in the slope of the caloric curve as well as the existence of a maximum in the specific heat are characteristics of the solid-to-liquid transition in clusters . These features are clearly seen in Figs. 2(a) and 3(a), indicating a melting-like transition in the Na13 cluster. This transition, which occurs over a broad range of energy involving one or more intermediate stages (such as isomerizations, coexistence of liquidlike and solidlike phases, partial and surface melting etc.), has been widely discussed in earlier studies of atomic clusters. Figure 4(a) shows the rms bond-length fluctuation, $`\delta `$, as a function of the cluster total energy, calculated using Eq. (5). It shows two abrupt variations at different energies that are in contrast with the analog of the Lindemann criterion  for bulk melting, where a single abrupt increase in $`\delta `$ is observed.

By performing thermal quenchings from configurations generated at different total energies during the MD simulation, the different melting mechanisms occurring in the Na13 cluster can be investigated. It is found that the first abrupt change in $`\delta `$ at low energy (temperature) is due to cluster isomerization involving only surface atoms, whereas the second increase at higher energy (temperature) corresponds to isomerizations where the central (bulk like) atom as well as the surface atoms are involved. The onset of the surface and volume isomerizations occur at $`T`$=149 and 226 K, respectively, whereas the temperature corresponding to the maximum of the specific heat (which indicates the transition to the liquidlike state) is at $`T`$=260 K. Similar isomerization mechanisms and transition temperature values have been obtained from MC simulations of Na13, using the same many-body model potential . A similar melting behavior has also been obtained for Ni13 and Al13 through MD simulations using the Gupta potential . Nevertheless, this melting mechanism is not exclusive of metal clusters interacting by a many-body Gupta potential since both surface and volume isomerizations have also been obtained for a 13-atom cluster modeled by a pair-wise additive potential with a soft repulsive core interaction .

### B Na20

The double icosahedron with a capped atom over the central pentagonal ring displayed in Fig. 1(b) corresponds to the lowest-energy structure of Na20. By heating up this cluster we obtain the caloric curve and specific heat displayed in Figs. 2(b) and 3(b), and the rms bond length fluctuation shown in Fig. 4(b). In addition to the slight change in the slope of the caloric curve at low energies, the specific heat shows a shoulder at low energy (temperature), before it reaches its maximum at higher energy (temperature). A very abrupt increase in the $`\delta `$ value is obtained at low energy that corresponds, according to the caloric curve, to a cluster temperature of 57 K. A second abrupt increase in $`\delta `$ at around 157 K is also obtained after further heating of Na20.

The microscopic features characterizing this melting behavior can be extracted by performing periodic quenchings using the cluster configurations along the MD run at a given cluster energy. The analysis of the quenched cluster structures indicate that the first abrupt increase in $`\delta `$, that is in correspondence with the shoulder in the specific heat at low energies, is related with isomerization transitions where the extra atom incorporates into the cluster surface, generating a fluctuating six-atom ring in the central region of the cluster structure. This isomerization is equivalent to a surface reconstruction of the icosahedral surface without diffusion. The second jump in $`\delta `$ is associated to a surface melting stage where the double icosahedron transforms into a more compact structure containing two internal atoms. At higher energies (temperatures) complete melting is observed, characterized by diffusive motion where the $`\delta `$ value levels off. The two abrupt increases in $`\delta `$ obtained for Na20 are equivalent to those described for Na13, except that in the larger cluster the first increase corresponds to an isomerization due to the outer atom lying over the cluster surface. These phenomena occurring in Na13 and Na20 before they fully melt, showing atomic diffusion, can be denominated as a premelting stage.

The premelting phenomenon was first obtained in MD simulations of Ni14 and Ni20  and some BeN, $`N`$=9,11,12,14, clusters . Later, other metal clusters with magic and non-magic number sizes also reveled a premelting stage (see, for example, Refs. ). The two-step melting transition described above was also obtained using the q-jumping MC method and the Gupta potential . A good agreement was found in the premelting and melting temperatures calculated with this approach . However, higher temperature values, for the two transitions, were found by first-principles orbital-free MD on Na20 , which might be due to the much shorter simulation time used in such calculation.

### C Na55

Figure 1(c) shows the lowest-energy configuration of Na55 corresponding to the two-shell Mackay icosahedron. The caloric curve and specific heat obtained by heating up this structure are displayed in Figs. 2(c) and 3(c), respectively. Figure 4(c) shows the rms-bond length fluctuation where an abrupt increase is observed at a total energy that, according to the caloric curve, indicates a melting temperature of 151 K, whereas the maximum in the specific heat is obtained at a slightly higher temperature of 166 K. In contrast to Na13 and Na20, the Na55 cluster shows a single abrupt change in $`\delta `$. However, by visualizing the dynamical trajectories of each atom in the cluster at the temperatures within the transition region, we found that the melting process develops in several stages. In the first one, the most external layer fluctuates between an incomplete and complete icosahedral surface by expelling and receiving atoms back and forth. At higher energies an exchange of atoms between the intermediate and most external layers is obseved, and with a further increase in energy, fully melting is found, where the central atom also contributes to the cluster diffusive motion. This complex melting mechanism has been recently studied in Al55 and other aluminum clusters by introducing dynamical degrees of freedom .

Similar melting stages have also been obtained by MC simulation of Na55 . However, in those calculations a slightly higher melting temperature of 175 K was reported. The first-principles orbital-free MD simulation of Na55 also reported  a melting transition at 190 K. Again, the smaller melting temperature obtained in the present work might be due to the much longer simulation times we have used in our MD simulations as compared to those used in Ref.. Nevertheless, none of the melting temperatures calculated by us and other authors  are in agreement with the experimental value (320 K) reported for Na55 . This discrepancy between theory and experiment has not yet been solved since it would require a more detailed modeling of the energy landscape of Na55, that includes not only information on the basins of attraction of the equilibrium structures but also on the topology around saddle points connecting the lowest-energy minima. Additionally, in a more detailed description of the potential energy surface it would be possible that the global minimum for Na55 does not correspond to the icosahedral structure but to an unknown special geometry with larger thermal stabiltity against melting. At present, it represents a theoretical challenge to characterize the potential energy surface of systems with such number of degrees of fredom using first-principles methods, however intense efforts are currently being performed to solve this problem. On the other hand, further experimental work is expected in the near future that confirms the relatively high value of the melting temperature of the Na55 cluster .

### D Na135, Na142, and Na147

The global minima of the larger sodium clusters investigated in this work correspond to icosahedral structures. The three-layer Mackay icosahedron shown in Fig. 1(f), is the lowest-energy isomer of Na147. The lowest-energy structures of the Na135, Na142 are incomplete icosahedra obtained by removing 12 and 5 vertex atoms, respectively, from the 147 Mackay icosahedron (see Figs. 1(d) and 1(e), respectively). Despite the existence of an incomplete surface layer in the lowest-energy structure of Na135, Na142, they show caloric curves, Figs. 2(d-e), and rms bond-length fluctuations, Figs. 4(d-e), very similar to those obtained by heating up the complete icosahedron structure of the 147-atom cluster, Figs. 2(f) and 4(f). The calculated melting temperatures obtained using the Lindemann criterion ($`\delta `$ $``$ 0.15) for Na135, Na142, and Na147 are 135 K, 190 K, and 171 K, respectively. These values are smaller than those obtained from the maximum of the specific heat (see Table I). By visualizing the atomic coordinates as a function of time at the energies where the $`\delta `$ values change abruptly, it is observed that the melting is initiated at the cluster surface. For the three sizes investigated, the atomic mobility increases with temperature starting from the most external layer and propagating into the internal layers. This stage, known as surface melting , precedes the complete cluster melting characterized by the diffusive motion of all the atoms in the cluster which is observed at temperatures where the specific heat is maximum.

Our calculated melting temperatures for these larger sodium clusters are lower than those obtained by other authors using the MC method with the same potential  and the orbital free MD simulations . This difference, as in the smaller cluster sizes, we assume is due to the much longer simulation times we have used in our MD calculations. On the other hand, although our results show that the largest melting temperature corresponds to the Na142 cluster, in agreement with the experimental data , the absolute values of our calculated melting temperatures are about 30 $`\%`$ lower than the experimental values . These results indicate that the many-body Gupta potential, which does not include electronic degrees of freedom, only provide a qualitative description of the melting of sodium clusters.

## IV Summary

The melting-like transition of NaN, N = 13, 20, 55, 135, 142, and 147, has been investigated through microcanonical MD simulations using a phenomenological many-body Gupta potential. The solid-to-liquid transition was studied by calculating caloric curves, rms bond-length-fluctuations and specific heats. The indicators of the cluster melting correspond to changes in the slope of the caloric curve, abrupt increases in the $`\delta `$ values, and the existence of maxima in the specific heats. Table I shows the melting temperatures calculated for all cluster sizes using those criteria. The main features coming out from these data are: (i ) The melting temperatures calculated from the maxima of the specific heats are systematically higher than the values obtained using the Lindemann criterion. (ii ) There is an irregular variation in the melting temperatures as a function of the cluster size, the highest value being the one corresponding to the Na142 cluster. These results are in qualitative agreement with the experimental data . (iii ) The calculated melting temperature for the Na55 cluster is about 40 $`\%`$ lower than the reported experimental value . (iv ) The melting transition in sodium clusters is a complex phenomenon that involves several stages in which the system undergoes different structural changes (isomerization, premelting and surface melting) before it shows a diffusive regime characteristic of the liquidlike phase.

A comparison of the present results with those obtained using the same many-body potential and the MC method , and with those generated from orbital-free MD simulations , indicate that in general, the melting temperatures calculated by heating up the lowest-energy isomer, are lower when much longer simulation time is employed. In this work, the simulation time was extended up to the nanosecond time regime where it is very likely that the time-averages of the physical quantities that characterize the cluster melting might be much better converged.

In order to obtain a better (quantitative) agreement with the experimental results on the melting of sodium clusters it would be necessary to either extend the simulation time in the first-principles MD calculations or to design a many-body potential that describes with a higher level of approximation the complex topology of the potential energy landscape of sodium clusters. Work in both directions is currently under progress.

###### Acknowledgements.

This work was supported by DGSCA-UNAM Supercomputing Center. JARN acknowledges a graduate fellowship from DGEP-UNAM.