Supercell band calculations and correlation for high-$T_C$ copper oxide superconductors

First principle band calculations based on local versions of density functional theory (DFT), together with results from nearly free-electron models, can describe many typical but unusual properties of the high-$T_C$ copper oxides. The methods and a few of the most important results are reviewed. Some additional calculations are presented and the problems with the commonly used approximate versions of DFT for oxides are discussed with a few ideas for corrections. It is concluded that rather modest corrections to the approximate DFT, without particular assumptions about strong correlation, can push the ground state towards anti-ferro magnetic (AFM) order. Spin fluctuations interacting with phonons are crucial for the mechanism of superconductivity in this scenario.

It is often assumed that strong correlation is important for an understanding of the high-T C problem [1,2,3]. The failure of the local approximations to the exchangecorrelation functional, such as the local (spin) density approximation, LDA or LSDA [4,5,6,7] to produce the anti-ferro magnetic (AFM) insulating state of many undoped transition metal oxides is generally quoted as the essential reason for discarding DFT calculations for high-T C copper oxides [3]. DFT is essentially exact, but the approximations to make it practical for use in real applications seem inappropriate for oxides. However, despite this problem there are several LDA results that fit to the observed high-T C properties, provided that doping and supercells are considered in order to account for imperfect lattice conditions such as stripes, phonons or spin waves [8,9,10]. Here is presented a short review of those DFT results together with a discussion of the correlation problem and an attempt to include further corrections to LSDA. Estimations of couplings λ caused by spin-phonon coupling (SPC) and pure spin fluctuations are also made.

II. NEARLY FREE-ELECTRON MODEL.
The spin polarized potential from an AFM order (like that of undoped La 2 CuO 4 ) can be generated from V AF M (x) = V 0 exp(−i Q · x), where Q = π/a 0 , is at the Brillouin Zone (BZ) boundary, so that there is one Cu site with spin up and one with spin down within a distance a 0 (a 0 is the lattice constant). The potential of the other spin has a phase shift of π. The band dispersion from a 1-dimensional (1D) nearly free-electron model (NFE) is obtained from a 2x2 eigenvalue problem, and it has a gap of 2V 0 at the zone boundary. Suppose now that an additional potential modulation V q (x) = V q exp(i q · x) exists, and that | q |≪| Q |, so that its periodicity or "wavelength" is much larger than 2a 0 . The product of these two modulations can be written V (x) = V q exp(−i( Q − q) · x), where the two amplitudes are combined into one coefficient, V q . Fig. 1 shows the real space configurations along the 100-direction. This potential describes AFM modulations in stripes, with magnetic nodes (Cu with zero moment) separated by a 0 Q/q. The values of V q for phonons or spin waves of different lengths are not known from the NFE model, but some values will be fed in from the ab-initio Linear Muffin-Tin Orbital (LMTO, [11]) calculations, as will be described below. The LMTO calculations are slow for large supercells. The cells need to be large to cover the periodicity of realistic phonon and/or spin-waves, and so far it has been possible to extend the cells in one direction only, usually the CuO bond direction. In contrast, the NFE model is very simple and it can easily be extended to two dimensions (2D). The band dispersion for k-points (k x , k y ) is now obtained from the eigenvalues of a 3x3 ma- trix [12]. Not only stripe order, but also "checkerboard" configuration can be modeled.
Examples of the density-of-states (DOS) for 3 different V q are shown in figure 2. The potential modulations are parallel to the CuO bond directions along x and y, so the gaps appear on the (k x , 0) and (0, k y ) lines, as shown in figure 3. However, almost no effect from V q appears on the band dispersion in the diagonal direction (along (k, k)). Therefore, the Fermi surface (FS) is not affected along the diagonal, but fluctuations, described through amplitude variations of V q for different positions and different time, will make the FS smeared in the two bond directions. An example of this is shown in figure 4.

