Enhancement of superexchange due to synergetic breathing and hopping in corner-sharing cuprates

Cuprates with corner-sharing CuO4 plaquettes have received much attention owing to the discoveries of high-temperature superconductivity and exotic states where spin and charge or spin and orbital degrees of freedom are separated. In these systems spins are strongly coupled antiferromagnetically via superexchange mechanisms, with high nearest-neighbour coupling varying among different compounds. The electronic properties of cuprates are also known to be highly sensitive to the presence, distance and displacement of apical oxygens perpendicular to the CuO2 planes. Here we present ab initio quantum chemistry calculations of the nearest-neighbour superexchange antiferromagnetic (AF) coupling J of two cuprates, Sr2CuO3 and La2CuO4. The former lacks apical oxygens, whilst the latter contain two apical oxygens per CuO2 unit completing a distorted octahedral environment around each Cu atom. Good agreement is obtained with experimental estimates for both systems. Analysis of the correlated wavefunctions together with extended superexchange models shows that there is an important synergetic effect of the Coulomb interaction and the O–Cu hopping, namely a correlated breathing-enhanced hopping mechanism. This is a new ingredient in superexchange models. Suppression of this mechanism leads to drastic reduction in the AF coupling, indicating that it is of primary importance in generating the strong interactions. We also find that J increases substantially as the distance between Cu and apical O is increased. Cuprates exhibit exotic states because of the interplay between spin, charge and orbital degrees of freedom. Ab initio calculations now show that a mechanism called orbital expansion plays a key role in the magnetic properties of cuprates.

and m are not too large. Additionally in the SCF process, all orbitals are self-consistently optimized in the field of the CAS WF, to yield the variational minimum. The CAS WF is then augmented using a number of second-order techniques, including n-electron perturbation theory (NEVPT2) 18,19 , multi-reference linearized coupled cluster (MR-LCC2) 19,20 or multi-reference configuration interaction with single and double excitations (MR-CISD) 11 . These methods capture the remaining (weak) correlation involving electrons and orbitals outside of the active space. Using such approaches, it is possible, for example, to justify that the lowest electron removal state of La 2 CuO 4 is the Zhang-Rice singlet state 21 and study it in detail 22 . We use a variety of second-order methods to gauge their reliability. As the active space is enlarged, the corresponding second-order corrections diminish. The key question that arises is: What is the 'minimal' active space necessary to obtain a qualitatively correct reference WF sufficient to compute J reliably? We find that the necessary active space needs to be far larger than previously imagined, including relatively high-energy Cu 4d and O 3p orbitals. Exclusion of these from the active space leads to a dramatic underestimation of J.
We analyse the reason for the strong dependence of J on the active space and, in particular, the importance of 4d orbitals. As mentioned above, the superexchange mechanism depends on O-Cu hopping. The Coulomb energy cost U eff of this hopping is strongly reduced by an expansion of the Cu 3d orbitals, referred to as breathing 23 , when an electron hops into a Cu 3d orbital. This breathing effect at the same time increases t eff , the Cu-O effective hopping integral 24 . In a similar way, the O 2p orbital breathes as the O occupancy is changed. In the superexchange mechanism, J depends on both U eff and t eff (Supplementary Note 4) and the breathing effects therefore strongly influence J, as we shall shortly show.
The breathing effect involves a single 3d → 4d excitation, leading to an expansion of the charge density when an electron is added to the d shell. There are also important double 3d3d → 4d4d excitations, Enhancement of superexchange due to synergetic breathing and hopping in corner-sharing cuprates Nikolay A. Bogdanov 1 ✉ , Giovanni Li Manni 1 , Sandeep Sharma 1,2 , Olle Gunnarsson 1 and Ali Alavi 1,3 ✉ Cuprates with corner-sharing CuO 4 plaquettes have received much attention owing to the discoveries of high-temperature superconductivity and exotic states where spin and charge or spin and orbital degrees of freedom are separated. In these systems spins are strongly coupled antiferromagnetically via superexchange mechanisms, with high nearest-neighbour coupling varying among different compounds. The electronic properties of cuprates are also known to be highly sensitive to the presence, distance and displacement of apical oxygens perpendicular to the CuO 2 planes. Here we present ab initio quantum chemistry calculations of the nearest-neighbour superexchange antiferromagnetic (AF) coupling J of two cuprates, Sr 2 CuO 3 and La 2 CuO 4 . The former lacks apical oxygens, whilst the latter contain two apical oxygens per CuO 2 unit completing a distorted octahedral environment around each Cu atom. Good agreement is obtained with experimental estimates for both systems. Analysis of the correlated wavefunctions together with extended superexchange models shows that there is an important synergetic effect of the Coulomb interaction and the O-Cu hopping, namely a correlated breathing-enhanced hopping mechanism. This is a new ingredient in superexchange models. Suppression of this mechanism leads to drastic reduction in the AF coupling, indicating that it is of primary importance in generating the strong interactions. We also find that J increases substantially as the distance between Cu and apical O is increased.
which provide radial (in-out) correlations 11 . For a fixed number of d electrons, these correlations lead to a contraction of the charge density, at least if the basis has sufficient flexibility to satisfy the virial theorem. Correlation and breathing compete, making their simultaneous description complicated. Both effects lead to occupancy of 4d orbitals, but they are otherwise very different, and the ab initio calculations need to have the flexibility to capture both effects in a balanced way.

