format: text

text:
# REFERENCES

Magnetic transition and polaron crossover in a two-site single polaron model including double exchange interaction

Jayita Chatterjee*** e-mail: moon@cmp.saha.ernet.in, Manidipa Mitra and A. N. Das

Saha Institute of Nuclear Physics
1/AF Bidhannagar, Calcutta 700064, India

PACS No.71.38. +i, 63.20.kr, 75.30. Vn

Abstract

A two-site double exchange model with a single polaron is studied using a perturbation expansion based on the modified Lang-Firsov transformation. The antiferromagnetic to ferromagnetic transition and the crossover from small to large polaron are investigated for different values of the antiferromagnetic interaction ($`J`$) between the core spins and the hopping ($`t`$) of the itinerant electron. Effect of the external magnetic field on the small to large polaron crossover and on the polaronic kinetic energy are studied. When the magnetic transition and the small to large polaron crossover coincide for some suitable range of $`J/t`$, the magnetic field has very pronounced effect on the dynamics of polarons.

1. Introduction

The origin of ferromagnetism in manganese perovskites $`La_{1x}X_xMnO_3`$ ($`X=Ba`$, $`Sr`$, $`Ca`$ etc.) is the double exchange mechanism . The discovery of anomalous magnetotransport phenomena  in these compounds has stimulated intensive studies in its magnetic as well as electrical properties. However, the simple double-exchange interaction alone is not sufficient to explain the experimental results . Depending on the doping, temperature and the radius of the dopant ion these oxides show various phases . A complete understanding of the properties of $`Mn`$ perovskites is still far from clear.

One of the keys to understand unusual physical properties is to find the role of coupling between the carriers and the underlying lattice. Several theoretical models have been proposed based on lattice-carrier coupling . Many recent experiments  indicate that the electron-phonon ($`e`$-ph) interaction shapes its properties very crucially. Moreover, small to large polaron crossover is reported by many experimental groups . There are models  which incorporate double exchange interaction in a polaronic model. Min and co-workers  have studied the role of $`e`$-ph interaction in systems where double exchange interaction is present. The combined model of spin double exchange and lattice polaron  is used to investigate the effect of small to large polaron crossover on the magnetic and transport properties under the mean-field approximation scheme. It concludes that the effect of polaron narrowing on the colossal magnetoresistance (CMR) is more pronounced in CMR manganites with low magnetic transition temperature ($`T_c`$) than in the high $`T_c`$ manganites. Within the mean field theory  the magnetic transition, metal-insulator transition and a large drop in magnetoresistance occur at the same temperature.

In an attempt to have a clearer view from a nearly exact calculation, we include double exchange interaction in a two site one polaron model and follow a perturbation expansion  based on a modified Lang-Firsov (MLF) phonon basis where the lattice distortions produced by the electron are treated as variational parameters . This method shows good convergence as well as nearly exact results  for almost the entire range of $`e`$-ph coupling for $`t/\omega _01`$, where $`\omega _0`$ is the phonon frequency. We investigate the ferromagnetic (FM) to antiferromagnetic (AFM) transition, the crossover from large to small polaron, the behavior of the effective hopping of the itinerant electron and the kinetic energy as a function of the $`e`$-ph coupling strength for the ground state of the system. Effect of the magnetic field on the large to small polaron crossover and on the polaronic kinetic energy are also studied.

The paper is organized as follows. In Sec. 2 we define the model hamiltonian describing different interactions and calculate different physical quantities which indicate the behaviour of small-large polaron crossover and AFM-FM transition as a function of $`e`$-ph coupling. In Sec. 3 we present the results obtained in our calculation and discussions. Sec. 4 contains the concluding remarks.

2. Formalism

The Holstein model  including double exchange interaction for the two-site single-polaron system is described by the Hamiltonian

$`H`$ $`=`$ $`{\displaystyle \underset{i,\sigma }{}}ϵn_{i\sigma }{\displaystyle \underset{\sigma }{}}t\mathrm{cos}({\displaystyle \frac{\theta }{2}})(c_{1\sigma }^{}c_{2\sigma }+c_{2\sigma }^{}c_{1\sigma })`$ (1)
$`+`$ $`g\omega _0{\displaystyle \underset{i,\sigma }{}}n_{i\sigma }(b_i+b_i^{})+\omega _0{\displaystyle \underset{i}{}}b_i^{}b_i+J{\displaystyle \underset{<ij>}{}}\stackrel{}{S}_i.\stackrel{}{S}_j`$ (2)
where $`i`$ =1 or 2, denotes the site. $`c_{i\sigma }`$ ($`c_{i\sigma }^{}`$) is the annihilation (creation) operator for the itinerant electron with spin $`\sigma `$ at site $`i`$, $`n_{i\sigma }`$ (=$`c_{i\sigma }^{}c_{i\sigma }`$) is the corresponding number operator and $`g`$ denotes the on-site $`e`$-ph coupling strength. $`b_i`$ and $`b_i^{}`$ are the annihilation and creation operators respectively for the phonons corresponding to interatomic vibrations at site $`i`$, $`\omega _0`$ is the corresponding phonon frequency and $`ϵ`$ is the bare-site energy. $`\stackrel{}{S}_i`$ and $`\stackrel{}{S}_j`$ are the core-spins at the site $`i`$ and $`j`$ respectively, $`\theta `$ is the angle between the core-spins $`\stackrel{}{S}_i`$ and $`\stackrel{}{S}_j`$. The transfer hopping integral ($`t`$) is modified through the relative angle $`\theta `$ as $`t\mathrm{cos}(\frac{\theta }{2})`$ because of the strong Hund’s coupling between the spins of the core eletrons and itinerant electron . $`J`$ is the superexchange antiferromagnetic interaction between the neighbouring core-spins. Since we will restrict ourselves to the single-electron case we will not consider the electron spin indices.

If one compares our simple model with the manganite (e.g. $`La_{1x}Ca_xMnO_3`$) system the core spin and the itinerant electron may be identified with the $`t_{2g}^3`$ (localized) electrons and the $`e_g^1`$ (mobile) electron of $`Mn^{3+}`$ of the manganite system respectively. It may be mentioned that for manganites the Jahn-Teller(JT) coupling is important whereas the Holstein phonon mode, in general, should correspond to the breathing mode in manganites . To include JT coupling one should consider two $`e_g`$ orbitals (say $`\alpha `$ and $`\beta `$) and their coupling to the JT phonon mode. For the single electron case, as considered here, if one neglects the interorbital hopping ($`t_{\alpha \beta }=0`$) then the Hamiltonian can be separated into two parts each corresponds to 1-orbital case and as a consequence the 2-orbital problem is reduced to an effective single orbital problem with an effective Holstein type interaction as in Hamiltonian (1). For the two electron case the above simplification cannot be done and a study of two-site two-electron problem considering two orbitals in the context of manganites is in progress. We have considered here one itinerant electron in a two-site 1-orbital system. So, the density of itinerant electron $`<n>`$ is equal to 0.5 which, in principle,corresponds to $`x=0.5`$ for manganite system. It may be mentioned that the studies of 1-orbital  and 2-orbital Kondo model  in the context of manganites using Monte Carlo techniques showed that the results for 2-orbital model at $`n=1`$ ($`x=0`$) are similar to those for single orbital model at $`x=0.5`$.

By introducing new phonon operators $`a=(b_1+b_2)/\sqrt{2}`$ and $`d=(b_1b_2)/\sqrt{2}`$, the Hamiltonian ($`H`$) is separated into two parts : one corresponding to the in-phase mode which does not couple with the electronic degrees of freedom and the other involving out-of-phase mode ($`H_d`$) which represents an effective $`e`$-ph system and cannot be solved analytically .

The MLF transformation with variable phonon basis is used so that a convergent perturbation expansion can be obtained. The transformed Hamiltonian is

$`\stackrel{~}{H_d}`$ $`=`$ $`e^RH_de^R`$ (3)
$`=`$ $`\omega _0d^{}d+{\displaystyle \underset{i}{}}ϵ_pn_it\mathrm{cos}({\displaystyle \frac{\theta }{2}})[c_1^{}c_2\mathrm{exp}(2\lambda (d^{}d))`$ (4)
$`+`$ $`c_2^{}c_1\mathrm{exp}(2\lambda (d^{}d))]+\omega _0(g_+\lambda )(n_1n_2)(d+d^{})+JS^2\mathrm{cos}\theta `$ (5)
where $`R=\lambda (n_1n_2)(d^{}d)`$, $`\lambda `$ is a variational parameter related to the displacement of the $`d`$ oscillator, $`g_+=g/\sqrt{2}`$ and $`ϵ_p=ϵ\omega _0(2g_+\lambda )\lambda `$.

For the perturbative expansion following Ref.  the basis set is chosen as $`|\pm ,N=\frac{1}{\sqrt{2}}(c_1^{}\pm c_2^{})`$ $`|0_e|N`$, where $`|+`$ and $`|`$ are the bonding and antibonding electronic states and $`|N`$ denotes the $`N`$th excited oscillator state within the MLF phonon basis. The diagonal part of the Hamiltonian $`\stackrel{~}{H_d}`$ in the chosen basis is treated as the unperturbed Hamiltonian ($`H_0`$) and the remaining part of the Hamiltonian $`H_1=\stackrel{~}{H_d}H_0`$, as a perturbation.

The unperturbed energy of the state $`|\pm ,N`$ is given by

$`E_{\pm ,N}^{(0)}`$ $`=`$ $`N,\pm |H_0|\pm ,N`$ (6)
$`=`$ $`N\omega _0+ϵ_pt_{eff}\left[{\displaystyle \underset{i=0}{\overset{N}{}}}{\displaystyle \frac{(2\lambda )^{2i}}{i!}}(1)^iN_{C_i}\right]+JS^2\mathrm{cos}\theta `$ (7)
where $`t_{eff}=t\mathrm{cos}\frac{\theta }{2}\mathrm{exp}(2\lambda ^2)`$ and the general off-diagonal matrix elements of $`H_1`$ between the two states $`|\pm ,N`$ and $`|\pm ,M`$ may be calculated for $`(NM)>0`$ as in Ref. .

Within the chosen basis, the unperturbed ground state is the $`|+|0`$ state with the unperturbed energy, $`E_0^{(0)}=ϵ_pt_{eff}+JS^2\mathrm{cos}\theta `$.

The first order correction to the ground state wave function is obtained as,

$`|\psi _0^{(1)}={\displaystyle \frac{[\omega _0(g_+\lambda )2\lambda t_{eff}]}{(E_0^{(0)}E_{+,1}^{(0)})}}|,1{\displaystyle \underset{N=2,3,4,..}{}}{\displaystyle \frac{t_{eff}(2\lambda )^N}{\sqrt{N!}(E_0^{(0)}E_{e,N}^{(0)})}}|e,N`$ (8)
where e = + or - for even and odd N respectively.

The first order correction to the energy ($`E_0^{(1)}`$) is zero since $`H_1`$ has no diagonal matrix element in the chosen basis. The second order correction to the ground state energy is given by

$`E_0^{(2)}`$ $`=`$ $`{\displaystyle \underset{N=1,3..}{}}{\displaystyle \frac{|\frac{t_{eff}(2\lambda )^N}{\sqrt{N!}}+\omega _0(g_+\lambda )\delta _{N,1}|^2}{(E_0^{(0)}E_{,N}^{(0)})}}`$ (9)
$`+`$ $`{\displaystyle \underset{N=2,4..}{}}{\displaystyle \frac{|\frac{t_{eff}(2\lambda )^N}{\sqrt{N!}}|^2}{(E_0^{(0)}E_{+,N}^{(0)})}}`$ (10)
Higher order corrections to the ground state wavefunction and energy are obtained following Ref. . For the study of the effect of the magnetic field ($`h`$) we include a term -$`\stackrel{~}{g}_i\mu _BhS_i\mathrm{cos}\frac{\theta }{2}`$ to the Hamiltonian in Eq. (1) (where $`\stackrel{~}{g}`$ is the Lande g factor) and as a result a term $`2\mu _{eff}h\mathrm{cos}\frac{\theta }{2}`$ is added to $`E_0^{(0)}`$ where $`\mu _{eff}`$ (= $`\stackrel{~}{g}S\mu _B`$) is the local moment of the core spins. In this paper we express the magnetic field ($`h`$) in a unit of $`\mu _{eff}=1`$.