III. AB-INITIO 1D-LMTO.
By 1D (1-dimensional) LMTO we mean ab-initio LMTO calculations for long (and narrow) supercells most often oriented along the CuO bond direction [13,14]. These calculations are based on the local version of DFT, the local (spin) density approximation, LDA or LSDA [4,5]. Cells with phonon distortions and/or spin waves within lengths of 4,8 and 12 lattice constants are typically considered in these calculations. The wavelengths of spin waves are twice as long as those of phonons. Hole doping, h, in La (2−h) Ba h CuO 4 (LBCO) is generally modeled by the virtual crystal approximation (VCA) where the nuclear and electronic La-charges (57.0) are reduced to (57-h/2) to account for a perfectly delocalized doping (h in holes per Cu). But some calculations with real La/Ba substitutions show that improved superconducting properties can be expected from periodic doping distributions [15]. The maximal phonon distortion, u i , for different sites, i, and AFM magnetic moments on Cu, m, depend on temperature, T , force constants and spin stiffness. Appropriate values of u i (T ) and m(T ) at T ≈ 100K are deduced from experiments and calculations  [ 16,17,18,19,20,21]. A striking result of these calculations (also made for HgBa 2 CuO 4 ) is that a gap (or a pseudogap) will open up in the DOS, very similar to what is found in the NFE models. The gap will appear near to E F for the undoped material if the wavelengths (Λ) are long, while for short phonons or spin waves the gap moves to lower energy [8]. In summary, Λ = 1/h, where Λ is in units of a 0 and h is the number of holes per Cu.
A second important result is that the spin waves will be stronger, and the pseudogaps will be deeper, if the waves coexist with phonon distortions of the correct wave length and phase [13]. Interactions between phonons and spin waves have also been suggested in order to ex-plain neutron scattering [22,23]. Phonons like the "halfbreathing" O-mode and modes with z-displacements of La and apical oxygens are most effective in this process of SPC. This mechanism offers direct explanations of phonon softening, q-dependent spin excitations and various isotope effects [10,14].
The rather few LMTO results permit to establish only a few values of V q for different q,i of the phonons and spin waves. A general trend of larger V q for long wave lengths seems clear, both for phonons and for spin waves. A procedure based on the partial character of the states above and below the gap in the NFE model confirms this trend for spin waves, and it permits to fill in some of the q-dependent points not calculated by LMTO. In all this gives some confidence in the q, i-variations of V q . However, one can expect more vivid variations of the spin part near q → 0, because LDA underestimates the transition towards AFM.
When these V q -values are fed into the 2D NFE model it gives an approximately linear variation of q as function of h, up to a saturation near q ∼ 0.125 for h larger than ∼ 0.13 [12]. Hence, the linear dependence for low h is in qualitative agreement with the result from 1D-LMTO, but the pace is different. This is probably because the gap opens both along x and y in the 2D NFE model, which makes the progression of the gap a bit slower than in 1D. The saturation is because the V q 's decrease for increasing q and it becomes impossible to open a gap at the low energy where E F should be for large h. However, if q x and q y are assumed to be different in the 2D NFE model it is possible to follow one gap towards lower energy. A second weaker gap moves to higher energy.

