Disentangling the role of bond lengths and orbital symmetries in controlling $T_c$ in YBa$_2$Cu$_3$O$_7$

Optimally doped YBCO (YBa$_{2}$Cu$_{3}$O$_{7}$) has a high critical temperature, at 92 K. It is largely believed that Cooper pairs form in YBCO and other cuprates because of spin fluctuations, the issue and the detailed mechanism is far from settled. In the present work, we employ a state-of-the-art \emph{ab initio} ability to compute both the low and high energy spin fluctuations in optimally doped YBCO. We benchmark our results against recent inelastic neutron scattering and resonant inelastic X-ray scattering measurements. Further, we use strain as an external parameter to modulate the spin fluctuations and superconductivity. We disentangle the roles of Barium-apical Oxygen hybridization, the interlayer coupling and orbital symmetries by applying an idealized strain, and also a strain with a fully relaxed structure. We show that shortening the distance between Cu layers is conducive for enhanced Fermi surface nesting, that increases spin fluctuations and drives up $T_{c}$. However, when the structure is fully relaxed electrons flow to the d$_{z^2}$ orbital as a consequence of a shortened Ba-O bond which is detrimental for superconductivity

Even while more than three decades have passed since their discovery, the origin of superconductivity in cuprates remain highly debated. It is largely believed that at least for several variants of cuprates it is primarily spin fluctuations that drive superconductivity [1,2]. However, the issue is far from settled. It is difficult to resolve because the phase diagram is very dense: small excursions in parameter space drives the material from one phase to a new phase. Multiple low energy scales are present and sometimes intertwined. Theoretically, the challenge has been to build a material specific ability that incorporates all such interactions in the right proportions which can also predict ways to disentangle their roles. Optimally doped YBCO (YBa 2 Cu 3 O 7 ) is one of the higher T c superconductors in the cuprate family. In this work we present a high-fidelity ab initio theory that is designed to realize this objective, and we use it to explore how spin fluctuations and superconductivity can be modified through strain. In a recent ultrafast experiment [3] on YBa 2 Cu 3 O 7 large amplitude distortions on apical Oxygens were induced to modulate the inter layer coupling. Here we apply a uniaxial strain along the caxis. We perform two distinct excursions under strain; (a) an ideal strain where all atoms displace in proportion to their height along the c axis, and (b) allowing the internal coordinates to relax under the strain. We show that strain modifies superconductivity in both cases, but in different ways, thus highlighting how a detailed understanding of the mechanism is essential in order to predict and control unconventional superconductivity.
High resolution inelastic neutron scattering (INS) data for spin fluctuations exists [4,5], but the data is available only up to 60 meV. State-of-the-art recent resonant inelastic X-ray scattering (RIXS) picks up the signatures of bosonic fluctuations of different kinds [6], including those whose mechanisms are intertwined. Thus to de-cipher what constitutes the primary component of the observed RIXS spectra is a challenge for both theorists and experimentalists alike. RIXS data has been taken from the Cu-L 3 edge for YBa 2 Cu 3 O 7 along the (0,π) line [6] with excitations observed up to 300 meV. However, measurements could not be performed for the important momentum region around (π, π), which likely drives superconductivity.
In the present letter, we calculate the magnetic susceptibilities and superconducting instability for optimally doped YBCO using a new high fidelity, ab initio approach [7,8]. For the one-particle Green's function it combines the quasiparticle self consistent GW (QSGW ) approximation [9] with CTQMC solver [10] based dynamical mean field theory (DMFT). This framework [11,12] is extended by computing the local vertex from the two-particle Green's function by DMFT [13,14], which is combined with nonlocal bubble diagrams to construct a Bethe-Salpeter equation [15,16]. The latter is solved to yield the essential two-particle spin and charge susceptibilities χ d and χ m -physical observables which provide an important benchmark. Moreover they supply ingredients needed for the Eliashberg equation, which yields eigenvalues and eigenfunctions that describe instabilities to superconductivity. We will denote QSGW ++ as a shorthand for the four-tier QSGW +DMFT+BSE+Eliashberg theory. The numerical implementation is discussed in Pashov et al. [8] and codes are available on the open source electron structure suite Questaal [17]. Some details are also given in the supplemental material.
QSGW ++ has high fidelity because QSGW captures non-local dynamic correlation particularly well in the charge channel [8,18], but cannot adequately capture effects of spin fluctuations. DMFT does an excellent job at the latter, which are strong but mostly controlled by a arXiv:2012.04897v1 [cond-mat.str-el] 9 Dec 2020 local effective interaction given by U and J. In this letter, we have used U = 8eV and J = 0.7eV which is similar to what is generally used for cuprates [19]. That it can well describe superconductivity has now been established in several materials [15,16].
In YBCO, we explore the full potential of our ability by performing rigorous bench-marking of our computed magnetic susceptibilities (resolved in momentum and energy) against the low and high energy spectral data from INS and RIXS respectively. We show that we reproduce all intricate structures in momentum and energy spaces observed in INS and RIXS from our theory in a parameter free manner. That it is possible to reproduce the RIXS spectra from the spin susceptibility alone indicates that RIXS is measuring an excitation primarily magnetic in nature in this compound.
The superconducting glue the present theory can characterize originates from some combination of spin and charge susceptibility [13,15]. That we are able to reliably recover experimental neutron and RIXS data provides some confirmation that we have an adequate foundation to describe superconductivity of this kind. We can use the same method to probe how spin fluctuations and superconductivity are affected when the system is perturbed, in particular how T c evolves with tensile strain. In a prior work [15], this machinery was used to explain in detail how T c evolves with tensile strain in Sr 2 RuO 4 where it could be benchmarked against experiments. Fig. 1 benchmarks the dynamical structure factor S(q, ω) = (1 − e −βω ) −1 Im χ mag (q, ω) computed by QSGW ++ (Fig. 1b) against direct INS measurements of S (Fig. 1a) in the vicinity of the antiferromagnetic point Q AF = (π, π); and also against RIXS which measures any bosonic excitation, including the magnetic structure factor (Fig. 1c). For the former, we can compare only up to maximum value reported, 60 meV. At least in this range of energy and momentum QSGW ++ is in good qualitative and quantitative agreement. Since the RIXS measurement [21] did not sample the region close to Q AF , we compare to RIXS along the line Q connecting Γ and (π, 0). A second branch appears at high energy. For both low-energy and high-energy excitations the maxima of the peaks along Q line shown by the dashed line are in a very good agreement with the experimental data (blue dots). It confirms the magnetic nature of the excitation measured in this RIXS response, a subject of some controversy [22][23][24].
Next we use the Eliashberg theory derived from spin and charge susceptibilities to estimate T c , and investigate how it is affected by strain. We use QSGW ++ to simulate an uniaxial strain on the c−direction and study its effect on the superconducting order. The uniaxial strain is carried out by reducing the c−axis up to 8% with a concomitant expansion the plane. The volume change ∆V to experiment [25]. For each strain, we consider two scenarios : an ideal strain where all internal displacements are fixed to their projection along the c axis, and another case where atoms are relaxed to the zero-force condition. Forces are computed within density functional theory. We will denote these scenarios as 'SS' (for single shot) and 'SO' (for structure optimized). SO corresponds to the actual mechanical response of YBCO subject to an 33 tensile strain. For convenience of presentation we define a strain with the opposite sign of the usual definition: . For both scenarios, we compute the variation of the critical temperature T c by comparing the superconducting instability computed by solving the linearized Eliashberg equation (see Appendix). Comparing these two scenarios distinguishes two competing effects: on the one hand, the ideal strain changes the topology of the Fermi surface in a way that favors superconducting order. On the other hand, subsequently allowing the internal coordinates to fully relax empties the d z 2 orbital, which is  I. Interlayer Cu-Cu spacing (first column), copper to apical oxygen distance (second column), and vertical component of AO to Ba distance (third column). We report distances for the pristine material (first row), and under uniaxial strain ( z ). Parameters for ideal (SS) and fully relaxed (SO) structures are shown (see text).
unfavorable for the superconducting order.   2 shows these results in more detail. T c is greatly increased in the SS scenario, and the middle and bottom panels show how T c is correlative to magnetic susceptibility at (π, π). Both the spin fluctuation and the superconductivity instability show a similar trend which confirms that the spin fluctuation are an essential contributor of the superconductivity. We can also separate the contributions from stretching the c axis from the contributions by reductions in the basal plane by varying the Poisson ratio. In one case we used ν=0 which freezes the lattice vector in the basal plane (red symbols in Fig. 2). T c changes only marginally with ν, which indicate that the main effect are coming from the reduction of the c-axis. In another scenario, we expand a and b axes about 2% to match YBa 2 Cu 3 O 7 epitaxially on an STO substrate (green symbol).
On the (0, 0)−(π, π) line, the Cu d x 2 −y 2 in the two planes of YBa 2 Cu 3 O 7−δ couple weakly through the apical O, splitting these otherwise degenerate states into a bond-antibond pair. The Fermi surface connected with these orbitals splits into two sheets. fully relaxed strain (SO). SS and SO differ in the degree that d x 2 −y 2 and d z 2 are coupled, as discussed in the text. This is manifest by the colorbar, which shows the off-diagonal component |G(Q, ω= 0) z 2 ,x 2 −y 2 |, an indicator of the hybridization between d z 2 and d x 2 −y 2 on the Fermi surface.
As shown in Fig. 3, the Fermi surface is formed of three bands. The two curved lines correspond to the bonding and antibonding d x 2 −y 2 bands noted above. The interlayer hybridization is strongest at the two antinodal points where either Q x =π or Q y =π. As strain is applied, the interlayer distance decreases (Table I). This increases the interlayer hybridization, which further splits the bonding and antibonding d x 2 −y 2 Fermi surfaces. The antibonding surface becomes flatter and thus more square. Making the arc more square improves the nesting of momentum transfer Q=(π, π) for electrons living at the antinodal point. In d-wave superconductivity, it increases the attractive interaction [26,27], i.e works constructively for d-wave superconductivity. As shown in Fig. 4, not only the magnetic susceptibility increases but the nesting vector moves closer to Q AF where the magnetic susceptibility is maximum. These two effects cumulatively explain the large enhancement of T c observed in the SS scenario.
In SO scenario, we observe a similar splitting of the bonding and antibonding d x 2 −y 2 Fermi surfaces, but with an important difference. In SO case, the d x 2 −y 2 hybridize with d z 2 . This hybridization increases because relaxation reduces the Ba-AO vertical distance, e.g. by 25% when the c axis is reduced by 4% (Table I). As a consequence, the AO environment is changed which affects the Cu d z 2 orbital. In the unstrained case, d z 2 sits at −1.48 eV below the Fermi level, and it is marginally changed in the SS case (for z =4% it resides at −1.41 eV) while in the SO case, d z 2 is pushed closer to the Fermi Static magnetic susceptibility χ(q, ω=0) on the (0, 0)−(π, π) line, for pristine YBCO (black line), and YBCO subject to SS and SO kinds of strain. SS and SO induce opposite effects on both the peak position and amplitude χ(q, ω=0). In the SS case Tc increase both because χ(q, ω=0) increases and the peak shifts closer to QAF . level (−1.18 eV). The Fermi surface mainly composed of d x 2 −y 2 state becomes strongly hybridized with d z 2 . This is apparent from colorbar in Fig. 3. The two orbitals hybridize close the antinodal point which is known to be unfavorable to d -wave superconductivity [28][29][30]. Indeed, not only almost filled band, as d z 2 is unfavorable for T c but opposite spin nearest-neighbor coupling to form Cooper pair is less favorable when two orbitals are active at the Fermi level. When d z 2 orbital are not included in correlated subspace and in the Eliashberg equation, this detructive effect disappears.
To recap, in the idealized scenario compression of the planes increases the interlayer hybridization which enhances Fermi surface nesting and, hence, T c However, with a proper relaxation (SO), the Ba-AO bond length decreases dramatically to force out-of-plane contributions to the planar physics and changes the orbital components of the Fermi surface. This contribution works destructively for T c .
To conclude, we have first shown that QSGW ++ is able to predict magnetic fluctuations in YBa 2 Cu 3 O 7 with high fidelity. Using this technique to explore the parameter space when YBa 2 Cu 3 O 7 is subject to strain, we find that T c can be dramatically altered, but how it is altered depends on the details of the displacements. Subject to an ideal strain with no internal relaxation, T c is increased owing to enhanced interlayer hybridization, which changes the shape of the Fermi surface and makes nesting more favorable. However, in a more realistic scenario, the Fermi surface also suffers from a competing d z 2 hybridization which is detrimental to T c . ACKNOWLEDGMENTS CW acknowledges insightful and stimulating discussions with Antoine Georges. This work was supported by the Simons Many-Electron collaboration. CW was supported by grant EP/R02992X/1 from the UK Engineering and Physical Sciences Research Council (EP-SRC). F. J. are supported by the EPSRC Centre for Doctoral Training in Cross-Disciplinary Approaches to Non-Equilibrium Systems (CANES, EP/L015854/1). F.J is supported by the Simons Many-Electron Collaboration. For computational resources, we thank PRACE for awarding us access to SuperMUC at GCS@LRZ, Germany and Irene-Rome hosted by TGCC, France and Cambridge Tier-2 system operated by the University of Cambridge Research Computing Service (www.hpc.cam.ac.uk) funded by EPSRC Tier-2 capital Grant No. EP/P020259/1 and ARCHER UK National Supercomputing Service.

SUPPLEMENTARY MATERIAL
The calculation in the unstrained case was perfomed using crystal struture reported in [31]. Paramagnetic DMFT is combined with non magnetic QSGW via local projection on Cu 3d on the Cu augmentation spheres to form the correlated subspace. DMFT loop was performed using CTQMC impurity solver [10] and [32]. The two particle susceptibility needed in BSE for the magnetic susceptibility was computed using Exact Diagonalisation impurity solver with 6 bath sites on a mesh of 50 bosonic frequencies and 500 fermionic frequencies. A benchmark was done with CTQMC solver to check the accuracy of the hybridization fit. The Eliashberg equation was solved using an impurity susceptibility computed with CTQMC impurity solver [33].