Now a proper choice of $`\lambda `$ is to be made so that the perturbative expansion becomes convergent. Our previous work  has shown that the $`\lambda `$, obtained by minimizing the unperturbed ground state energy, gives satisfactory convergence to the perturbation series for $`t/\omega _01`$ and the convergence becomes very rapid and excellent with decreasing value of $`t/\omega _0`$. Here also we will follow the similar procedure. Minimizing the unperturbed ground state energy $`E_0^{(0)}`$ with respect to $`\lambda `$ we obtain

$`\lambda `$ $`=`$ $`{\displaystyle \frac{\omega _0g_+}{\omega _0+2t_{eff}}}`$ (11)
Minimization of the unperturbed ground state energy with respect to $`\theta `$ gives an approximate value ($`\theta _{\mathrm{MLF}}`$) of $`\theta `$

$`\mathrm{cos}{\displaystyle \frac{\theta _{\mathrm{MLF}}}{2}}`$ $`=`$ $`[{\displaystyle \frac{t\mathrm{exp}(2\lambda ^2)}{4JS^2}}+{\displaystyle \frac{\mu _{eff}h}{2JS^2}}](\mathrm{for}\mathrm{nonzero}\mathrm{solution}\mathrm{of}\theta _{\mathrm{MLF}})`$ (12)
However, exact $`\theta `$ should be evaluated from the minimization of the exact ground state energy. For each value of $`g_+`$ we calculate the energy upto the sixth order in perturbation  and find out for which value of $`\theta `$ the energy (including perturbation corrections) is minimum.

We have also calculated $`t_{eff}^{KE}=E_{Kin}=<\psi _G|`$ $`t\mathrm{cos}(\frac{\theta }{2})[c_1^{}c_2\mathrm{exp}`$ $`(2\lambda (d^{}d))`$ $`+c_2^{}c_1\mathrm{exp}(2\lambda (d^{}d))]`$ $`|\psi _G>`$ where $`\psi _G`$ is the ground state wave-funcion which we have calculated upto the fifth order corrections in perturbation. $`t_{eff}^{KE}`$ describes the kinetic energy of the system and it reduces to the effective hopping ($`t_{eff}`$) for both the weak and strong coupling limits .

The static correlation functions $`n_1u_1_0`$ and $`n_1u_2_0`$, where $`u_1`$ and $`u_2`$ are the lattice deformations at sites 1 and 2 respectively, produced by an electron at site 1, are the standard measure of polaronic character. The correlation functions for this two-site system may be written as

$`n_1u_1_0`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[(g_++\lambda )+{\displaystyle \frac{A_0}{N_G}}\right]`$ (13)
$`n_1u_2_0`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[(g_+\lambda ){\displaystyle \frac{A_0}{N_G}}\right]`$ (14)
$`\mathrm{where}A_0\psi _G|n_1(d+d^{})|\psi _G`$ (15)
where $`N_G`$ is the normalization factor to the wavefunction. The variation of the physical quantity $`\lambda _{corr}/g_+=<n_1(u_1u_2)>_0/g_+`$ with the $`e`$-ph coupling strength manifests the nature of large to small polaron crossover. Its value becomes 1 in the extreme small polaronic limit where the lattice distortion is very local and decreases significantly in the small to large polaron crossover region as $`g_+`$ decreases.

3. Results and Discussions

In the combined model of polaron and double exchange interaction, the nature of variation of $`\theta `$, $`t_{eff}^{KE}`$ and $`\lambda _{corr}/g_+`$ with $`g_+`$ reveal different transitions and crossover regions. We focus our attention to the parameter space where $`t>J`$ as it corresponds to real systems . We present the results for $`t=1.0`$ and $`JS^2=0.05,0.15`$, 0.25 and $`0.5`$ (in a scale of $`\omega _0=1.0`$). The bare-site energy $`ϵ`$ is taken as zero.

Fig. 1 shows the change of relative orientation ($`\theta `$) of two core spins as a function of $`e`$-ph coupling strength ($`g_+`$) for the ground state. For small values of $`g_+`$ two core spins are either aligned parallel or canted depending on the strength of $`t/JS^2`$. With increasing $`g_+`$ the angle ($`\theta `$) increases and finally a transition from the FM or canted AFM state to the AFM state occurs. Nature of this transition depends on the ratio of $`t`$ to $`JS^2`$. For $`t`$=1.0, $`JS^2`$=0.05 the FM to AFM transition is very sharp. With increasing value of $`JS^2`$ the transition occurs at a lower value of $`g_+`$ and it becomes broader. For higher values of $`JS^2`$ (=0.25, 0.5), the crossover from the canted state to the AFM state is very smooth. It may be mentioned that Yunoki and Moreo  studied the 1-orbital Kondo model with AF interaction ($`J`$) between the core spins (in absence of $`e`$-ph interaction) in the context of manganites and obtained FM state at $`n=0.5`$ for low values of $`J`$. With increasing $`J`$ the ground state at quarter filling ($`n=0.5`$) appears to be a spin spiral state and ultimately an AF state. For the two orbital model in presence of $`e`$-ph coupling a similar FM state is observed for low values of $`e`$-ph coupling for $`n=1`$ and an AF state for large values of $`e`$-ph coupling.

With increasing $`g_+`$ a large to small polaron crossover occurs in the system in addition to the magnetic transition. This is clearly seen from Figs. 2 and 3 where the variations of $`t_{eff}^{KE}`$, $`\lambda _{corr}/g_+`$ and $`\theta `$ with $`g_+`$ are shown. The large to small polaron crossover is identified by the sharp fall in $`t_{eff}^{KE}`$ and rise in $`\lambda _{corr}`$. For $`t=1.0`$ and $`JS^2=0.15`$ the FM to AFM transition and the large to small polaron crossover occur simultaneously (Fig. 2) around $`g_+=1.2`$. However, for $`JS^2`$=0.05 the magnetic transition occurs at a higher values of $`g_+`$ than that coresponding to large to small polaron crossover (Fig. 3). In this situation the FM state has two regions: one with appreciable value of $`t_{eff}^{KE}`$ and the other with much reduced value of $`t_{eff}^{KE}`$ (Fig. 3). In the former case, charge carrier is the delocalized large polaron with appreciable hopping and this regime is expected to be metallic in the thermodynamic limit. In the latter region the hopping would be heavily suppressed by the small polaron formation and it may result in an insulating FM state in some cases. Thus we find that there is a possibility of FM metal-FM insulator-AFM insulator transition for low values of $`JS^2/t`$. It is to be noted that if one follows the usual MLF method (zeroth order of perturbation within our approach) the FM-AFM transition coincides with large to small polaron crossover for any value of $`JS^2/t`$. This is demonstrated in Fig. 3 for $`JS^2/t`$=0.05, where we find the coincidence of the FM-AFM transion with the large-small polaron crossover within the usual MLF approach whereas the nearly exact calculations yield different result. The nearly exact results, obtained from the MLF perturbation method, makes the large-small polaron crossover smooth in favour of the conclusion of Löwen .

In Fig. 4 we have shown the effect of the external magnetic field on $`\theta `$, $`\lambda _{corr}/g_+`$ and $`t_{eff}^{KE}`$ for $`JS^2=0.15`$. The magnetic field favours the FM state, hence the FM-AFM transition and associated large-small polaron crossover (for $`JS^2/t=0.15`$) take place at a higher value of $`g_+`$ with increasing field. From Fig. 4 it is clear that the magnetic field reduces the value of the local distortion ($`\lambda _{corr}/g_+`$) and enhances $`t_{eff}^{KE}`$ significantly in the large-small polaron crossover region provided the crossover is associated with a FM-AFM transition.

For $`JS^2/t=0.05`$ the system remains in the FM state in the large-small polaron crossover region, consequently the magnetic field has no effect on $`t_{eff}^{KE}`$ in this crossover region. However, the field has an effect on $`\theta `$, hence on $`t_{eff}^{KE}`$ at the FM-AFM transition.

In Fig. 5 we plot the change in $`t_{eff}^{KE}`$ due to the magnetic field as a function of $`g_+`$. This quantity may be related to the magnetoresistance for a system in the thermodynamic limit. For polarons the conductiviy at low temperatures is dominated by the tunneling mobility which is proportional to $`t_{eff}^2`$ within the zeroth order of perturbation . In general the $`t_{eff}^{KE}`$ is a measure of delocalization of the electron. A reduction of $`t_{eff}^{KE}`$ will cause a reduction in mobility. In Fig. 5 the change in $`t_{eff}^{KE}`$ due to the field is very sharp and prominent at the FM-AFM transtion when it is associated with the large-small polaron crossover for $`JS^2/t=0.15`$, while it shows a broad smaller peak for $`JS^2/t=0.05`$ where the FM-AFM transition is not associated with the polaron crossover. For the former case the magnetic field would have pronounced effect on the transport properties.

4. Conclusions

From our studies on the two-site double exchange model with a single polaron as a function of $`e`$-ph coupling strength we conclude that the nature of the FM-AFM transition depends on the relative values of $`J`$ and $`t`$. The transition is sharper for smaller values of $`J/t`$. For high values of $`J/t`$ a canted state is stable instead of a FM state for weak $`e`$-ph coupling ($`g_+`$) and the crossover from the canted state to the AFM state is very smooth with increasing $`g_+`$. For suitable values of $`J/t`$ (=0.15) the FM-AFM transition coincides with the large-small polaron crossover. For this case the external magnetic field has very prominent effect on the polaronic local distortion and kinetic energy in the transition region (as in Figs. 4(b) and 4(c) respectively). Within the MLF method (zeroth order of perturbation) the polaron crossover is always associated with FM-AFM transition for any value of $`J/t`$ whereas our results based on convergent perturbation expansion yield that the coincidence of magnetic transition and polaron crossover will depend on the value of $`t`$ and $`J`$. For low values of $`J/t`$ a crossover from the FM large polaronic state to the FM small polaronic state with reduced hopping occurs and then a transition to the AFM state takes place at a higher value of $`g_+`$.

Figure captions :

FIG. 1. Variation of the relative angle ($`\theta `$) between the core-spins, with $`g_+`$ for $`t=1.0`$ and $`JS^2=0.05`$, 0.15, 0.25 and 0.5 (in units of $`\omega _0`$=1).

FIG. 2. Variations of $`\theta `$, $`t_{eff}^{KE}`$ and $`\lambda _{corr}/g_+`$ with $`g_+`$ for $`t=1.0`$, $`JS^2=0.15`$ and $`h=0`$.

FIG. 3. Variations of $`\theta `$ (with and without perturbation), $`t_{eff}^{KE}`$, $`t_{eff}`$, $`\lambda _{corr}/g_+`$ and $`\lambda /g_+`$ with $`g_+`$ for $`t=1.0`$, $`JS^2=0.05`$ and $`h=0`$. Note that $`\theta `$(MLF), $`t_{eff}`$ and $`\lambda /g_+`$ are obtained within the MLF method without considering perturbation corrections.

FIG. 4. Variation of (a) $`\theta `$, (b) $`\lambda _{corr}/g_+`$ and (c) $`t_{eff}^{KE}`$ with $`g_+`$ for different values of the magnetic field $`h=`$0, 0.02 and 0.04 (in units of $`\mu _{eff}`$=1) for $`t=1.0`$ and $`JS^2=0.15`$.

FIG. 5. Variation of the change in $`t_{eff}^{KE}`$ due to the magnetic field ($`\frac{(t_{eff}^{KE}(h)t_{eff}^{KE}(0))}{h}`$) as a function of $`g_+`$ for $`t=1.0`$, $`JS^2=0.05`$ and 0.15 and $`h`$=0.04.


text:
# Quantum phonons and the charge density wave transition temperature: a dynamical mean field study

## I Introduction

The present-day understanding of electron-lattice interactions in metals is based on the Migdal-Eliashberg (ME) theory, which exploits the smallness of typical phonon frequencies, $`\mathrm{\Omega }`$, and temperatures, $`T`$, relative to typical electronic energies, $`t`$. This theory may be regarded as an expansion in powers of the adiabatic parameter $`\gamma =\mathrm{max}(\mathrm{\Omega },T)/t`$, in which only the leading nontrivial term is retained. In many cases the leading term is of order $`\gamma ^0`$ and thus within ME theory most quantities of physical interest are independent of ion mass. A notable exception is the superconducting transition temperature $`T_c`$, which is controlled by a logarithmic divergence with an upper cutoff set by the phonon frequency, and involves an interaction (the Coulomb pseudopotential) which has a logarithmic suppression with a lower cutoff set by the phonon frequency. The resulting, rather complicated, dependence of $`T_c`$ on ion mass was determined by McMillan.