IV. SUPERCONDUCTIVITY
The emerging picture from these band results is that moderately strong spin fluctuations, which exist for many combinations of q and ω, will be enhanced through SPC to some particular phonon distortions. The selection depends on the atomic character of the phonon mode as well as on doping and wave length (i.e. on q). The latter is because the pseudogap appears at E F so that a maximum amount of kinetic energy can be gained from the SPC mode. This is for low T when the states below the gap is of one spin and well separated from the (unoccupied) state of the other spin above the gap. But for increasing T there will be mixed occupations of the states around E F through the Fermi-Dirac function (and through thermal disorder), which will decrease the spin density. This will decrease the spin polarization of the potential, which in turn will decrease the spin density even more, and at some temperature T * the support of the pseudo gap from the spin wave will collapse [10]. The fact that V q is largest at low doping favors larger T * when h → 0. A high DOS at E F is important for a high superconducting T C , and therefore is the pseudogap in competition with superconductivity in under-doped cuprates [12].
It was suggested that non-adiabatic electron-phonon coupling could be enhanced in the cuprates, because of the low Fermi velocity in the z-direction [24]. This velocity is comparable to the vibrational velocity and the electronic screening appeared to be insufficient. However, the screening can be made by other electrons moving within the planes. Moreover, phonon frequencies calculated within LDA without assumptions of incomplete screening, agree satisfactory with experimental frequencies [18,19,20]. Instead, excitations of virtual phonons coupled to spin waves can be important for the mechanism of superconductivity, but SPC makes the separation between pure electron-phonon coupling, λ ep , and λ caused by spin fluctuations, λ sf , less clear.
The important observation is that atomic phonon distortions will trigger enhancement of spin waves. If so, the possibility for larger λ ep is open, because instead for the common approach to ignore spin effects in electronphonon coupling, there are larger matrix elements when the spin polarized part of the potential is involved. This is readily imagined for a system which is nonmagnetic when phonon distortions are absent, but magnetic when the distortions are present. An estimation of pure λ ep and the coupling parameter for SPC, λ SP C , in LBCO has been presented earlier [22]. With a total N (E F ) ≈ 0.9(eV · cell · spin) −1 , and distortion amplitudes and potential shifts as in ref. [12] this leads to λ ′ s of the order 0.36 and 0.6 for pure phonons and SPC, respectively [22]. These values appear sufficiently large for a large T C in simple BCS-type formulations, but the precise relation for T C depends also on other parameters [25]. The third mechanism is that spin fluctuations work without coupling to phonon excitations. Still, phonon distortions might be present, and will in that case be important for enhancements of rapid (high energy) spin fluctuations. For instance, it was found that the local exchange enhancement (for AFM) on Cu sites are largest when the surrounding atoms (La or oxygens) have been pushed away from the Cu, and this leads to SPC between phonons and spin waves with equal q and ω [12]. But strong high-ω spin excitations should also be possible on rows of Cu with enhanced exchange, perpendicular to the phonon. Ab-initio calculations are not easy for such configurations because of the large size of the required unit cells, but one can use some of the existing results for first estimations of λ sf for rapid spin fluctuation (so rapid that the phonon distortions appear to be static compared to the spin fluctuation). The enhancements depend on the different type of distortions as was discussed above. In this case we calculate λ sf = N dV /dm 2 /(d 2 E/dm 2 ), where V is the potential, E is the total energy of the spin wave and m is the magnetic moment (per Cu) [26]. The difference in free energy between non-polarized and an induced (by magnetic field) AFM wave can be written E m = E 0 + κm 2 , so that d 2 E/dm 2 = 2κ, which permits to calculate λ sf = N · (∆V ) 2 /κ for "harmonic" spin fluctuations. With parameters from the band results this makes λ sf equal to 0.03 when no lattice distortions are present, and 0.12 and 0.11 when the lattice contain distortions, see Table 1.
The frequency has not been calculated explicitly, but a shortcut via exchange enhancements can be used for a simple estimation, as for FM fluctuations [26,27]. Thus, ω sf ≈ 1/(4N S), where S = ∆V /∆H is the local Stoner enhancement (on Cu) and ∆H is the external field ( 5 mRy in these calculations). This makes ω sf 100-200 meV depending on the type of lattice distortion. Fig. 5 shows the combined results for direct λ SP C and indirect λ sf as function of frequency. The general shape is similar as the electron-boson coupling function that has been extracted from recent optical spectra on Hg-and Bi-based copper oxides [28].
The amplitudes of λ SP C and λ sf increase further if both of the latter distortions (plane oxygens and La) are considered in the calculation. The matrix element increases to about 14.5 mRy. Despite this relative increase it is seen that pure spin fluctuations give not as large λ ′ s as for SPC. This is especially clear when all lattice modes are considered. However, these pure spin fluctuations can be important for superconductivity, since they appear at large ω. In contrast, SPC is essentially limited to below the highest phonon frequency, i.e. to 50-70 meV. The relatively small absolute values of λ sf are partly caused by the use of LDA. Improved DF schemes which could predict an AFM insulating state for the undoped systems, would lead to larger exchange enhancements, larger λ sf and lower ω sf in doped systems.