Ab initio calculations.
To study the electronic structure of cuprates, we employ the embedded cluster model. With this approach, accurate high-level calculations are performed for a small representative unit of the solid, while its environment is treated in a more approximate manner 25 . The details of the employed model are presented in Methods section.
We first perform CASSCF calculations with two singly occupied Cu 3d x 2 −y 2 orbitals in the active space CASSCF (2,2), similar to the one-band Hubbard model. Such minimal active-space calculations account for the unscreened Anderson superexchange mechanism (d 9 −d 9 and d 8 −d 10 configurations) and give a qualitatively correct AF J coupling. The value of the J obtained this way is, however, only ~20% of experimental data [26][27][28][29] (Table 1). As can be seen, the second-order corrections nearly double J, but are clearly insufficient. The uniform behaviour of the different dynamical correlation methods suggests that the extended Hilbert space (CASSCF(2,2) reference WF plus second-order perturbation) is inadequate to qualitatively describe the system. Enlarging the reference WF represents a natural remedy to this problem. The only exception is the difference-dedicated configuration interaction (DDCI) method, which gives J values very close to experiment on top of CASSCF(2,2) reference 8,30 . However, the DDCI is essentially a subspace of the MR-CISD, and significant differences of J calculated by these two methods imply that the description of electronic structure given by DDCI is far from being complete.
Because an electron hopping from the bridging O σ-bonding 2p y to the Cu 3d x 2 −y 2 orbital plays a crucial role in the superexchange (see, for example, ref. 4 ), this orbital is an obvious candidate to add into the active space. Such CASSCF(4,3) calculations roughly correspond to an unscreened multi-band Hubbard model. However, the obtained magnetic couplings turn out to be less than 1 meV higher compared with CASSCF(2,2). The reason is that, despite the inclusion of important ligand-hole determinants (d 9 -p 5 -d 10 and d 10 -p 4 -d 10 ), their energies are too high to be effective because the orbital optimization is primarily driven by the dominant d 9 -p 6 -d 9 configuration 5,7,30 . When we include the effect of further excited determinants at second-order level on top of the CASSCF(4,3) WF, J is still much smaller than the experimental value, indicating that important correlation effects are still missing.
To effectively lower the energy of d 9 -p 5 -d 10 , d 8 -p 6 -d 10 and d 10 -p 4 -d 10 determinants, one has to take into account orbital relaxation that comes along with them 5,7 . This can be done by adding a proper set of orbitals previously kept empty to the active space, namely Cu 4d and O 3p. Having additional d orbitals in the active space has been shown to be necessary to describe multiplet splittings for the late transition metals of the first row (see, for example, refs. [31][32][33]. Because of the variationality of the orbital optimization within the CASSCF procedure, the active orbitals are allowed to change, and a balanced choice of active space is required to ensure convergence. Such a balanced active space can be constructed with Cu 3d and 4d orbitals with e g character plus the bridging oxygen 2p y and 3p y orbitals 7 (Fig. 1, top three rows). Results for such CASSCF(8,10) calculations are shown in the third block of Table 1. The extension of the active space leads to a systematic differential effect with J increasing significantly at all levels of theory. Results close to experiment were reported using this active space together with a different formulation of the perturbation theory 7 , but our calculation results are about 80% of experimental values. To achieve a balanced description of all relevant effects, we consider all copper 3d and 4d orbitals, together with the bridging oxygen 2p and 3p orbitals, resulting in CASSCF (24,26). This active space yields a diagonalization problem in a space of ~10 14 Slater determinants, which we treat with DMRG and FCIQMC as approximate solvers 13,14 . With the additional many-body contribution from the Cu t 2g and π-bonding O orbitals taken into account in the large CASSCF, we find further stabilization of the singlet compared with the triplet. Second-order correction on top of the CASSCF(24,26) reference finally brings J close to the experimental values (see corresponding block in Table 1). The active orbitals are shown in Fig. 1; note that both 3d and 4d orbitals have significant amplitudes at the bridging oxygen p orbitals.
To verify that J values are converged with respect to the active space size, we performed even larger computations. We further correlated the peripheral O x 2 − y 2 -like 2p and 3p orbitals, the latter being strongly mixed with Cu 4s (Fig. 1). The last block in Table 1 shows the results of these CASSCF (28,30) calculations. The obtained J values are indeed close to the CASSCF(24,26) ones, with a slightly larger fraction being captured by CASSCF itself and less by the perturbation theory correction.
In a simple theory of superexchange 1,2 , a model of Cu 2 O is treated with one non-degenerate orbital on each atom. As discussed above, including only these orbitals in CASSCF (4,3) underestimates J by almost one order of magnitude. We now discuss why it is necessary to consider the large active space.  23 . Second, the hopping between the two sites is enhanced, as the Cu 3d orbital expands 24 . As in the CASSCF calculations, we use a fixed orthogonal basis set for all intermediate states. Therefore the breathing effect of a 3d orbital is described as a mixing of the 3d and 4d orbitals. The system can effectively expand or contract an effective 3d orbital, being a linear combination of a 3d and a 4d orbital, depending on their relative sign. To illustrate how this happens, we consider a Cu 2 dimer, including just one 3d and one 4d level on each atom, as indicated in Fig. 2. The levels have spin but no orbital degeneracy. We use the Hamiltonian ) . (1) The first index on c Liσ refers to the site, and the second labels the orbital. That is, i = 1(2) refers to a 3d (4d) orbital. The hopping between the Cu atoms is described by t ij . We include the direct on-site Coulomb integrals U 11 , U 12 and U 22 , describing 3d-3d, 3d-4d and 4d-4d interaction, respectively. K i refers to a Coulomb integral with three equal orbitals and the fourth different: These integrals are crucial for the breathing effect. If, for example, the 3d orbital on an atom is doubly occupied, the last term in equation (1) can excite a single electron from the 3d orbital ϕ 1 to the 4d orbital ϕ 2 . The effect of breathing is already evident in the atom (Supplementary Note 3). For simplicity, we here put t 12 = t 21 = − √ t 11 t 22 , U 12 = √ U 11 U 22 and K 1 /K 2 = √ U 11 /U 22 . We have used ε 2 − ε 1 = 24 eV, U 11 = 13 eV, U 22 = 10 eV, K 1 = −8 eV, t 11 = −0.5 eV and t 22 = −0.8 eV. Table 2 presents the singlet-triplet splitting obtained by solving the Hamiltonian in equation (1). It illustrates how inclusion of the integral K 1 strongly increases the splitting, because of breathing effects. To understand these results better, we consider a simpler calculation within only three configurations for the singlet state: where |vac⟩ is the vacuum state with no electrons. These configurations are shown schematically in Fig. 2. |1⟩ corresponds to the d 9 -p 6 -d 9 configuration mentioned above, while |2⟩ and |3⟩ resemble the d 8 -p 6 -d 10 configuration without and with 4d occupation, respectively. The Hamiltonian in equation (1) within the basis given by equation (3) reads Diagonalizing this matrix, we obtain the second column of Table 2. These results agree rather well with the full calculation for the model in equation (1), although the basis set in equation (3) is incomplete. The splitting is smaller because the higher-energy configurations have been neglected. We can now use Löwdin folding, focusing on the upper 2 × 2 corner of (H − z) −1 where ΔE = ε 2 − ε 1 + U 12 and we have introduced the approximation z ≈ 2ε 1 in some places. The matrix in equation (5) shows rather clearly that there is an interference between breathing and hopping from the 3d orbital on one site and the 4d orbital on the other site. The effective value of U has now been reduced, while the effective hopping has been increased, since K 1 < 0 and t 11 and t 12 have the same sign. For the triplet case, the basis state |2⟩ does not exist, and these renormalization effects are not present. The singlet-triplet splitting is then Distance from Cu site r (Å) This illustrates the importance of the renormalization of both U 11 and t 11 by the breathing effect.
The presented model only includes non-degenerate d orbitals. Including the full five-fold degeneracy increases the renormalization of U by approximately a factor of five. The model has only 3d-3d hopping, but it illustrates the breathing effects. The more realistic case of 3d-2p-3d hopping results in a more complex expression for J where renormalization of t and U cannot be disentangled (Supplementary Note 4).