A series of recent experimental papers reporting isotope effects (i.e. ion-mass dependence) for a variety of electronic properties  motivates re-examination of this issue. Of particular interest here are reports of a pronounced increase of the charge ordering transition temperatures of various members of the ‘CMR manganite’ family of manganese perovskites when $`{}_{}{}^{18}\mathrm{O}`$ is substituted for $`{}_{}{}^{16}\mathrm{O}`$. This leads us to re-examine the theory of charge density wave (CDW) instabilities in order to determine how the transition temperature $`T_{\mathrm{co}}`$ is affected by ion mass. A theory applicable to materials of current interest must go beyond the Migdal-Eliashberg approximation, must apply to the case of $`d=3`$ spatial dimensions and must be able to treat the physically relevant limit $`\gamma 1`$.

The $`\mathrm{\Omega }`$-dependence of $`T_{\mathrm{co}}`$ has been studied. In one spatial dimension pioneering numerical work of Hirsch and Fradkin followed by analytical studies by Bourbonnais and Caron found a strong suppression of the zero temperature CDW order parameter with increasing $`\mathrm{\Omega }`$. In $`d=2`$  and $`d=\mathrm{}`$ the quantum phonon problem has been attacked by a direct numerical approach using quantum Monte Carlo (QMC) methods. However, numerical issues related to the mismatch between phonon and electron frequency scales however restricted the QMC studies to the antiadiabatic limit, $`\mathrm{\Omega }t`$. Also, in Refs.  the isotope effect on the CDW transition temperature was not studied in detail, although values for $`T_{\mathrm{co}}`$ for a given phonon frequency were calculated. In $`d=\mathrm{}`$ and for weak coupling strengths the CDW transition temperature is strongly reduced with decreasing ion mass and a reduction ratio of $`T_{\mathrm{co}}(\mathrm{\Omega }=0.5t)/T_{\mathrm{co}}(\mathrm{\Omega }=0)0.25`$ has been reported.

Especially in $`d=\mathrm{}`$ considerable effort has been expended to construct weak-coupling perturbation theories which continuously connect the adiabatic to the antiadiabatic limit. One important outcome of this work was the suggestion that future research should focus on generalizations of the ME theory that work with dressed phonons (renormalized ME) but include higher order non-adiabatic effects such as vertex corrections. Our method is such an approach. A large number of studies based on the Lang-Firsov method have been published (see e.g. Ref. ). For one dimensional systems the Lang-Firsov results are in agreement with previous work. However, this is an uncontrolled approximation optimized for the single polaron problem and its applicability to metallic densities is not established.

In this paper we construct an essentially analytic theory of the charge density wave transition temperature in the physically relevant limit $`\gamma 1`$. We build on previous work which showed how to construct an expansion in $`\gamma `$ within the dynamical mean-field approximation and which analysed the classical limit ($`\mathrm{\Omega }=0`$) of the charge density wave problem. Here we show how to calculate $`T_{\mathrm{co}}`$ including the first non-Migdal term in $`\gamma `$ for arbitrary values of $`\mathrm{\Omega }/T_{\mathrm{co}}`$. We find that $`T_{\mathrm{co}}`$ depends very strongly on $`\mathrm{\Omega }`$, becoming very small as $`\mathrm{\Omega }`$ is increased (at fixed coupling) beyond the $`\mathrm{\Omega }=0`$ value of $`T_{\mathrm{co}}`$, but that the specific behavior depends strongly on whether spin-1/2 or spinless electrons are considered. For spinless electrons, relevant to half-metals such as CMR manganites, a divergent isotope exponent is obtained in the weak electron-phonon coupling limit.

## II Model and Formalism

In this paper we study the simplest model of electrons interacting with phonons, namely the Holstein model $`H_{\mathrm{Hol}}=H_{\mathrm{el}}+H_{\mathrm{ph}}+H_{\mathrm{el}\mathrm{ph}}`$ with

$`H_{\mathrm{el}}`$ $`=`$ $`{\displaystyle \underset{ij}{}}t_{ij}c_{i\sigma }^{}c_{j\sigma }^{}\mu {\displaystyle \underset{i}{}}\left(c_{i\sigma }^{}c_{i\sigma }^{}n\right),`$ (1)
$`H_{\mathrm{ph}}`$ $`=`$ $`{\displaystyle \frac{1}{2\mathrm{\Lambda }}}{\displaystyle \underset{i}{}}\left(r_i^{\mathrm{\hspace{0.17em}2}}+\dot{r}_i^2/\mathrm{\Omega }^2\right),`$ (2)
$`H_{\mathrm{el}\mathrm{ph}}`$ $`=`$ $`{\displaystyle \underset{i}{}}r_i\left(c_{i\sigma }^{}c_{i\sigma }^{}n\right)`$ (3)
where we have absorbed the electron-phonon coupling into the phonon coordinate $`r`$ which thus has dimension of energy as does the phonon stiffness parameter $`\mathrm{\Lambda }`$. We further specialize to a mean electron density per spin direction of $`n=1/2`$; this implies $`\mu =0`$. We also assume a bipartite lattice, which for our purposes we define as a lattice possessing a dispersion $`ϵ_\stackrel{}{k}`$ (Fourier transform of $`t_{ij}`$) and a wavevector $`\stackrel{}{Q}`$ such that $`ϵ_{\stackrel{}{k}\stackrel{}{Q}}=ϵ_\stackrel{}{k}`$ for all $`\stackrel{}{k}`$ in the Brillouin zone. An equivalent definition is that the lattice may be divided into two sublattices $`A`$ and $`B`$, such that $`t_{ij}`$ connects $`A`$ sites only to $`B`$ sites. The important feature of the dispersion is the bare density of states per site per spin $`\rho (ϵ)=\frac{1}{N}_\stackrel{}{k}\delta (ϵϵ_\stackrel{}{k})`$. In this paper we use the semicircular form $`\rho (ϵ)=\frac{1}{2\pi t^2}\sqrt{4t^2ϵ^2}`$. It is useful to introduce the dimensionless electron-phonon coupling $`\lambda `$ and adiabatic parameter $`\gamma `$ by

$$\lambda =\mathrm{\Lambda }\rho _0,\gamma =\mathrm{\Omega }\rho _0.$$
(4)
Here $`\rho _0=1/\pi t`$ denotes the bare density of states at the chemical potential.

We now establish notation by reviewing the derivation of the $`T_{\mathrm{co}}`$ equation within dynamical mean-field theory (DMFT). The fundamental assumption of DMFT is the momentum independence of the electron self energy $`\mathrm{\Sigma }_n=\mathrm{\Sigma }(i\omega _n)`$ ($`\omega _n=2\pi T(n+1/2)`$). Consequently, in the absence of spatial symmetry breaking, $`\mathrm{\Sigma }_n`$ can be derived from an effective impurity model specified by the action

$`S(\{r_{nm}\};\{c_n\})`$ $`=`$ $`(T/2\mathrm{\Lambda }){\displaystyle \underset{n}{}}r_{n0}\left(\omega _{n0}^2/\mathrm{\Omega }^2+1\right)r_{0n}`$ (6)
$`n{\displaystyle \underset{n}{}}r_{n0}(2s+1)\mathrm{Tr}\mathrm{ln}\left[c_n\delta _{nm}Tr_{nm}\right].`$
In this action $`r_{nm}`$ are bosonic fields describing the phonons and $`c_n`$ mean field functions describing the fermions (which have been integrated out). We have indexed the bosonic fields by the difference of two fermionic Matsubara frequencies $`\omega _{nm}=\omega _n\omega _m`$. The factor $`(2s+1)`$ takes the values one and two for spinless ($`s=0`$) and spin-1/2 ($`s=1/2`$) electrons, respectively. The partition function is the functional integral over the bosonic fields

$$Z=𝒟[r]\mathrm{exp}(S)$$
(7)
and is a functional of the effective field $`c`$ alone. It is useful to define the impurity Green’s function $`𝒢`$ and self energy $`\mathrm{\Sigma }`$ by

$$𝒢_n\frac{1}{2s+1}\frac{\delta \mathrm{ln}Z}{\delta c_n}\frac{1}{2s+1}\frac{1}{c_n\mathrm{\Sigma }_n(\{c_n\})}.$$
(8)
The effective field is fixed by equating the local Green’s functions of the original lattice model and the effective impurity model. For a semicircular density of states one obtains

$$c_n=i\omega _n+\mu t^2𝒢_n\left(\{c_n\}\right).$$
(9)
For temperatures below the CDW transition the self energy acquires two components for the two sublattices $`A`$ and $`B`$ which calls for two effective impurity models given by the effective fields $`a`$ and $`b`$ respectively. The resulting self-consistency equations are:

$`a_n`$ $`=`$ $`i\omega _n+\mu t^2𝒢_n\left(\{b_n\}\right)`$ (10)
$`b_n`$ $`=`$ $`i\omega _n+\mu t^2𝒢_n\left(\{a_n\}\right).`$ (11)
Near the transition (assumed to be second order) we can expand $`a_n=c_n+ϵ_n^{}`$ and $`b_n=c_nϵ_n^{}`$ leading to $`ϵ_n^{}=t^2_m\frac{𝒢_n}{c_m}ϵ_m^{}`$. Making use of Eq. (8) we obtain

$$ϵ_n^{}=\frac{1}{1+(c_n\mathrm{\Sigma }_n)^2/t^2}\underset{m}{}\frac{\mathrm{\Sigma }_n}{c_m}ϵ_m^{}.$$
(12)
It is convenient to rewrite Eq. (12) in the more symmetric form

$$ϵ_n=\underset{m}{}\chi _n^{1/2}\mathrm{\Gamma }_{nm}\chi _m^{1/2}ϵ_m.$$
(13)
Here we have introduced $`\chi _n=𝒢_n/(i\omega _n\mathrm{\Sigma }_n)`$, the irreducible two particle vertex $`\mathrm{\Gamma }_{nm}=_\sigma ^{}\mathrm{\Sigma }_n^\sigma /𝒢_m^\sigma ^{}`$ and $`ϵ_n=ϵ_n^{}\chi _n^{1/2}`$. Eq. (13) is the statement that the matrix $`\chi _n^{1/2}\mathrm{\Gamma }_{nm}\chi _m^{1/2}`$ has eigenvalue one and $`T_{\mathrm{co}}`$ is the temperature at which this eigenvalue problem has a solution. It can be easily shown that $`\chi _n`$ is the ‘bare’ CDW susceptibility of the original lattice model $`\chi _n=1/N_\stackrel{}{k}G_n(\stackrel{}{k}+\stackrel{}{Q})G_n(\stackrel{}{k})`$.

## III Solution in the adiabatic limit

As written, Eq. (13) is general, but not very useful because not only must the matrix be diagonalized numerically, but also determination of the matrix elements for arbitrary $`\lambda `$ and $`\gamma `$ requires performing the functional integral over phonon fields numerically (i.e. via quantum Monte Carlo) and this becomes prohibitively lengthly for physically relevant regimes of $`\mathrm{\Omega }`$ and $`T`$. Here we turn this difficulty to advantage, exploiting the smallness of $`\mathrm{\Omega }`$ and $`T`$ to construct an analytic expansion in powers of $`\gamma `$. The expansion builds on the observation that the biggest contribution to a diagram (e.g. for the electronic self energy or the vertex) with $`n`$ phonon loops is of order $`[\mathrm{max}(T,\mathrm{\Omega })/t]^n`$. This can be seen as follows. A generic diagram $`D^{(n)}`$ with $`n`$ phonon loops contains $`n`$ internal sums over bosonic Matsubara frequencies $`\omega _1,\mathrm{},\omega _n`$ (here $`\omega _n=\omega _{n0}`$) each associated with one (or more) free phonon propagators given by