V. CORRELATION AND CORRECTIONS TO
LDA.
The question is if the real electron-electron correlation in the oxides is much larger than what is included in LDA. The latter includes exchange and correlation (XC) for an electron gas of varying density with the electron gas radius r s = 0.62ρ −1/3 as parameter. The density ρ contains one electron within r s . The largest values of r s (≈ 2a.u.) are found in the low density valence-electron region between the atoms, but they are still smaller than any typical atomic radius, R A in transition metals, their oxides and these cuprates. This gives a major argument against using a strong on-site correlation parameter to an atom, because the effect of an additional electron will readily be screened out, and is not noticed beyond r s . In low-density materials however, such as the alkali metals, r s ≈ R A and on-site correlation could be advocated. But LDA seems to describe alkali metals well, which provides an additional argument in favor of local approximations of DFT.
The fact is that LSDA does not produce an AFM insulating ground state of several oxides, so the problem is real and it may be related to the spin-polarized part of the potential. Small errors with important consequences of this type are known for LSDA. A well-known example is for Fe, where LSDA predicts nonmagnetic (NM) facecentered cubic (fcc) lattice as the ground state, whereas the generalized gradient approximation (GGA, [6]) correctly finds the ferro magnetic (FM) body-centered cubic (bcc) ground state [29]. The difference between LDA and GGA is not large concerning the bands, magnetic moments and other properties. The band properties of bcc Fe are quite good in calculations using both these DFT functionals. But the order of the total energies for FM-bcc and NM-fcc is reversed in the two cases.
An additional indication that the AFM insulating state of the cuprates is not far away in LSDA is the fact that the AFM moment depends on the number of k-points. A sufficient number of k-points is required for convergence, but early works noticed that the AFM ground state became stable if a coarse k-point mesh was used [30,31]. A similar conclusion for the appearance of weak ferromagnetism in highly doped LBCO can be understood from DOS effects and the Stoner criterion [32]. Thus, as for Fe, a small detail can have large consequences. It might be that the total energy difference between two very different ground states (fcc nonmagnetic and bcc-magnetic for Fe, and NM metallic and AFM insulating for cuprates) are very close, so that small errors in the computational details lead to the 'wrong' ground state.
Another indication of a nearby AFM state is that LMTO calculations with off-center linearization energies can stabilize AFM in HgBa 2 CuO 4 in calculations with sufficient number of k-points. The normal procedure in LMTO is to chose linearization energies at the center of the occupied sub bands, but if they are chosen at the bottom of the bands, it leads to a slight localization of the states. An AFM order needs less hybridization energy and AFM can be stabilized [33]. This is rather ad-hoc, since LMTO should not be made in such a way. A theoretical justification is need for a better DFT scheme.
A first workable non-local density functional is the GGA (generalized gradient correction [6]), where the gradient, or the first derivative of the charge density, dρ(r)/dr, brings out corrections to the LSDA. The gradients are mainly radial in atoms and solids because of the steep drop in charge density for increasing r close to an atom. Only valence electrons have large amplitudes in the interstitial region between the atoms, where the gradients are moderate. However, higher gradients such as the second derivative of the density can be large in the interstitial, a feature not captured by the GGA. Nevertheless, by using GGA in the LMTO calculations for LBCO there is some improvement, since the required critical magnetic field for having a gap is reduced by about 25 percent compared to LSDA.
Densities with quadratic dependence as function of r around a fixed point at r = 0 are used in a 2-particle model for calculations of correlation including corrections due to the second derivatives. A non-interacting density with quadratic gradients have the form The Schrödinger equation is solved for a density of electrons surrounding a fixed electron at r = 0, where it is taken into account that the effective mass is 1/2 [34,35]. The interacting potential can be written The exchange-correlation within the surrounding electron cloud is taken into account through LDA (µ xc ), and the external potential V ext can include relative kinetic energy variations. The two last terms in V (r) are not very important for the results for gradient variations. The strongest term is the unscreened Coulomb repulsion that diverges at r = 0. The second term, the centrifugal term is included only for exchange (with ℓ = 1, higher ℓ will have higher energy), when the Pauli principle requires that the density for equal spin must be zero at r = 0 [35]. The solution Ψ(E, r) of the Schrödinger equation for this potential is used to calculate the XC energy as the difference in Coulomb energy between interacting and non-inteacting densities; µ = (Ψ(E, r) 2 − ρ(r))/rd 3 r.
The energy E is determined from a required boundary condition of Ψ(E, r s ). The variation of the correlation potential on the second gradient is quite easy to understand. If the gradient is positive, so that there is a tendency for a 'hole' at r=0, it will be less effective for the Coulomb repulsion to make the correlation hole deeper at this point. Therefore, µ c will decrease for positive density gradients near the interstitial region. The gradient can change closer to the high density regions near the nuclei and make the correlation larger than for a constant density, but correlation becomes negligible in comparison to exchange when the density is high. The C-correction is parametrized through the electron gas parameter r s and Q = ρ(r s )/ρ(0), so that C = 1−5 √ Q/4+ (r s ) * Q/2 for Q > 1 and as the inverse of this parametrization if Q < 1. The expression is only applied for Q between -0.25 and 0.8, and C − 1 is used as a scaling factor of the correlation part from LDA. This correction makes the potential a little more repulsive in the interstitial region. The Cu-d states becomes more localized, which can promote AFM order, see Table II.
One fundamental theorem of DFT is that of "vrepresentability", that essentially says that there is a one-to-one correspondence between the density ρ(r) and potential V (r) [4]. The only possibility is that two equivalent charge densities can give an uninteresting shift in the potential, which can be absorbed by the shift of the energy spectrum. Imagine now that an external potential is applied to an electron gas of varying density ρ(r). The free electron wave functions will oscillate more if more kinetic energy is given to the electron gas, or less if kinetic energy is subtracted, even if ρ remains constant. The exchange part of the potential (X) will be modified because of the ability of the wave function to readjust itself around a second electron of the same spin (to create an exchange "hole"), as can be imagined from Slater exchange [7] or the two-particle model [35]. A similar change occurs for correlation (C), which is caused by the immediate Coulomb repulsion between all electron pairs. If ρ is constant in space there is only an uninteresting potential shift, but for varying charge densities there can be more subtle effects. The LDA is derived from the principle of density variations at the Fermi energy, and the relation between E F (which is a measure of the kinetic energy) and density is given by; This relation can be compared to E F (r) = E F − V (r) within real atoms, where E F is from the band calculation with potential V (r). Except for r very close to the nuclei there is not too large a difference between E (1) F (r) and E (2) F (r). In some regions the two values can be quite different. If so, it implies that the density is not in equilibrium with the kinetic energy as was used in the derivation of the LDA potential. Oxygen sites in oxides are often negatively charged, which suggests that the kinetic en- ergy should be smaller than what comes out from eq. 1. These effects are moderated by the fact that the charge transfer is a result of an attractive potential. It turns out that positive and negative deviations from eq. 1 exist in shells within all atoms.
Thus, band calculations with correction for the kinetic energy will be attempted where the local exchange potential is written where X(f e , r s ) is a scaling function of the normal Kohn-Sham potential µ KS , and f e = E (1) F (r s ). A derivation of the scaling function X(f e , r s ) can be done directly from the Slater functions by use of energy shifts in the arguments of the plane waves [7]. However, negative shifts leading to localized waves have to be avoided, and the resulting X-function seems too sensitive to small variations of f e . As an alternative we determine the scaling function from a two-particle model as in ref. [34,35], with a renormalization to make X(1, r s ) = 1. The result is displayed in fig. 6 for the appropriate range of r s and f e . The real value of f e for r ≤ 0.05R RW S increases and can be larger than 50 near the nuclei. This is extreme and f e is here cut off at 3. In preliminary calculations for undoped LBCO we do a rescaling of the exchange due to kinetic energy and of correlation due to second gradients. The comparison is made with standard LSDA for an AFM unit cell where a staggered magnetic field is applied to the Cu sites, and it is found that a field of ± 8.8 mRy is the limit for having a zero gap (E g ∼ 0.25mRy). The magnetic moment is then ±0.18µ B per Cu. The amplitude of the required magnetic field to obtain a zero gap is reduced when the calculations include the correction factors, see Table II.
The absolute values can depend on details of the band calculations [33], but the trend towards a stability of an AFM state is clear. Corrections of this type will be interesting for seeing enhanced spin fluctuations for long wave length spin waves in doped cuprates. Further enhancements of λ sf can also be expected [12]. Potential corrections must ultimately be tested for other types of materials in order to verify that some properties will not deteriorate. Here for the cuprates the corrections permit to proceed in the modeling of more properties from the band results.

VI. CONCLUSION
Results of band calculations for supercells with frozen phonons and spin waves suggest that a Fermi-liquid state can cause pseudo gaps and dynamic stripes. Together with a NFE model it is possible to describe the doping dependence of many normal state properties of the cuprates. The coupling between phonon distortions and spin fluctuations seems to be crucial for the mechanism of superconductivity, so that the spin-polarized part of λ is most enhanced by simultaneous excitations of phonons and spin waves. Two different mechanisms for superconductivity mediated by spin fluctuations are possible. The largest coupling parameter is when a phonon is excited together with the spin fluctuation. Lower couplings, λ sf at larger energies, are independent of the phonon excitation, but these spin fluctuations can nevertheless profit from possible phonon distortions of the lattice. Still, absolute numbers are too small when using LDA, as is also concluded from the absence of AFM stability in LDA calculations for undoped systems. However, it is argued that corrections due to higher order density gradients and kinetic energy are able to bring the band calculations closer to AFM. This is shown in calculations for undoped LBCO by yet very approximate corrections due to non locality and kinetic energy. Refinements of such corrections will be of interest for application to supercells, including doping, phonon distortions and spin waves, since it can be expected that realistic λ ′ s will be obtained.