Reduction of U.
The calculated bare on-site Coulomb integral between two 3d electrons is very large (~28 eV), leading to drastically suppressed charge fluctuations in the simplest model. For this reason, the CASSCF(2,2) and CASSCF(4,3) calculations give a very small J. However, by increasing the active space size, this energy cost can be strongly reduced. Crucial effects are the change of the effective radial extent of the 3d orbital (breathing) and rearrangements of the non-3d charge density as the number of 3d electrons varies (screening) 23 , which are captured in the CASSCF (24,26) calculation with second-order correction.
To disentangle these different effects, we performed a series of simpler, constrained calculations 34 . We put all hopping integrals from d (3d or 4d) basis functions on the Cu atoms equal to zero. We can then prescribe the total occupancy of d orbitals on each Cu atoms. We performed two calculations, one with the configuration d 9 -p 6 -d 9 and one with d 8 -p 6 -d 10 . In both cases, the system is allowed to relax fully, except that hopping to or from d orbitals is suppressed. We then obtain that the energy of the d 8 -p 6 -d 10 state is higher than the d 9 -p 6 -d 9 by 10-13 eV depending on initial conditions of constrained calculations (Supplementary Note 2). This means that bare U ≈ 28 eV has been reduced to U eff ≈ 10 eV. Experimental 35 and theoretical estimates [36][37][38][39] suggest that U eff is reduced even further (~8 eV). This may be owing to more long-ranged effects left out in our finite-size cluster calculation. Figure 3 shows charge differences due to breathing and screening for the d 8 -p 6 -d 10 calculation, discussed above. A calculation was first performed for the d 9 -p 6 -d 9 state, then a d electron was moved from one Cu atom to the other, keeping all orbitals unchanged. The corresponding densities at two copper sites are denoted ρ (d 8 ) and ρ (d 10 ) . This d 8 -p 6 -d 10 state is then allowed to relax self-consistently, giving the densities ρ d 8 and ρ d 10. The solid red curve in Fig. 3b shows the change in the charge density ρ d 10 − ρ (d 10 ) , illustrating how charge is moved from the inner part of Cu to the outer part (breathing). The dashed red curve shows the radial integral of the charge density difference, revealing that more charge is removed from the inner part than is added to the outer part. Since the number of d electrons is the same in the two calculations, non-d charge has been moved away from the Cu atom with the d 10 configurations as a response to the addition of one d electron (screening). Adding a d electron to a Cu atom thus only leads to an increase of the net charge by about half an electron, because of screening, which substantially reduces the energy cost.
Dependence of J on position of apical oxygens. As seen in Table 1, the magnetic coupling in Sr 2 CuO 3 is nearly two times larger than that in La 2 CuO 4 . In both cases, the computation of J is done using only two magnetic centres, so this difference should not be attributed to the dimensionality of the two materials. The other structural difference is the presence of apical oxygen ions in La 2 CuO 4 , which changes the local multiplet splittings, mainly the position of 3d z 2 levels 26,40-42 . The relative energy of the 3d z 2 orbital is believed to be connected to the shape of the Fermi surface and the value of the critical temperature in doped cuprates [43][44][45][46] .
There is experimental evidence that J also changes depending on the local geometry 47 . However, because different compounds have to be used experimentally, local distortions are accompanied by changes of Cu-O distances and type of adjacent metal ions. Therefore, it is instructive to investigate the dependence of J on the distance to apical oxygen ions in La 2 CuO 4 compound with an accurate computational method. We varied the apical O's positions within the cluster while keeping the electrostatic potential untouched and computed magnetic couplings using the procedure described above. The results of these calculations are presented in Fig. 4. It can be seen that, with increase of the distance to apical oxygen, the NN J grows. Moreover, the growth is faster when more electron correlation is taken into account. One obvious effect that leads to an increase of J is the lowering of the 4d z 2 orbital energy and corresponding enhancement of the orbital breathing. We observe 13% growth of the occupation of 4d z 2 orbitals upon 0.8 Å displacement of apical oxygens at the CASSCF (24,26) level.
The computational strategy presented and justified here was recently used to predict the superexchange strength in infinite-layer nickelate compounds 48 preceding the consistent experimental studies 49,50 .

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/ s41567-021-01439-1.