$$𝒟_n^0=\mathrm{\Lambda }\frac{\mathrm{\Omega }^2}{\mathrm{\Omega }^2+\omega _n^2}.$$
(14)
If $`F(i\omega _1,\mathrm{},i\omega _n)`$ depicts the appropriate product of Green’s functions and phonon propagators we can write $`D^{(n)}=T^n_{1,\mathrm{},n}F(i\omega _1,\mathrm{},i\omega _n)`$. In the classical limit, $`\mathrm{\Omega }/T0`$, the free phonon propagators become delta functions and therefore $`D^{(n)}T^n`$ as noted in Ref. . In the quantum case it is more appropriate to rescale the temperature $`\stackrel{~}{T}=T/\gamma `$ leading to $`D^{(n)}=\gamma ^n\stackrel{~}{T}^n_{1,\mathrm{},n}F(i\gamma \stackrel{~}{\omega }_1,\mathrm{},i\gamma \stackrel{~}{\omega }_n)`$. The rescaled bosonic Matsubara frequencies $`\stackrel{~}{\omega }`$ become continuous variables in the extreme quantum limit $`\stackrel{~}{T}0`$ and because $`F`$ has a small frequency expansion we conclude that the biggest contributions to the $`n`$ phonon loop diagram are of order $`\gamma ^n`$ as argued in Ref. . The diagrammatic expansion shown in Fig. 1 includes all terms with the minimal number of internal phonon loops, i.e. includes the leading nontrivial order for the electronic self energy and the vertex. Note that the leading order in $`\gamma `$ contains both, the direct (first term in Fig. 1 (b)) and the exchange diagrams for the vertex. The conventionally defined ME approach only sums the bubble diagrams for the CDW susceptibility, i.e. the exchange diagrams for the vertex (which we find to be of the same order) are neglected. Freericks and co-workers have considered a variety of different vertex corrections to the conventionally defined ME approximation but they organized these in a coupling constant expansion and not in the small frequency expansion we considered here. Specific vertex corrections are of higher order in $`\gamma `$ than the ones treated in this work.

Adding the electron bubbles in Fig. 1 leads to a ‘dressed’ phonon propagator $`𝒟_n=[(𝒟_n^0)^1\mathrm{\Pi }_n]^1`$ with a phonon self energy $`\mathrm{\Pi }_n`$ given by

$$\mathrm{\Pi }_n=(2s+1)T\underset{m}{}𝒢_{n+m}𝒢_m.$$
(15)
The occurrence of the free phonon propagator restricts the bosonic Matsubara frequencies $`\omega _n`$ to values on the scale of $`\mathrm{\Omega }`$. Following the arguments given in Ref.  we assume that the dominant contributions to $`\mathrm{\Pi }_n`$ come from frequencies of order $`t`$ and consequently set the index $`n`$ under the sum to zero. To leading order in $`\gamma `$ we can replace the full by the free Green’s function and find $`\mathrm{\Pi }_n=1/\mathrm{\Lambda }_c+𝒪(T^2)`$ with $`\mathrm{\Lambda }_c=3/[4(2s+1)\rho _0]`$. The ‘dressed’ phonon propagator reads

$$𝒟_n=\overline{\mathrm{\Lambda }}\frac{\overline{\mathrm{\Omega }}^2}{\overline{\mathrm{\Omega }}^2+\omega _n^2}$$
(16)
with the renormalized phonon stiffness and phonon frequency

$`\overline{\mathrm{\Lambda }}`$ $`=`$ $`\mathrm{\Lambda }/(1\mathrm{\Lambda }/\mathrm{\Lambda }_c)`$ (17)
$`\overline{\mathrm{\Omega }}`$ $`=`$ $`\mathrm{\Omega }(1\mathrm{\Lambda }/\mathrm{\Lambda }_c)^{1/2}.`$ (18)
Due to this influence of the electrons on the phonons the given expansion is actually an expansion in $`\overline{\lambda }=\overline{\mathrm{\Lambda }}\rho _0`$ and $`\overline{\gamma }=\overline{\mathrm{\Omega }}\rho _0`$. Maybe a more familiar way to express the dressed electron phonon coupling is $`\overline{\lambda }=\rho _0\mathrm{\Lambda }𝒟(0)`$. The importance of not neglecting the renormalization of phonons by electrons was stressed by several authors, both in studies of $`d=2`$ and $`d=\mathrm{}`$ spatial dimensions. In the limit $`\mathrm{\Lambda }\mathrm{\Lambda }_c`$ the renormalized phonon stiffness $`\overline{\mathrm{\Lambda }}\mathrm{}`$ and the renormalized phonon frequency $`\overline{\mathrm{\Omega }}0`$ signaling that the lattice becomes unstable against local distortions. In the absence of long range order the ground state expectation value of the lattice distortion $`r`$ would change from $`r=0`$ for $`\mathrm{\Lambda }<\mathrm{\Lambda }_c`$ to $`r0`$ for $`\mathrm{\Lambda }>\mathrm{\Lambda }_c`$. This polaronic instability has been extensively discussed in Ref.  and it is important to note that the coupling constant is strongly enhanced and the adiabatic parameter is strongly suppressed near the instability.

Within the small $`\gamma `$ limit we obtain

$`\mathrm{\Gamma }_{nm}`$ $`=`$ $`{\displaystyle \frac{T}{\rho _0}}\left[(2s+1)\lambda +\overline{\lambda }{\displaystyle \frac{1}{1+\stackrel{~}{\omega }_{nm}^2}}\right]`$ (20)
$`+𝒪\left(\mathrm{max}[T^2/t^2,\overline{\gamma }^2]\right)`$
and

$`\mathrm{\Sigma }_n`$ $`=`$ $`{\displaystyle \frac{\overline{\lambda }T}{\rho _0}}{\displaystyle \underset{m}{}}𝒢_m{\displaystyle \frac{1}{1+\stackrel{~}{\omega }_{nm}^2}}`$ (22)
$`+𝒪\left(\mathrm{max}[T^2/t^2,\overline{\gamma }^2]\right).`$
Here, we have rescaled the bosonic frequencies $`\stackrel{~}{\omega }_{nm}=\omega _{nm}/\overline{\mathrm{\Omega }}`$. In Eq. (22) as well as in $`\chi _n`$, we can replace the full Green’s function by the bare one $`𝒢_n𝒢_n^0=[i\omega _ni\mathrm{sgn}(\omega _n)\sqrt{\omega _n^2+4t^2}]/2t^2`$ to leading order in $`\overline{\gamma }`$. Previous studies asserted that for small phonon frequencies vertex corrections such as the second term in Eq. (20) can be neglected. We see that this is not quite correct: rather, the vertex corrections for $`\overline{\gamma }0`$ affect only the diagonal elements of the matrix in Eq. (13), but in this diagonal region they are of the same order as the contribution from the first part of Eq. (20).

The importance of the vertex corrections has been demonstrated in the classical limit. For $`\overline{\gamma }0`$ ($`\overline{\mathrm{\Omega }}0`$) we can expand

$$\frac{1}{1+\stackrel{~}{\omega }_{nm}^2}\delta _{nm}+\frac{1}{\stackrel{~}{\omega }_{nm}^2}\left(1\delta _{nm}\right)\delta _{nm},$$
(23)
i.e. the phonon propagator becomes diagonal for $`\overline{\gamma }=0`$. In this limit both vertex and self energy corrections to the BCS case modify the upper cutoff of the logarithmic divergence in the CDW susceptibility leading to a suppression of the CDW transition temperature, which has been interpreted in terms of inelastic scattering of electrons by phonons.

Returning now to the quantum problem we insert Eq. (20) into the eigenvalue equation (13) yielding

$`ϵ_n`$ $`=`$ $`\chi _n^{1/2}(2s+1){\displaystyle \frac{\lambda T}{\rho _0}}{\displaystyle \underset{m}{}}\chi _m^{1/2}ϵ_m`$ (25)
$`\chi _n^{1/2}{\displaystyle \frac{\overline{\lambda }T}{\rho _0}}{\displaystyle \underset{m}{}}{\displaystyle \frac{1}{1+\stackrel{~}{\omega }_{nm}^2}}\chi _m^{1/2}ϵ_m.`$
For generic values of $`\mathrm{\Omega }/T_{\mathrm{co}}(0)`$ we solve Eq. (25) numerically by restricting the sums to a finite range $`m<N`$ and diagonalizing the resulting matrix exactly numerically. Convergence is improved by using an analytical approximation for Eq. (25) for large frequencies $`\omega _n>\omega _N`$. In the large frequency limit the sums in Eq. (22) and Eq. (25) including the phonon propagator equal their $`n=m`$ values times the factor $`T_m1/(1+\stackrel{~}{\omega }_m^2)=(\overline{\mathrm{\Omega }}/2)\mathrm{coth}(\overline{\mathrm{\Omega }}/2T)`$. Therefore, the coupling constant $`(2s+1)\lambda `$ can be replaced by the effective coupling constant

$$\lambda ^{}=(2s+1)\lambda /[1(2s+1)\lambda I)]$$
(26)
with

$$I=\frac{T}{\rho _0}\underset{|n|>N}{}\frac{1}{(\chi _n^0)^1+(\overline{\lambda }\overline{\mathrm{\Omega }}/\rho _0)\mathrm{coth}\left(\frac{\overline{\mathrm{\Omega }}}{2T}\right)}$$
(27)
where we have introduced $`\chi _n^0=𝒢_n^0/i\omega _n`$. For $`n,m=0,\mathrm{},N`$ the leading eigenvalue of the matrix

$`M_{nm}`$ $`=`$ $`{\displaystyle \frac{2\lambda ^{}T}{\rho _0}}\chi _n^{1/2}\chi _m^{1/2}{\displaystyle \frac{\overline{\lambda }T}{\rho _0}}\chi _n^{1/2}\chi _m^{1/2}`$ (29)
$`\times \left\{{\displaystyle \frac{1}{1+[2\pi (nm)T/\overline{\mathrm{\Omega }}]^2}}+{\displaystyle \frac{1}{1+[2\pi (n+m)T/\overline{\mathrm{\Omega }}]^2}}\right\}`$
has to be found. $`T_{\mathrm{co}}`$ is given by the condition that the leading eigenvalue equals one. We have checked the convergence with $`N`$; typically we used values ranging from $`N=256`$ to $`N=4096`$ yielding a relative error $`\delta T_{\mathrm{co}}/T_{\mathrm{co}}<0.01`$. For later use we give the weak coupling $`\overline{\lambda }0`$ result for the effective coupling constant. In this limit the charge density wave transition temperature becomes very small and the Matsubara sum can be evaluated as an integral yielding $`1/\lambda ^{}1/[(2s+1)\lambda ]\mathrm{ln}(2t/\omega _N)`$.

## IV Results

The CDW transition temperature resulting from a numerical diagonalization of (29) is shown in Fig. 2. Decreasing the ionic mass, i.e. increasing the phonon frequency $`\mathrm{\Omega }`$, is seen to suppress $`T_{\mathrm{co}}`$ considerably for both spinless and spin-1/2 electrons. The suppression becomes very pronounced when $`\mathrm{\Omega }(2s+1)\pi T_{\mathrm{co}}(\mathrm{\Omega }=0)`$. This is in qualitative agreement with the strong suppression of the zero temperature order parameter reported in $`d=1`$ where the CDW state is destroyed when the phonon frequency is of the order of the adiabatic ($`\gamma =0`$) mean-field gap value. We note that the scaled $`T_{\mathrm{co}}`$ curves fan out at larger phonon frequencies: in the extreme quantum limit the suppression of $`T_{\mathrm{co}}(\mathrm{\Omega })`$ relative to $`T_{\mathrm{co}}(\mathrm{\Omega }=0)`$ depends strongly on coupling. This is an indication that quantum non-Migdal effects (high energy electron self energy $`\mathrm{\Sigma }_n1/2\overline{\mathrm{\Lambda }}\overline{\mathrm{\Omega }}𝒢_n^0`$) are different from classical non-Migdal effects ($`\mathrm{\Sigma }_n\overline{\mathrm{\Lambda }}T𝒢_n^0`$).

$`T_{\mathrm{co}}(0)`$ increases considerably for spin $`s=1/2`$, essentially because the spin sum increases the bare susceptibility by a factor of two. Nevertheless, the change in spin $`s`$ cannot be accounted for by a change of $`\lambda `$ (i.e. $`\rho _0`$) only. Reducing $`\lambda `$ by a factor of two when changing the spin from $`s=0`$ to $`s=1/2`$ does not bring the transition temperature curves close to each other as shown in the inset of Fig. 2 (see the circles connected with a solid, $`s=0`$, or dashed, $`s=1/2`$, line). The rather different suppression of $`T_{\mathrm{co}}`$ for $`s=0`$ and $`s=1/2`$ was also found in studies in $`d=1`$. Even on the scale $`T_{\mathrm{co}}(0)`$ the large frequency tails are qualitatively different in the two cases. We provide some analytic insight in Section V.

Within the spin-1/2 Holstein model sizable transition temperatures might be obtained even in the antiadiabatic limit and for modest coupling strengths making the CDW transition accessible to QMC simulations. Freericks and Jarell report a suppression of $`T_{\mathrm{co}}/T_{\mathrm{co}}(0)0.25`$ for $`\mathrm{\Omega }/T_{\mathrm{co}}(0)7.1`$ (corresponding to $`\gamma =0.282`$). This result, obtained numerically in the antiadiabatic limit, indicates an even stronger suppression of $`T_{\mathrm{co}}`$ than obtained here in the adiabatic limit ($`\gamma `$ values used in our calculation vary from 0 to 0.08) but is of the same order. The comparison of absolute numbers is difficult because the work Ref.  does not use the semicircular density of states we have chosen.

It is important to note that $`T_{\mathrm{co}}(0)`$ is not the mean-field CDW transition temperature $`T_{\mathrm{co}}^{\mathrm{MF}}=(8t\gamma /e\pi )\mathrm{exp}[1/(2s+1)\lambda ]`$. Even at zero phonon frequency the charge density wave transition temperature is suppressed by thermal fluctuations of the phonon fields, e.g. $`T_{\mathrm{co}}(0)/T_{\mathrm{co}}^{\mathrm{MF}}=0.33,0.40,0.49`$ for $`s=0`$ and $`\lambda =0.258,0.204,0.156`$ respectively. In Fig. 3 $`T_{\mathrm{co}}`$ is plotted as function of the coupling constant $`\lambda `$. The suppression of the charge density wave transition temperature relative to its ‘BCS’ or static mean field value due to thermal fluctuations (white circle) and due to additional quantum fluctuations (white squares) is clearly visible. At weak coupling quantum fluctuations of the phonon fields are so effective in reducing $`T_{\mathrm{co}}`$ that it seems unlikely that CDW materials can be found with a Debye temperature $`\mathrm{\Theta }_\mathrm{D}`$ much larger than the mean field transition temperature $`T_{\mathrm{co}}^{\mathrm{MF}}`$. Higher order corrections in $`\overline{\gamma }`$ are accompanied by higher powers of $`\overline{\lambda }`$ and can no longer be neglected when the coupling constant is large. The comparison with the full numerical solution of the classical DMFT equation given in Ref.  reveals that the expansion suggested here breaks down already at the surprisingly small value of $`\overline{\gamma }\overline{\lambda }10^2`$. The inset of Fig. 3 will be discussed in Section VI.

The dependence of the CDW transition temperature on the ionic mass is usually discussed in terms of the isotope exponent $`\alpha =\frac{1}{2}\frac{\mathrm{\Omega }}{T_{\mathrm{co}}}\frac{\mathrm{d}T_{\mathrm{co}}}{\mathrm{d}\mathrm{\Omega }}`$ which is shown in Fig. 4. Because of the rapidly increasing factor $`\mathrm{\Omega }/T_{\mathrm{co}}`$ the negative isotope exponent $`\alpha `$ monotonically increases to very large values for spinless electrons. In the case $`s=1/2`$, however, the isotope exponent flattens off as function of $`\mathrm{\Omega }/t`$ due to the different large frequency behavior of $`T_{\mathrm{co}}`$.

## V Analytic expressions

In the adiabatic and extreme quantum limits we can derive explicit analytic expressions for $`T_{\mathrm{co}}`$. First, we consider the limit $`\mathrm{\Omega }/T_{\mathrm{co}}(\mathrm{\Omega }=0)0`$, i.e. we expand (25) around the classical solution obtained in Ref. . An appropriate Ansatz for the eigenvector $`ϵ_n`$ is

$$ϵ_n=𝒞\frac{\chi _n^{1/2}}{\left[1+(\overline{\lambda }T/\rho _0)\chi _n\right]}.$$
(30)
Note that, $`ϵ_n`$ is not quite the classical solution because it still depends on $`\overline{\gamma }`$ via $`\chi _n`$. Therefore, after inserting (30) into (25) we still have to expand in $`\overline{\gamma }`$ which finally leads to

$$\frac{1}{2s+1}=\frac{\lambda T}{\rho _o}\underset{n}{}\frac{\chi _n^{\mathrm{class}}}{1+(\overline{\lambda }T/\rho _0)\chi _n^{\mathrm{class}}}\lambda K\left(\frac{\mathrm{\Omega }}{T_{\mathrm{co}}(0)}\right)^2$$
(31)
with $`\chi _n^0=(\sqrt{(2t/\omega _n)^2+1}1)/2t^2`$. $`K`$ is a numerical coefficient independent of temperature but weakly dependent on $`\overline{\lambda }`$. Its magnitude is relatively small, e.g. $`K=1.775\times 10^2`$ for $`\lambda =0.156`$ ($`s=0`$).

The first part of Eq. (31) is the classical $`T_{\mathrm{co}}`$ equation. The sum has a logarithmic $`T`$ dependence. Therefore, expanding $`T=T_{\mathrm{co}}(0)\delta T`$ leads to a term $`[\lambda /T_{\mathrm{co}}(0)]\delta T`$. Consequently, the CDW transition temperature in the limit $`\mathrm{\Omega }/T_{\mathrm{co}}(0)1`$ is given by

$$\frac{T_{\mathrm{co}}}{T_{\mathrm{co}}(0)}=1K\left(\frac{\mathrm{\Omega }}{T_{\mathrm{co}}(0)}\right)^2.$$
(32)
The quadratic drop is clearly seen in Fig. 2. We can use Eq. (32) to derive the isotope exponent in the classical limit. We obtain

$$\alpha =K\left(\frac{\mathrm{\Omega }}{T_{\mathrm{co}}(0)}\right)^2=\frac{T_{\mathrm{co}}}{T_{\mathrm{co}}(0)}1.$$
(33)
The isotope exponent as function of $`t`$ according to this result is depicted in Fig. 4 as a dashed line. The exponential decrease of $`T_{\mathrm{co}}(0)`$ with increasing kinetic energy $`t`$ of the electrons is responsible for the sharp upturn of the negative isotope exponent.

In the extreme quantum limit the temperature is not a relevant high energy scale. Therefore, the upper cutoff frequency introduced above should be given by $`\omega _N=\eta \overline{\mathrm{\Omega }}`$ to yield results to logarithmic precision with $`\eta `$ being a numerical factor of order one. We can replace the phonon propagator by

$$\frac{1}{1+\stackrel{~}{\omega }_{nm}^2}1\stackrel{~}{\omega }_{nm}^2.$$
(34)
Moreover, to the same precision all Matsubara sums can be replaced by integrals and the temperature occurs only as a lower cutoff in the theory. The $`T_{\mathrm{co}}`$ equation (25) reads simply

$$1=\left(\lambda ^{}\overline{\lambda }\right)t_T^{\eta \overline{\mathrm{\Omega }}}\chi (\omega )d\omega .$$
(35)
For small frequencies the self energy (22) is given by $`\mathrm{\Sigma }_n=i\omega _n\overline{\lambda }`$. Therefore, $`\chi (\omega )=\chi _0(\omega )/(1+\overline{\lambda })`$ which, inserted in Eq. (35), gives the CDW transition temperature in the limit $`\mathrm{\Omega }/T_{\mathrm{co}}1`$

$$T_{\mathrm{co}}=\eta \overline{\mathrm{\Omega }}\mathrm{exp}\left(\frac{1+\overline{\lambda }}{\lambda ^{}\overline{\lambda }}\right)\eta \overline{\mathrm{\Omega }}\mathrm{exp}\left(\frac{1}{\lambda _{\mathrm{eff}}}\right).$$
(36)
Note that, $`\overline{\mathrm{\Omega }}`$ enters the calculation as the upper cutoff of the integral (35) and as lower cutoff in the effective interaction strength $`\lambda ^{}`$ (see (26)) in close resemblance to the well known McMillan formula for the superconducting transition temperature. Because the effective coupling constant $`\lambda _{\mathrm{eff}}`$ depends logarithmically on $`\overline{\mathrm{\Omega }}`$ (via $`\lambda ^{}`$) we obtain for the isotope exponent in the extreme quantum limit

$$\alpha =\frac{1}{2}\left[1\frac{1+\overline{\lambda }}{(1\overline{\lambda }/\lambda ^{})^2}\right].$$
(37)
The extreme quantum limit for fixed values of the phonon frequency $`\mathrm{\Omega }`$ is realized for small values of $`\lambda `$ (or large values of $`t`$) due to the fast drop of $`T_{\mathrm{co}}(0)`$ with decreasing coupling strength (see Fig. 3). For spinless electrons the isotope exponent becomes $`\alpha 1/(2\lambda ^2I^2)`$, i.e. $`\alpha `$ diverges in the limit $`\lambda 0`$ with power law behavior. The divergence of the negative isotope exponent with increasing kinetic energy of the electrons is shown in Fig. 4 as a dotted line. The isotope exponent for spin-1/2 electrons, however, remains finite, equalling $`\alpha =3/2`$ for $`\lambda 0`$. The difference between the cases $`s=0`$ and $`s=1/2`$ is due to a near-cancellation of the vertex parts for $`s=0`$. Consequently, the leading order of $`\lambda _{\mathrm{eff}}`$ is $`𝒪(\lambda ^2)`$ for $`s=0`$ and $`𝒪(\lambda )`$ for $`s=1/2`$.

In the ME theory for the superconducting transition temperature the effect of Coulomb interaction is considered by introducing a pseudopotential $`\mu ^{}`$. Therefore, it is tempting to absorb Coulomb effects on the charge density wave transition temperature in a similar way in $`\lambda ^{}`$. Indeed, a repulsive Coulomb interaction is an electron-electron interaction similar to the attractive vertex part, first term in Fig. 1(b), and will change the value of $`\lambda ^{}`$ simply by replacing $`(2s+1)\lambda (2s+1)\lambda U`$ (see also Ref. ). $`\lambda ^{}`$ acts as a pseudopotential with a non vanishing value for $`U=0`$. Nevertheless, Coulomb interaction will also change the electron self energy and the ‘vertex corrections’ (second and following terms in Fig. 1(b)) which is a more subtle effect due to the different frequency dependence of electron-phonon and Coulomb interactions and is left for future work.

## VI Relation to traditional approaches

It is traditional to write the CDW transition temperature in the form $`T_{\mathrm{co}}=W_{\mathrm{eff}}\mathrm{exp}\left(1/\lambda _{\mathrm{eff}}\right)`$ defining an effective bandwidth and coupling constant. The effective parameters $`W_{\mathrm{eff}}`$ and $`\lambda _{\mathrm{eff}}`$ are usually obtained by fitting to experimental data and therefore the functional dependence on the microscopic properties like electron-phonon coupling and adiabatic parameter remains hidden. The isotope exponent, however, explores this dependence and tests our physical understanding of the CDW transition.

A ‘BCS’ or static mean field approach would result in $`T_{\mathrm{co}}^{\mathrm{MF}}=2t\mathrm{exp}[1/(2s+1)\lambda ]`$, i.e. $`T_{\mathrm{co}}`$ depends exponentially on $`\lambda `$ but is independent of $`\gamma `$. For strong coupling the ‘BCS’ theory is clearly wrong because the Holstein model in this limit can be mapped on an effective pseudospin model and $`T_{\mathrm{co}}`$ has to decrease with increasing $`\lambda `$. Here we have demonstrated that also the weak coupling behavior is not predicted correctly and that instead of an ‘BCS’ isotope exponent $`\alpha =0`$ an enhanced ($`s=1/2`$) or even divergent ($`s=0`$) negative isotope effect is obtained.

Corrections to the BCS approximation have been considered. In the conventional Migdal-Eliashberg (ME) theory of electron-phonon effects in solids electron self energies arising from the interaction with the phonons are included but all vertex corrections are neglected. Other authors have argued that one should include also the feedback of the electrons on the phonon system but still neglecting the vertex corrections, thus defining a renormalized Migdal-Eliashberg theory (RME). In our notation this corresponds to dropping the second ($`\overline{\lambda }`$) term in Eq. (25). The RME approximation predicts an incorrect isotope effect because $`T_{\mathrm{co}}`$ depends on ion mass only via the dependence of $`\chi _m`$ on the normal state self energy, and comparison to our results shows that this dependence is not the leading one. At very weak coupling the RME theory reduces to the BCS theory (up to logarithmic corrections) and predicts $`\alpha =0`$. Freericks and co-workers added extra vertex corrections to the RME theory but these are higher order in $`\gamma `$ then the effects we studied in this paper.

We have plotted the charge density wave transition temperature for various coupling strengths $`\lambda `$ and fixed adiabatic parameter $`\gamma =0.1/\pi `$ and $`s=1/2`$ in the inset of Fig. 3 for both the RME approach and the present work. The charge density wave transition temperature obtained within RME theory is always larger than the one obtained from the systematic $`\overline{\gamma }`$ expansion presented here. A similar plot for a larger value of $`\gamma =0.5/\sqrt{\pi }`$ has been presented in Ref. , where the RME results were compared with QMC data, showing that the RME theory indeed overestimates the transition temperature for weak coupling strength at least in the antiadiabatic limit. Note that, the RME results presented in Ref.  are not precisely identical to ours because we have included the phonon self energy only to leading order in $`\overline{\gamma }`$, but we believe that this does not effect the conclusions drawn here. It is well known that $`T_{\mathrm{co}}(\lambda )`$ is maximized at a $`\lambda _{\mathrm{max}}<\mathrm{}`$. \[$`T_{\mathrm{co}}0`$ as $`\lambda 0`$ because the coupling vanishes; $`T_{\mathrm{co}}0`$ as $`\lambda \mathrm{}`$ because in this limit a picture of almost completely localized particles with intersite interactions $`t^2/\lambda `$ applies.\] The present approach underestimates both $`\lambda _{\mathrm{max}}`$ and $`T_{\mathrm{max}}`$, as may be seen (Fig. 3) by comparison to the classical treatment of Ref.  (to which the present treatment reduces as $`T_{\mathrm{co}}`$ becomes greater than $`\mathrm{\Omega }`$) and by comparison to QMC data (Ref. ) and indeed is not expected to apply for $`\lambda \lambda _{\mathrm{max}}`$.

A theory neglecting the renormalization of the phonon propagator by electrons predicts that $`T_{\mathrm{co}}`$ is maximized at $`\lambda =\mathrm{}`$. As an example we have plotted in the inset of Fig. 3 the conserving Hartree-Fock approximation discussed by Freericks et al. This corresponds to using Eq. (25) but with $`\overline{\lambda }`$ replaced by $`\lambda `$. In the weak coupling limit the difference between $`\overline{\lambda }`$ and $`\lambda `$ is negligible, so the CHF theory gives a good estimate of the CDW transition temperatures. The comparison with QMC data given in Ref.  demonstrate that the CHF theory, and therefore the theory given here, can (at least for $`s=1/2`$ and weak coupling strengths) be extended to the antiadiabatic limit.

Finally, we point out that the systematic expansion in $`\overline{\gamma }`$ discussed in this work reproduces the results of the renormalized ME theory when applied to the calculation of the superconducting transition temperature $`T_c`$. Following Freericks and Scalapino we find $`T_c\eta \overline{\mathrm{\Omega }}\mathrm{exp}[(1+\overline{\lambda })/\overline{\lambda }]`$ in the extreme quantum limit. For the half-filled Holstein model the charge density wave transition temperature is always larger than the superconducting transition temperature if $`\mathrm{\Omega }<\mathrm{}`$; only as $`\mathrm{\Omega }\mathrm{}`$ do they become equal. The superconducting transition temperature is controlled by a logarithmic divergence with upper cutoff set by the phonon frequency whereas the corresponding upper cutoff in the theory for $`T_{\mathrm{co}}`$ is enhanced towards $`2t`$ by the first diagram in Fig. 1(b).

## VII Conclusion and Comparison to Experiment

We have used the dynamical mean field method to obtain a systematic and tractable treatment of the isotope effect on charge density wave transitions. Our method is applicable to the physically relevant case in which both phonon frequency and transition temperature are small relative to the electronic bandwidth, and in combination with the results of Ref.  may be generalized to the crucial case in which other interactions beyond the electron-phonon one are important.

Quantum phonon effects are found to suppress the charge density wave transition temperature severely. This is in agreement with previous work on one dimensional models demonstrating that the suppression of $`T_{\mathrm{co}}`$ is not due to the strong quantum fluctuations characteristic of one spatial dimension. The isotope exponent of $`T_{\mathrm{co}}`$ for spin-1/2 electrons is enhanced compared to the ME prediction. In the case of spinless electrons the isotope exponent may even be divergent. Although beyond the scope of the theory, the consideration of phonon self energy effects lead to an optimal charge density wave transition temperature for values of the coupling strength slightly below the polaronic instability at $`\lambda _c`$ but peak position and height are underestimated.

Our results suggest that it will be difficult to find CDW systems with mean field transition temperatures much less than typical Debye frequencies. It is of interest to examine data with this in mind. For quasi-one-dimensional CDW materials the mean field CDW transition temperature is typically larger than room temperature and therefore larger than or of the same order as the Debye temperature $`\mathrm{\Theta }_\mathrm{D}`$. In $`\mathrm{TaS}_3`$, $`\mathrm{K}_{0.3}\mathrm{MoO}_3`$ and $`(\mathrm{TaSe}_4)_2\mathrm{I}`$ one deduces experimentally values of $`T_{\mathrm{co}}^{\mathrm{MF}}540K`$, $`T_{\mathrm{co}}^{\mathrm{MF}}330`$ and $`T_{\mathrm{co}}^{\mathrm{MF}}1040K`$ respectively. Quasi-two-dimensional CDW materials include the dichalcogenides and the molybdenum bronzes and oxides. For these materials it is more difficult to deduce mean field values of the CDW transition temperature because the CDW transition is not accompanied by a metal insulator transition. For $`2H\mathrm{TaSe}_2`$ one might argue from measurements of the optical conductivity for a value of $`T_{\mathrm{co}}^{\mathrm{MF}}580K`$ similar to the quasi-one-dimensional materials. Charge order transitions are also observed in many two and three dimensional transition metal oxides including some members of the high-temperature superconductor and ‘CMR’ family. Frequently, the Coulomb repulsion of the charge carriers is responsible for the charge order transition and therefore beyond the scope of this paper. Nevertheless, the phenomenology is very similar and e.g. for $`\mathrm{La}_{1.67}\mathrm{Sr}_{0.33}\mathrm{NiO}_4`$ one might deduce $`T_{\mathrm{co}}^{\mathrm{MF}}=750K`$. As was suggested by our theory the case of $`\mathrm{\Theta }_\mathrm{D}T_{\mathrm{co}}^{\mathrm{MF}}`$ seems to be extremely rare.

Our results bear at least a qualitative similarity to data obtained on ‘CMR’ manganites. Recently, a colossal isotope shift of the CDW transition temperature in $`\mathrm{Nd}_{0.5}\mathrm{Sr}_{0.5}\mathrm{MnO}_3`$ was reported. With magnetic field the negative isotope exponent varies from $`1.2`$ to $`4.5`$. Within the double exchange picture of the manganites an increase of the applied magnetic field means an increase of the kinetic energy $`t`$. Due to the large Hunds rule coupling of the conduction electron spins to the core spins of the Mn ions the conduction electrons are fully spin polarized. Therefore we should compare to the case $`s=0`$. To demonstrate the effect of an increasing $`t`$ in our calculation we draw a line for fixed $`\mathrm{\Lambda }`$ and $`\mathrm{\Omega }`$ in Fig. 4. Indeed, we find a monotonically increasing negative isotope effect with increasing $`t`$.

Acknowledgements. We thank A. Deppeler and M. Gugros for useful discussions. SB acknowledges the DFG, the Rutgers University Center for Materials Theory and NSF DMR0081075 for financial support; AJM acknowledges NSF DMR0081075.


text:
# Ab Initio Molecular Dynamics Simulation of Liquid GaxAs1-x Alloys

## I Introduction

The study of liquid metals and alloys has drawn considerable attention recently, in particular due to the possibility of carrying out first principles or ab initio calculations for these systems. In these calculations the electronic structure is evaluated quantum mechanically using Density Functional Theory (DFT) and the corresponding forces are used to move the ions according to classical molecular dynamics. Using this approach it is possible to calculate both the atomic and electronic structure consistently and to see how changes in one are correlated with changes in the other. By now, a number of groups have used these methods to calculate both the thermodynamic and transport properties of a variety of liquid metals and alloys.

Liquid semiconductors are of particular interest from the point of view of ab initio calculations. By “liquid semiconductors,” we mean liquids of materials which are semiconducting in their solid phases, such as Si, Ge, GaAs, and CdTe. Somewhat surprisingly, most of these are reasonably good metals in their liquid phases. For example, Si, Ge, and GaAs all have conductivities near melting which lie in the metallic range, and which tend to decrease with increasing temperatures, as is characteristic of metals. This metallic behavior is correlated with an increase in coordination number on melting, the liquid is thus more close-packed than the solid and has a higher density. By contrast, $`\mathrm{}`$-CdTe is poorly conducting in its liquid state and its conductivity increases with increasing temperature, characteristic of semiconductors.

Recent ab initio calculations for several of these materials give behavior which is in good agreement with these experiments. Godlevsky et al have found, in agreement with experiment, that stoichiometric GaAs is metallic, whereas stoichiometric CdTe is a reasonable insulator. These differences in the electronic properties were related to the differences in the structural properties occurring within the melt. An earlier calculation by Zhang et al studied stoichiometric $`\mathrm{}`$-GaAs using the Car-Parrinello version of ab initio molecular dynamics. In this calculation too, it was found that $`\mathrm{}`$-GaAs is a metallic, weakly ionic liquid, with a larger coordination number than the insulating solid phase.

In this paper, we describe a numerical study of $`\mathrm{}`$-GaxAs1-x over a range of concentrations, using ab initio techniques. The liquid is miscible over the whole concentration range, unlike the solid, which exists only at stoichiometry. Such a study is of interest for a variety of reasons. The properties of of $`\mathrm{}`$-Ga and $`\mathrm{}`$-As stand in marked contrast to each other. $`\mathrm{}`$-Ga is a close-packed liquid metal with a coordination number of $``$ 9 and its electronic density of states is almost free-electron like. In contrast, $`\mathrm{}`$-As has the same coordination number as in the crystalline phase ($``$ 3), and is a narrow band-gap semiconductor in the liquid state. Thus while we expect $`\mathrm{}`$-Ga to show metallic bonding, $`\mathrm{}`$-As is expected to retain the covalent character of the bond upon melting. Previous ab initio calculations for these liquids  have indeed confirmed this picture. Thus as the stoichiometry is varied for $`\mathrm{}`$-GaxAs1-x, we would expect interesting changes both structurally and electronically. It is also of interest to compare the structures of $`\mathrm{}`$-GaxAs1-x for $`x=0.2`$ and $`x=0.8`$ with the structure for the corresponding pure liquids (As and Ga respectively) to see how changes in structural properties are correlated with changes in electronic properties.

We turn now to the body of the paper. In section II, we briefly summarize our approach and method of calculation. Our results are presented extensively in Section III, together with some analysis connecting the calculations to a qualitative picture of the liquid state. Finally, in Section IV, we give a short concluding discussion.

## II Method and Computational details

Our method of carrying out the ab initio simulations is similar to that of our previous work. We use the plane-wave pseudopotential approach with generalized norm-conserving pseudopotentials in the Kleinman-Bylander form, treating the $`d`$-wave part as the local component. The exchange-correlation potential is computed using the local-density approximation (LDA), using the Ceperley-Alder result for the exchange-correlation energy, as parametrized by Perdew and Zunger. The details of the code can be found in the literature

Our liquid-state molecular dynamics (MD) simulations were carried out at a temperature $`T=1600`$ K, just above the melting point of GaAs ($`T=1515`$ K). We have considered five concentrations of GaxAs1-x: x = 0.2, 0.4, 0.5, 0.6, and 0.8. In each case, we used a cubic 64-atom supercell with periodic boundary conditions. For this size cell, the actual numbers of Ga atoms in the five samples were 13, 26, 32, 38, and 51. (Since the minority component in the 20% and 80% samples, the statistics will be quite poor for them and the results should be given greatest significance for the majority components in those cases.) The atomic densities for the five concentrations were obtained from the measured density of $`\mathrm{}`$-GaAs at this temperature, together with Vegard’s Law (i. e. linear interpolation of atomic volumes). We use a 10-Ry cutoff for the energies of the plane waves included in the wave function expansion, and $`\mathrm{\Gamma }`$-point sampling for the supercell Brillouin zone. For the electronic structure, we used Fermi-surface broadening corresponding to an electronic subsystem temperature of $`k_BT^{el}=0.1`$eV. In calculating the electronic wave functions, at each concentration we include eight empty bands. We control the ionic temperature using the Nosé-Hoover thermostat . The equations of motion are integrated by means of the Verlet algorithm, using an ionic time step of 125 a.u. ($``$ 3 fs). For each concentration, the samples were equilibrated for about 0.2ps, following which simulations were carried out for more than 3ps.

These simulations start from an initial configuration which, for the stoichiometric liquid, was generated from a classical molecular-dynamics (CMD) simulation using potentials of the Stillinger-Weber form. Although our particular potentials were originally derived for liquid GexSi1-x alloys, they should give a reasonable starting configuration for the Ga1-xAsx alloys which we study. We use these potentials but with Ga and As masses rather than those of Si and Ge. This model system was melted using classical molecular dynamics at high temperatures, then cooled down gradually to the temperature of interest. The configuration thus obtained from CMD was equilibrated for 0.2 ps (about double the expected relaxation time for this system). For the nonstoichiometric liquids, we used the same starting configuration as in the stoichiometric case, but with Ga atoms randomly substituted for the appropriate number of As atoms (or vice versa), so as to give the correct concentration.

## III Results

Fig. 1 shows the three partial pair correlation functions gGaGa(r), gAsAs(r) and gGaAs for liquid GaxAs1-x at the concentration $`x=0.5`$ and a temperature T = 1600 K. A number of features deserve mention. First, at $`x=0.5`$, the principal peaks of all three partial pair correlation functions occur at about the same separation, namely 2.5Å. This indicates the non-ionic character of the bonds; in ionic liquids the partial pair correlation functions for like atoms are out of phase with the corresponding function for unlike atoms . Our results for gαβ(r) are in close agreement with the recent calculation by Godlevsky et al ; an earlier calculation by Zhang et al  also gives the same features although these authors get a much stronger principal peak for gGaAs(r) than is found in either our calculations of that of Ref. . While the principal peaks are in phase, there are some differences among the partial g(r)’s. For example, at $`x=0.5`$, gAsAs(r) has a slightly higher and narrower first peak, and a stronger second peak, than gGaGa(r), while gGaGa(r) has a broad first peak and no obvious peak beyond that. We have also calculated the coordination numbers for the first shell of neighbors, defined as the integral of $`4\pi \rho r^2`$g(r) from zero out to the minimum after the first maximum in g(r). \[Here g(r) is the total pair correlation function, which does not distinguish between the two species, normalized so that it approaches unity at large r.\] The coordination we calculate in this way at $`T=1600`$ K is 5.8, in good agreement with the experimental estimate of 5.5 $`\pm `$ 0.5. Note also that this value is larger than the diamond-structure value of $`4`$ but significantly smaller than the value expected in a close-packed liquid, which would be in the range of 9 or 10. This value indicates the persistence of covalent bonding in $`\mathrm{}`$-GaAs.

Fig. 2 shows gGaGa(r) at $`x=0.8`$ and gAsAs(r) at $`x=0.2`$. The latter shows more short-range order than the former - specifically, a sharper main peak and a broad second peak, rather than a single broad principal peak. We believe that there are several causes for these differences. First, pure Ga has a much lower melting temperature than either Ga0.5As0.5 or As. Thus, at the same temperature of 1600 K, we expect less short-range order for gGaGa(r) at $`x=0.8`$ than for gAsAs(r) at $`x=0.2`$, as seen in our calculations. In addition, we expect the $`x=0.2`$ sample to show some residue of the complex local structure seen in pure $`\mathrm{}`$-As, which is, in turn, quite similar to that of the crystalline phase. Indeed, our calculated g$`{}_{AsAs}{}^{}(r)`$ at $`x=0.2`$ has some of the same features as those in the calculated g(r) for pure $`\mathrm{}`$-As. However, there are also some observable differences, which may be related to the fact that, experimentally, pure $`\mathrm{}`$-As is semiconducting while $`\mathrm{}`$-Ga0.2As0.8 is calculated to be metallic (as shown later when we calculate the electronic structure). First, if we integrate gAsAs(r) for the $`x=0.2`$ sample out to the first minimum beyond the main peak, we obtain a coordination number of $`3.2`$ As atoms for the first shell of nearest neighbors surrounding an As atom. This is slightly larger than the value $`3`$ reported for pure $`\mathrm{}`$-As (at T=1150 K) experimentally and in previous calculations. A more important difference is that pure $`\mathrm{}`$-As has a sharp second peak in g(r), whereas for gAsAs(r) in our simulations this second peak is very broad. While some of this broadening may result from the higher temperature ($`T=1600K`$) in our simulations for $`x=0.2`$, these results suggest that the local structure seen in $`\mathrm{}`$-As is preserved only out to the first shell of neighbors in $`\mathrm{}`$-Ga0.2As0.8.

On the other hand, the features we see for gGaGa(r) at $`x=0.8`$ are qualitatively similar and consistent with those seen for liquid Ga at lower temperatures . Using the procedure indicated above, we get a Ga coordination number of 5.8 for the first shell. If we integrate the total g(r) up to the first shell (i. e., including both As and Ga neighbors of Ga), we obtain a Ga coordination number of 6.9. While these values are smaller than the coordination number of 9.0 reported for pure $`\mathrm{}`$-Ga, we attribute the difference to the lowering of the first peak due to the higher temperature of our simulations ($`T=1600`$ K) compared to those for pure $`\mathrm{}`$-Ga ($`T=702`$ K and $`T=982`$ K) .

Figs. 3-5 show information about the various partial alloy structure factors at the same temperature. The partial structure factors $`S_{ij}(𝐤)`$ are defined in one of the standard ways:

$$S_{ij}(𝐤)=(N_iN_j)^{1/2}\underset{i}{}\underset{j}{}e^{i𝐤(𝐑_i𝐑_j)}(N_iN_j)^{1/2}\delta _{𝐤,0}$$
(1)
where $`i`$ and $`j`$ denote the two components of the binary alloy. Fig. 3 shows the calculated total structure factor S(k), as weighted by neutron scattering factors at a concentration $`x=0.5`$. S(k) is defined by

$$S(k)=\frac{b_{i}^{}{}_{}{}^{2}S_{ii}(k)+2b_ib_jS_{ij}(k)+b_{j}^{}{}_{}{}^{2}S_{jj}(k)}{b_{i}^{}{}_{}{}^{2}+b_{j}^{}{}_{}{}^{2}}$$
(2)
where $`b_i`$ and $`b_j`$ are the corresponding experimental neutron-scattering lengths ( $`b_{Ga}`$ = 7.2 and $`b_{As}=6.7`$) For comparison, we also show in Fig. 3 the quantity S(k) as measured by neutron diffraction. As can be seen from the Figure, the two agree quite well. In particular, the calculations convincingly reproduce the experimentally observed shoulder on the high-$`k`$ side of the principal peak in S(k).

Fig. 4 shows the three partial alloy structure factors SGaGa(k), SAsAs(k), and SGaAs(k) at x = 0.5 and T = 1600K. The structure factor for like pairs is always positive, with a conspicuous first peak, while the structure factor between opposite pairs is negative for small k, becoming positive at k values corresponding to the peaks in the other two partial structure factors. It is of interest to compare these results to those found in other model calculations. For a mixture of hard spheres of packing fraction 0.45 (characteristic of the liquid near melting) and ratio of hard sphere diameters of 0.9, the cross-correlation function is negative at small k and has a peak near that of the two same-species functions, as in our calculations. By contrast, for liquid NaCl, which is strongly ionic), the cross-correlation function has a strong negative peak at the same k as the peaks of the same-species partial structure factors. Thus, our results are more similar to the hard-sphere structure factors, suggesting that liquid GaAs is at most only weakly ionic.

Finally, Fig. 5 shows SGaGa(k) at $`x=0.8`$ and SAsAs(k) at $`x=0.2`$. Once again, like the real-space correlation functions at the same concentrations, the Ga-Ga structure factors show slightly less correlation (i. e., a lower principal peak and a less conspicuous second peak) than do the corresponding As-As structure factors. One possible reason for this behavior, as for the corresponding real-space correlation functions, is that the $`x=0.8`$ sample is further from melting than is the $`x=0.2`$ liquid.

SAsAs(k) for $`x=0.2`$ shows characteristic differences from that for pure $`\mathrm{}`$-As, which are analogous to those discussed earlier for the g(r)’s. Specifically, S(k) for pure $`\mathrm{}`$-As has a split principal peak with maxima at $`k=2.45`$ $`\AA ^1`$ and $`k=3.74`$ $`\AA ^1`$. By contrast, SAsAs(k) for $`x=0.2`$ has a peak at 2.5 $`\AA `$, but the second peak is reduced to only a shoulder at about 3.5 $`\AA `$. The fact that the second of the split peaks is smoothed to a shoulder at $`x=0.2`$ indicates a change in the local structure which is also reflected in the reduction of the second peak of partial g(r) as noted earlier.

We next discuss the $`x=0.8`$ sample, comparing our results with those for pure $`\mathrm{}`$-Ga at lower temperatures . The partial structure factor show the same qualitative features as the structure factors for the pure liquid, but there are some quantitative differences. Our results for $`x=0.8`$ (at $`1600`$ K) show a first peak in SGaGa(k) with a maximum of only about 1.3,which is lower than the experimental one seen in pure Ga at $`959`$ K ($``$ 1.7). We attribute the lowering of the first peak to the increased temperature, as seen experimentally in most liquid metals.

Further information about the short-range order in the liquid alloy may be obtained from bond angle distribution functions, shown in Figs. 6 and 7. These functions are defined by analogy with our previous work in liquid Ge. Namely, one considers a group of three atoms. Of these, one is denoted as the central atom; the other two atoms (denoted as “side atoms”), with the central atom, define a bond angle $`\theta `$. $`g_3(\theta ,r_c)`$ is the distribution of bond angles formed by all such groups of three atoms, such that both the side atoms lie within a cutoff distance $`r_c`$ of the central atom. Fig. 6 shows $`g_3(\theta ,r_c)`$ for the Ga-Ga-Ga angles at $`x=0.8`$, and for the As-As-As angles at $`x=0.2`$, each for two different choices of the cutoff radius $`r_c`$. Fig. 7 shows the same functions for Ga-Ga-Ga and for As-As-As at $`x=0.5`$.

The peaks in these distribution functions give hints about about the short-range bond-order in the liquids. For example, a peak near $`\theta =60^o`$ corresponds to a relatively close-packed arrangement of the corresponding atomic group, with many nearest neighbors. In contrast, a peak near $`100^o`$ indicates a more tetrahedral structure, typical of covalent bonding. Thus, the upper part of Figs. 6 and 7 suggest that the Ga ions form a rather close-packed arrangement at $`x=0.8`$ and $`x=0.5`$, since there is a strong peak near $`60^o`$ for both values of the cutoff $`r_c`$. By contrast, the lower parts of Figs. 6 and 7 suggest that, as expected, the As atoms have a more open arrangement at $`x=0.2`$, since there are two peaks in $`g_3`$ at this concentration: a side peak near $`50^o`$, and another noticeable peak near $`97^o`$. We also observe the peak at $`97^o`$ for As atoms at $`x=0.5`$, however it is broader and less pronounced than at $`x=0.2`$. In both Figs. 6 and 7, the Ga-Ga-Ga bond angle distribution depends little on $`r_c`$.

At $`x=0.2`$, the As-As-As distributions show noticeable $`5060^o`$ peaks only at the larger cutoff values. At the smaller cutoff radius, the $`5060^o`$ peak is missing. Thus, at short distances, the As clusters tend to maintain the local version of the structure they have in the pure liquid phase (and the crystalline phase) which shows a strong peak at $`97^o`$, but at larger cutoffs the local structure differs from that of pure $`\mathrm{}`$-As. We have made the same observation in connection with gAsAs(r) at the same concentrations, and with the shoulder in SAsAs(k). The $`97^o`$ peak implies some tetrahedral order persisting to $`x=0.2`$, though this peak is less pronounced than in pure $`\mathrm{}`$-As.

We have also calculated the electronic properties of $`\mathrm{}`$-GaxAs1-x. We calculate the single-particle electronic density of states $`N(E)`$ in the standard way, by using the expression

$$N(E)=\underset{𝐤,E_𝐤}{}w_𝐤g(EE_𝐤).$$
(3)
In this expression $`E_𝐤`$ denotes one of the energy eigenvalues of the single-particle wave functions at a particular $`𝐤`$ point within the supercell Brillouin zone, wk is the weight of that $`𝐤`$ point (defined below), and g(E) is a Gaussian smoothing function of width $`\sigma =0.2eV`$. Our calculation is carried out by sampling the supercell Brillouin zone at eight special $`𝐤`$ points, using the same choice of special points and weights as that of Holender et al in their calculations for pure $`\mathrm{}`$-Ga, which have been well tested and found to be an adequate representation of the supercell Brillouin zone. This choice is, however, convenient rather than unique; we expect that other choice having the same number of k points would have given similar results, as has been found in pure $`\mathrm{}`$-Ga. For each $`𝐤`$-point we include 40 conduction band states, and for each concentration, we obtain our final results by averaging over twelve representative liquid state configurations.

The resulting calculated density of states N(E) is shown for the four concentrations $`x=0.2,0.4,0.6`$, and $`0.8`$ in Fig. 8. \[We do not show our calculated N(E) for $`x=0.5`$, but it interpolates smoothly between $`x=0.4`$ and $`x=0.6`$.\] As in our previous studies, the alloy has a clearly metallic density of states for all concentrations $`x`$. However, just as in our previous results for GaxGe1-x, the density of states becomes more and more free-electron like as the concentration $`x`$ of the metallic component (Ga in this case) increases. Pure $`\mathrm{}`$-As is semiconducting and has been calculated to have deep minima in the electronic density of states at the Fermi energy ($`E_F`$), and also at an energy of $``$ -7 eV (measured from $`E_F`$) , which separates the $`s`$ and $`p`$ bands. Liquid Ga, on the other hand, has an almost free-electron like density of states. We see these features reflected in our simulations; as the figures indicate, the electronic density of states (DOS) has a pseudogap in the As rich phase which progressively fills up as the Ga concentration is increased, so that for $`x=0.8`$ it is hardly noticeable. But even at low Ga concentration ($`x=0.2`$), there is no minimum in the density of states at the Fermi energy. As for the pseudogap, we find that its position changes monotonically to lower energies (relative to $`E_F`$) with increasing Ga concentration. A similar pseudogap is reported in calculations for pure $`\mathrm{}`$-Ge , for the same reason, i. e., a partial separation between $`s`$-like and $`p`$-like bands.

We have also computed the frequency-dependent electrical conductivity $`\sigma (\omega )`$ for $`x=0.2`$, $`0.4`$, $`0.5`$, $`0.6`$, and $`0.8`$, at frequencies ranging up to 2eV. $`\sigma (\omega )`$ is given by the standard Kubo-Greenwood expression

$$\sigma (\omega )=\frac{2\pi e^2}{3m^2\omega \mathrm{\Omega }}\underset{i}{}\underset{j}{}\underset{\alpha }{}(f_jf_i)|\psi _i|\widehat{p_\alpha }|\psi _j|^2\delta (E_jE_i\mathrm{}\omega ).$$
(4)
where $`\psi _i`$ and $`\psi _j`$ are the single particle Kohn-Sham wave functions with Fermi occupancies $`f_i`$ and $`f_j`$ and energy eigenvalues $`E_i`$ and $`E_j`$. Once again, we calculate the conductivity using the same set of eight special k points used for N(E), and averaging over the same twelve representative ionic configurations, including 40 conducting band states for each $`𝐤`$. While there is some inaccuracy in using differences in Kohn-Sham eigenvalues, as we do in this equation, it is difficult to do better at present in liquid metals and semiconductors, and this same approach does give good agreement with experiment in a number of other liquid metals, especially at zero frequency but also for stoichiometric $`\mathrm{}`$-GaAs at finite frequencies.

The frequency-dependent conductivity is shown in Fig. 9 for $`x=0.2`$ and $`x=0.8`$, and its calculated zero-frequency limit is listed in Table I for the five concentrations x = 0.2, 0.4, 0.5, 0.6 and 0.8, all at $`T=1600K`$.

Several features of the conductivity graphs, and of the tabulated d. c. limits, deserve mention. First, the calculated value of the d. c. conductivity at x = 0.5 is very close to the experimental value: $`0.84\times 10^4`$ ohm-1cm-1, compared to the experimental value of $`0.79\times 10^4`$ ohm-1cm-1. Secondly, the conductivity has a weak minimum near x = 0.5. This is consistent with expectations based on second-order perturbation theory which would predict that alloy scattering (due to concentration fluctuations) would be a maximum near x = 0.5 Thirdly, the frequency-dependence of $`\sigma (\omega )`$ becomes more metallic as $`x`$ increases. Specifically, at x = 0.8, $`\sigma (\omega )`$ clearly decreases with increasing frequency, characteristic of a Drude metal, while for the highest As concentration (x = 0.2), the conductivity is nearly frequency-independent in the range of calculation. This behavior is closer to the nonmetallic behavior in which the conductivity increases with increasing frequency.

Finally, we have computed one important atomic transport coefficient, namely, the atomic self-diffusion coefficients $`D_{ii}`$ for the majority species in the two liquids Ga0.2As0.8 and Ga0.8As0.2. In both cases, the $`D_{ii}`$’s can be extracted from a plot of the mean-square atomic displacement versus time, which approaches a straight line in the limit of large time. The expression is

$$D_{ii}=\underset{t\mathrm{}}{lim}\frac{|𝐑_i(t)𝐑_i(0)|^2}{6t},$$
(5)
where $`𝐑_i(t)`$ denotes the position of an ion of species $`i`$ at time $`t`$, and the triangular brackets denote an average over all atoms of species $`i`$ and over initial times. Plots of the mean-square atomic displacement as a function of time are shown in Fig. 10 for both types of atoms at $`x=0.2`$ and $`0.8`$ as indicated above; ¿From these we obtain $`D_{Ga}=2.1\times 10^4`$ cm2/sec at $`x=0.8`$; $`D_{As}=1.7\times 10^4`$ cm2/sec at $`x=0.2`$.

The value for $`D_{As}`$ at $`x=0.2`$ is about three times larger than that obtained by Li  for pure $`\mathrm{}`$-As ($`D0.6\times 10^4`$ cm2/sec). This difference is probably due to a combination of several factors. First, the calculations by Li are carried out at a temperature $``$ 450 K lower than ours. Secondly, pure $`\mathrm{}`$-As seems to have more covalent bonding than Ga0.2As0.8, which probably impedes atomic motion, giving a lower atomic diffusion coefficient for As. Thus, in short, this behavior seems to be consistent with the rest of our picture, which is that $`\mathrm{}`$-GaxAs1-x rapidly acquires metallic conductivity, and corresponding structural properties for $`x`$ as small as 0.2.

## IV Discussion and Conclusions

The most striking results of these calculations is that GaxAs1-x remains metallic at all concentrations between $`x=0.2`$ and $`x=0.8`$, including the As rich value $`x=0.2`$. Thus even low Ga concentrations are sufficient to render $`\mathrm{}`$-GaxAs1-x metallic (recall that pure $`\mathrm{}`$-As is semiconducting). This is reflected in the electronic density of states which shows no minimum at the Fermi energy for all concentrations studied. The liquid structure is also consistent with metallic behavior at all concentrations between $`x=0.2`$ and $`0.8`$, although there are some clear deviations from the behavior seen in simple metallic alloys. Specifically, although the coordination number at all concentrations is larger than the value of four that might be expected in a predominantly covalent liquid, it is still smaller than that of a typical close-packed hard-sphere mixture. At $`x=0.2`$ and $`0.8`$, the pair correlation functions and structure factors resemble those of the corresponding pure liquids, except that the split first peak in $`S(k)`$ of $`\mathrm{}`$-As becomes a single peak with a weak shoulder in As.

It has been suggested  that semiconducting properties persist in a liquid only if the liquid has the same short-range order as the crystalline phase. In the alloys we study, the liquid state has a short range order which is distinctly different from the solid. We also find that all these alloys are metallic, in agreement with the suggestion of Ref. . By contrast, again in agreement with the picture advanced in , another compound semiconductor, stoichiometric CdTe, appears to preserve the crystalline short-range order, and also to retain its semiconducting characteristics in the liquid state. Such semiconducting behavior was indeed found in the ab initio calculation by Godlevsky et al. It would be of interest to extend their calculations off stoichiometry, where metallic behavior is likely.

The other characteristics of $`\mathrm{}`$-GaxAs1-x reflect its fundamentally liquid-metallic character. For example, the resistivity is predicted to exhibit typical Faber-Ziman behavior: a weak positive deviation from a linear concentration dependence near stoichiometry, which is caused by alloy scattering. The calculated values of the partial atomic diffusion coefficients are comparable to the diffusion coefficients of those other liquid semiconductors which are metallic in their liquid states. We find no evidence of a strong reduction in this value because of formation of clusters near stoichiometry; such cluster formation is not expected for GaxAs1-x because of the small electronegativity differences between the two species. One might speculate that, in other liquid semiconductors, such as stoichiometric $`\mathrm{}`$-CdTe, which remain poorly conducting in the liquid state, the local structure is much more ionic near melting and the atomic diffusion coefficients are correspondingly lower.

In summary, our calculations show that GaxAs1-x is a reasonable metal at all concentrations between $`x=0.2`$ and $`x=0.8`$. In particular, there is no evidence of strong compound formation in the liquid state near $`x=0.5`$. The electrical conductivity shows a concentration dependence typical of a liquid metallic alloy, with evidence of weak scattering from concentration fluctuations which reaches a maximum near $`x=0.5`$. The electronic density of states shows no minimum at the Fermi energy; instead, it has a pseudogap between $`s`$-like and $`p`$-like occupied state which persists at all concentrations, though it is considerably weaker in the Ga-rich alloys. The atomic diffusion coefficient is calculated to be similar to that of other liquid semiconductors which are metallic in their liquid state. Finally, the liquid structure shows some indications of deviation from the behavior of simple liquid metal alloys. The principal evidence of deviation from behavior characteristic of a simple liquid alloy formation comes from the calculated pair correlation functions, structure factors, and bond angle distribution functions, all of which show some weak indications of departure from close-packed behavior: smaller coordination functions than in typical hard-sphere liquids, and a weak residue of tetrahedral bonding.

## V Acknowledgements

We are grateful for support by NASA, Division of Microgravity Sciences, under grants number NCC3-555 and NCC8-152. We also thank M. Scheffler and the research group at the Fritz-Haber Institute (Theory Department) for use of their ab initio code FHI96MD, and D. Matthiesen, A. Chait, and E. Almaas for many valuable conversations. These calculations were carried out using the SP2 at the Ohio Supercomputer Center and the RS/6000 cluster at the NASA Lewis Research Center.

FIGURE CAPTIONS

1. (a) Pair correlation functions gGaGa(r), gAsAs(r), and gGaAs(r) for $`\mathrm{}`$-GaxAs1-x at a temperature of 1600 K and concentration $`x=0.5`$. The graphs are vertically offset by one unit each for clarity.
2. (top) gGaGa(r) for x = 0.8; (bottom) gAsAs(r) at x = 0.2 and $`T=1600`$ K.
3. Full line: Calculated neutron structure factor S(k) = $`_{i,j=Ga}^{As}`$fij(k)Sij(k), where $`f_{ij}(k)`$ is the neutron form factor, for $`x=0.5`$. Open circles: measured S(k) as obtained by neutron diffraction \[ref. \].
4. Partial structure factors SGaGa(k), SAsAs(k), and SGaAs(k) for $`x=0.5`$. SAsAs(k) is vertically offset for clarity.
5. Calculated partial structure factors (top) SGaGa(k) at x = 0.8; (bottom) SAsAs(k) at x = 0.2.
6. Calculated bond angle distribution functions $`g(\theta ,r_c)`$ for (top) groups of three Ga atoms at a concentration $`x=0.8`$, and (bottom) groups of three As atoms at a concentration $`x=0.2`$, both for two different cutoff radii, $`r_c=3.4\AA `$ and $`r_c=3.8\AA `$, as defined in the text.
7. Same as for Fig. 6 but for x = 0.5 and (top) g$`{}_{\text{GaGaGa}}{}^{}(\theta ,r_c)`$; (bottom) g$`{}_{\text{AsAsAs}}{}^{}(\theta ,r_c)`$;
8. Single-particle electronic density of states $`N(E)`$ for GaxAs1-x at $`x=0.2`$, $`0.4`$, $`0.6`$, and $`0.8`$. The Fermi energy in each case is shown as a dashed vertical line.
9. Calculated electrical conductivity $`\sigma (\omega )`$ for GaxAs1-x at $`x=0.2`$ and $`x=0.8`$.
10. Calculated mean-square displacement $`|𝐑_i(0)𝐑_i(t)|^2`$, plotted as a function of $`t`$ at x = 0.2 for As and at x = 0.8 for Ga.