Competition between d-wave superconductivity and magnetism in uniaxially strained Sr2RuO4

The pairing symmetry of Sr$_2$RuO$_4$ is a long-standing fundamental question in the physics of superconducting materials with strong electronic correlations. We use the functional renormalization group to investigate the behavior of superconductivity under uniaxial strain in a two-dimensional realistic model of Sr$_2$RuO$_4$ obtained with density functional theory and incorporating the effect of spin-orbit coupling. We find a dominant $d_{x^2-y^2}$ superconductor mostly hosted by the $d_{xy}$-orbital, with no other closely competing superconducting state. Within this framework we reproduce the experimentally observed enhancement of the critical temperature under strain and propose a simple mechanism driven by the density of states to explain our findings. We also investigate the competition between superconductivity and spin-density wave ordering as a function of interaction strength. By comparing theory and experiment, we discuss constraints on a possible degenerate partner of the $d_{x^2-y^2}$ superconducting state.


INTRODUCTION
Almost 30 years after the discovery of superconductivity in Sr 2 RuO 4 (SRO) [1], the symmetry of its superconducting order parameter (SCOP) remains an open question.Initially, its similarities with 3 Helium made it a prime candidate for spin-triplet pairing [2], corroborated by various measurements [3][4][5][6][7][8][9][10][11].Along with observations of time-reversal symmetry breaking (TRSB) supporting a two-component order parameter [12,13], the superconducting (SC) state was believed for a long time to be a chiral p-wave spin-triplet.However, conflicting evidence presented in various studies remained to be explained [14,15].First, the presence of nodal excitations is unexpected in a chiral p-wave SC [16][17][18][19].Second, the low critical field H c2 exhibited by SRO is typical for Pauli-limited spin-singlet SC [20] and the transition into the normal state upon applying a magnetic field appears to be first-order [21,22], with indications of a Fulde-Ferrell-Larkin-Ovchinnikov state for a certain parameter range, strongly pointing to a singlet SCOP [23].Third, none of the topologically protected edge states predicted in chiral p-wave states were observed in experiments [24][25][26][27], although this could be explained by a high Chern number [28].
In recent years, the chiral p-wave picture has basically been dismissed.First, the careful replications of key nuclear magnetic resonance experiments previously interpreted as supporting spin-triplet pairing have highlighted a heating effect and instead concluded that the SCOP corresponds to spin-singlet pairs [29][30][31].Second, applying uniaxial strain along the x principal crystallo-graphic axis was found to enhance the critical temperature (T c ) [32,33] and to lower the Fermi liquid coherence scale [34].This enhancement was shown to be maximal where the FS undergoes a Lifshitz transition, corresponding to a van Hove singularity (vHs) in the density of states (DOS), at a time-reversal invariant momentum point, inconsistent with odd-parity SCOPs like p-wave [35].Nowadays, a consensus appears to be crystallizing around the spin-singlet and even-parity natures of the SCOP, yet the debate is still ongoing.While ultrasounds experiments support the conclusion of a twocomponent order parameter [36,37] inferred by the observation of TRSB and the splitting between T c and the TRSB transition temperature [38,39], there are no twotemperature signature in bulk thermodynamical experiments such as specific heat and elastocalorimetry as well as scanning SQUID microscopy [40][41][42][43].As a result of this plethora of experimental evidence, SRO can be seen both as a critical playground for testing new theories with the goal of potentially unifying some of these contradicting observations and as a testbed to verify whether our interpretations of specific experiments are valid.Either way, it constitutes an ideal system to considerably advance our understanding of the mechanisms for unconventional superconductivity [44].
A general overview of possible ordering states in terms of their irreducible representations is given in Ref. 63, 73.In this paper, we investigate the leading superconducting instabilities of SRO using functional renormalization group (FRG) calculations [74], applied to a realistic model of the electronic structure derived from density functional theory (DFT) [75] that includes spin-orbit coupling (SOC).Note that previous studies of SRO using FRG were performed on tight-binding models, fitted to photoemission spectroscopy measurements [46,48,59,[76][77][78].In order to compare to experiments, we study the effect of uniaxial strain, tracking the evolution of T c as well as the type of ordering.We find a phase diagram with two different magnetic orders that compete with a single SCOP transforming like the B 1g irrep (often labelled as d x 2 −y 2 -wave).This competition is found to depend sensitively on the choice of interaction parameters.We show that a proper range of parameters lead to an increase of the superconducting T c in good agreement with experiments.

RESULTS
Electronic structure.-To describe the low energy electronic structure of SRO for the different strain values, we perform ab initio DFT calculations downfolded onto the t 2g orbitals of the ruthenium atoms using maximally localized Wannier functions as detailed in the methods section.There is an extensive literature using tightbinding Hamiltonians from ab-initio electronic structure calculations based on DFT for SRO, presenting overall very consistent results at the DFT level.Beyond DFT, it was shown that replacing the local SOC parameter of 100 meV obtained withing DFT with a value twice as large λ SOC ≃ 200 meV yielded a considerably better agreement with the experimentally observed Fermi surface [79], in consistent with the correlation-induced enhancement predicted theoretically [80,81].Indeed, this parameter adjustment leads to improved agreement with the FS both in the unstrained case and with the experimentally observed critical strain [79][80][81][82][83][84][85][86].A comparison to a selection of previously published tight-binding models is shown in App. A. Note that the addition of SOC breaks the SU (2) spin symmetry, but preserves an orbitally dependent SU (2) (so called pseudospin) symmetry [68].We keep the SOC fixed for all strain values [87].
In order to account for strong electronic correlations in this multi-orbital system, we use for most parts of this paper the O(3) symmetric Hubbard-Kanamori parametrization of the interaction Hamiltonian [88], which involves two key energy scales: the on-site Hubbard repulsion U and the Hund coupling J -see methods.As done routinely in FRG calculations [46,48,[76][77][78][89][90][91][92][93][94][95][96][97][98][99][100], we neglect the flow of the self-energy in our calcu-  Furthermore, we labelled the α, β and γ sheets on the FS for ϵxx = 0.In the unstrained case, the partial DOS for the dxz and dyz orbitals is identical, which is emphasized by the color mixing.
lations.Hence, the interaction parameters (U, J) should be considered as effective interactions with significance within our FRG framework rather than having a firstprinciple meaning.In this perspective, it is important to explore how the various instabilities are tuned by varying the interaction parameters.The guiding principles for constructing the model utilized follow two rules • We start from state of the art DFT with added spin-orbit coupling λ SOC such that we do have the right Fermi surface.
• The choice of our interaction parameters is restricted in such a manner that, from the Wannier model defined above, we recover the correct peak positions of the interacting spin-spin susceptibility as measured in experiments.
The FS and the density of states (DOS) obtained from this downfolded t 2g model are displayed on Fig. 1.The left column corresponds to the unstrained system (ϵ xx = 0) and the right columns to the optimally uniaxially strained system for which the Fermi level is at the vHs (ϵ xx = ϵ vHs xx ).Note that ϵ vHs xx does not include quasiparticle renormalization and therefore is not the same value as in experiments.The D 4h space group symmetry of the unstrained system is lowered down to D 2h by uniaxial strain and the B 1g irreducible representation of D 4h , of greatest relevance to our study, turns into the A g irreducible representation of D 2h .
Note the slightly unusual presentation of the FS in Fig. 1: this is due to the transformation from a tetragonal basis into a x-y plane which has to be done in this fashion to ensure periodicity of the downfolded model in the twodimensional primitive cell.Due to this, we have not a single but two k z values in the first primitive cell, i.e. the Z-point is located at the corner of the black square.All the results that follow are obtained for the quasi-2D model restricted to k z = 0 and k z = 2π/c, shown on Fig. 1.We verify that they are agreeing with the results from full 3D calculations in App.F.
The lowering of the symmetry under uniaxial strain lifts the degeneracy between the d xz and the d yz orbitals of the ruthenium atom, as seen in the DOS in Fig. 1.It also splits the d xy van-Hove singularity into two parts: one drifting away from the FS (x-direction) and one drifting towards the FS and crosses it at the Lifschitz transition (ϵ vHs xx ∼ 0.8% strain).On the FS shown in Fig. 1, we also highlight the dominant nesting vectors of the bare particle-hole susceptibility (χ 0 PH ), see App.B. First, q 1 = (2π(3a) −1 , π(3b) −1 ) (and all those related by symmetry) connects the α and β sheets of the FS.Second, q 2 = (π(2a) −1 , π(2b) −1 ) is connecting two van-Hove singularities and should become relevant at large interactions.Third, q 3 = (π(3a) −1 , π(3a) −1 ) also connects the α and β sheet of the FS.These vectors are consistent with the dominant spin fluctuations observed in neutron scattering experiments [101,102].Note that there is a family of nesting vectors connecting α and β sheet, all close to q 1 .With these insights from the non-interacting FS and DOS in mind, we proceed with the phase diagrams as a function of U and J for the unstrained and the ϵ xx = ϵ vHs xx cases.The results are presented on Fig. 2. The background color corresponds to the energy scale Λ c (expressed in meV) at which a divergence of the corresponding coupling is observed.The fastest divergent coupling corresponds to the dominant instability, which can either be superconductivity (in which case Λ c is expected to be proportional to the Berezinskii-Kosterlitz-Thouless [103-105] critical temperature T BKT ) or a spin-density wave (SDW) (in which case Λ c can be interpreted as the characteristic scale associated with the growth of the correlation length) [106] At the lowest U and J values, we find no divergence down to the lowest energy scale resolvable with our momentum resolution, and thus conclude that the system remains in the Fermi liquid (FL) state down to that scale.Apart from this unique point, we find three types of instabilities.Up to moderate U and low but finite J, we find a d x 2 −y 2 superconducting instability (corresponding to B 1g symmetry for the unstrained system, turning into A g for the strained one).Upon increasing U or J, we find that the dominant instability becomes a SDW with ordering vector q 1 .At even larger U and J, the system undergoes a high-temperature transition to another SDW phase characterized by the ordering vector q 2 .These ordering vectors are visible in both non-interacting and interacting susceptibilities.Since we do not incorporate the effect of the self-energy, we cannot observe the shift of the q 3 peak observed in Ref. [55].
Phase diagram and magnetic orderings.-The q 1 -SDW is driven by strong nesting between the α and β sheets.A corresponding peak in the spin-spin susceptibility has been well discussed both in the context of experimental observations [6,107] and theoretical discussions [54,67].It should be noted that this vector is connecting two different values of k z when backfolded in the three dimensional Brillouin zone.Its in plane analog, q 3 = (π(3a) −1 , π(3a) −1 ), was found to be subleading in earlier three-dimensional studies using the random phase approximation (RPA) [54].The q 3 peak is also found in DMFT calculations including vertex corrections [55,108].Here, we find the q 1 ordering to be the leading one, with the q 3 ordering also diverging but with smaller absolute magnitude.The increase of Λ c can be understood in terms of the Stoner criterion being fulfilled at a larger scale for larger U or J.At higher energy scales, the vHs are strongly smeared.This effect increases the importance of the q 2 ordering vector connecting two vHs points, leading to the emergence of the q 2 -SDW phase.
When applying uniaxial strain, the parameter range where we find a SDW is increased.This can be understood from the increase of the DOS at the Fermi level, which leads to a larger χ 0 PH and thereby a smaller interaction is required to fulfill the Stoner criterion.Beyond this effect, straining does not affect the structures of the phases and the q 2 -phase is still observable in the same parameter region, as the changes of the FS due to strain have counteracting effects: while in the y-direction the FS touches the vHs, it drifts further away from it in the x-direction.
Note that as we increase the strain beyond the Lifschitz transition, we do not find the SDW that is observed in experiments [38,41].The emergence of this phase has been understood as the removal of all curvature of the γ sheet between the upper/lower vHs and the X/X ′ points, which leads to strong nesting along this direction [109].We do not observe this phase at any investigated strain value.As was pointed out in Ref [110], strain suppresses quantum fluctuations preventing the ordered magnetic state, such that at high strain the phase transition emerges.Therefore, observing this transition in a diagrammatically unbiased approach is currently out of reach numerically, as it would require a full inclusion of the frequency dependence as well as multi-loop contributions [111].We hope to overcome this shortcoming in the future.
Superconductivity. -In the following, the superconducting phase is analyzed using a linearized gap equation on the FS.As shown in Fig. 3, we find a gap that transforms according to a B 1g for ϵ xx = 0 (A g for ϵ xx = ϵ vHs xx ) representation of the D 4h (D 2h ) point groups.In the band basis, this state has a dominant overlap with the d x 2 −y 2 harmonic and its main weight stems from the d xy orbital.Such a type of superconductor has been observed in several other studies [53, 54, 58, 59, 66-68, 109, 112].Although spin-orbit mixing distributes pairing over the different FS sheets, the dominant d xy orbital character that we find does lead to a larger pairing on the γ sheet.It should be noted however that experiments indicate that the gap function is sizeable on all three Fermi surface sheets.This potential difficulty could be resolved by employing an analysis based on the combination of FRG and mean-field theory [113].
The spectrum of the pair-pair susceptibility at Λ c contains the information of all possible subleading SCOPs.By analyzing this spectral distribution, we find a clear separation of the eigenvalue of the d x 2 −y 2 superconducting state by at least one order of magnitude from all eigenvalues of other SCOPs, for all parameters investigated.While this excludes any immediate degenerate state, no statement about the proximity of different symmetry states or individual critical temperatures can be drawn from FRG, because within this method the dominant instability is signalled by a divergent coupling and susceptibility [114].However, from the hierarchy standpoint, we still can extract tendencies towards different orderings as discussed in App.E. This hierarchy reveals that the p-wave pairing state [46,48,78] is always clearly subleading by a large margin and both extended s-wave and g-wave are consistently closer to the d-wave.We simulate the full 3D model at a selected point to check for consistency of our 2D simulations.The results of these 3D simulations are shown in App.F and agree with the 2D results.
The superconducting phase is generated by a spinfluctuation mechanism.The couplings U and J are cru-cial tuning knobs determining the onset of the phase and also control the transition to the neighbouring magnetic phase.When increasing U , the transition to a SDW is understood from the underlying ladder-type diagrams diverging as soon as U becomes larger than the critical value.Below that critical U , the still strong spin-fluctuations can drive a superconducting instability.However, increasing J has a more complex effect since it affects two different physical processes, which we discuss in terms of two distinct couplings, J ss and J dd in Eq. (1).J ss promotes spin-flip and pair-hopping processes, thus reducing the tendency to order magnetically while also increasing pair-correlations.J dd decreases inter-orbital density-density interactions, which reduces the inter-orbital repulsion between electrons on the same site.We unravel which of the two effects is most relevant for a) superconductivity and b) the magnetic transition.This is achieved by varying the two quantities independently, first in a simple RPA calculation and then in a full FRG calculation.
For the simple RPA calculation, we calculate χ 0 PH at Λ = 11.6 meV.We chose U = 0.3 eV to circumvent the Stoner instability and vary J ss and J dd between 0.0U and 0.3U independently.The dominant components of the bare susceptibility are presented in Fig. 8 of App.D. In general, we observe that varying J ss has barely any effect on χ RPA PH .J dd , on the other hand, increases the interorbital components by a significant amount.Therefore we expect J ss to have a weaker impact on the superconducting transition.Physically this is expected since J ss hampers the spin-fluctuations which are required to obtain an effective attraction required by the superconducting state.To support this claim and understand better the underlying interference mechanism, we developed a simple 2-band toy model in App.D.
In the full FRG simulation, we verify these conclusions, i.e. increasing J ss leads to a transition only at much larger values than the one for J dd .See Fig. 9 of App.D. Interestingly, J ss will generate a stronger admixture of higher order angular momentum superconductivity hosted by the d yz and d xz orbitals.These are however still sub-leading to the d x 2 −y 2 superconducting state.
Influence of strain.-Finally, we compare our results with experiments.We do so by examining the effect of strain from ϵ xx = 0.0% to ϵ xx = 1.3% on the leading instability of different (U, J) combinations.The general behavior of T c is consistent with earlier studies [76,109,115], while the predicted phases partially differ.The different critical scales Λ c can be interpreted as an estimate for T c of the instability.The results for all superconducting data points are summarized in Fig. 4.
For systems that start with a large initial critical scale at zero strain (Λ c (ϵ xx = 0)), no significant enhancement with respect to strain is found.The enhancement of T c is much larger when Λ(ϵ xx = 0) is smaller.This effect can be understood by looking at the DOS: large energy scales, or large temperatures, correspond to smeared out features in the DOS.Thus, the shift of the vHs due to ) group which requires some components to be purely imaginary due to the transformation behavior of the spin-orbitals under rotational symmetries [68].Since the gap function is only known up to a prefactor, we rescale it from −1 to 1, the sign and value is encoded in the colorbar.The first Brillouin-zone is marked by the grey square.Effect of strain ϵxx on the critical scale Λc for different values of U and J (upper left) and values of U and J for each line given in the upper right.DOS in the dxy orbital depending on ϵ (lower left) and theoretically predicted enhancement of Λc due to uniaxial strain as a function of its value for ϵxx = 0 (lower right).Each dotted line in the upper left plot corresponds to one U − J combination given in the upper right.There is a clear correlation between Λc(ϵxx = 0) and the ratio of increase in Tc which can be seen in both the upper left and lower right panels.The experimental data points in the lower right plot are extracted from [33].A proposed mechanism that explains this enhancement is detailed in the text.
strain is irrelevant since the vHs is not resolved, i.e. the DOS at the Fermi level does not change under strain.The lower Λ c , the sharper the vHs will become.Therefore, its shift enhances the DOS at the FS more strongly which in turn leads to a larger increase of T c .Thus, a lower Λ c (ϵ xx = 0) yields an enhancement of T c with ϵ xx which is both larger and taking place over a narrower range of strain.Once the vHs has crossed the Fermi level, T c is expected to go down again since we rapidly reduce the DOS at the Fermi level when straining further.Note that the nesting which generates the attractive interaction is not strongly asymmetric under strain, as can be seen in Fig. 1.Therefore, it cannot explain the asymmetry of Λ C observed in Fig. 4.However, the density of states is asymmetric, see Fig. 4, giving rise to the observed asymmetry in Λ C .
To compare our results to experiments, we evaluate Λ c (ϵ vHs xx )/Λ c (ϵ xx = 0) and plot it versus Λ c (ϵ xx = 0), hence measuring the increase of the critical scale depending on the initial one.We extract the corresponding experimental values from Ref. 33 by calculating the ratio of the maximal T c and the T c at ϵ xx = 0.These results are summarized in Fig. 4. We observe that the experiments indeed fit to the data predicted by FRG and we can extract a line of U and J combinations along which the experiment is reproduced.We find that the values on the line are around U = 1.1, 1.4 eV and J = 0.143 U, 0.1 U .Again, we emphasize that these should be considered as effective values valid within our FRG formalism.

DISCUSSION
In summary, we studied SRO starting from a firstprinciples description of its electronic structure and using a diagrammatically unbiased FRG approach.Using this framework, we investigated the influence of uniaxial strain as well as different contributions of Hund's coupling.We identified that the inter-orbital interaction reduction due to the density-density term J dd is the main driving force favoring superconducting order, which we found to be a pseudospin-singlet d x 2 −y 2 .Lastly, we showed that the experimental increase of T c as a function of strain can be recovered on a quantitative level from FRG simulations and from a comparison to these experiments we extracted effective values of the interaction parameters.
Our results highlight the dominance of a single d x 2 −y 2 SCOP that transforms like the B 1g representation (A g under uniaxial strain).We note that, while this SCOP agrees with many experimental measurements, it cannot explain the evidence for two-components and timereversal symmetry breaking.From the experimentally observed behavior of the time-reversal symmetry condition, we can infer that a partner of our found SCOP is required to remain invariant under moving the vHs through the Fermi level.This condition would be for example fulfilled by states with nodal lines along the x direction [60,61] or odd-frequency superconductors [68,69].
An interesting direction for future studies would be to investigate the effect of interaction terms consistent with D 4h symmetry but breaking full cubic symmetry.This could potentially influence the competition between different low energy orders [80].There are also many potential routes towards a more accurate investigation of the superconducting state.First, even though SRO is nearly perfectly layered, including the third spatial dimension increases the number of allowed ordering types [73].Secondly, taking into account the frequency-dependence of the gap function and self-energy would allow to gauge the relevance of the proposed odd-frequency state.Finally, an important advance would be to start the FRG flow from a correlated starting point i.e. from dressed quasiparticles and their effective interactions.This could for example be achieved by starting the FRG from a DMFT [116] description of the normal state.

Density functional
Theory.-We use density functional theory [117][118][119] and the Quantum ESPRESSO DFT package [120,121] with the PBE exchangecorrelation functional [122] to calculate the electronic structure.Cell parameters and internal coordinates of the crystal structure in the I4/mmm space group are relaxed in the conventional cell until all force components are smaller than 1 mRy/a 0 (a 0 : Bohr radius) and all components of the stress tensor are smaller than 0.5 kbar, yielding a relaxed in-plane (out-of-plane) lattice constant of a 0 = 3.878 Å(c = 12.900 Å).To calculate the strained structures, we fix one in-plane lattice constant of the conventional cell to the strained value, a new = (1 − s)a 0 , and relax the two orthogonal cell parameters as well as the internal coordinates as described above.After relaxation, we use the corresponding primitive unit cells containing one ruthenium atom each, i.e. three t 2g orbitals.We use scalar-relativistic ultrasoft pseudopotentials from the GBRV library [123], with the 4s and 4p (2s) semicore states for both strontium and ruthenium (for oxygen) atoms included in the valence.In the scalarrelativistic approximation, the spin-orbit coupling term is dropped [124].The energy cutoffs for the wave functions and charge density are set to 60 Ry and 720 Ry, respectively.We use a 12×12×12 Monkhorst-Pack k-point grid to sample the Brillouin zone, and a smearing of 0.01 Ry utilizing the Methfessel-Paxton scheme.To describe the low-energy physics we construct three rutheniumcentered t 2g -like maximally localized Wannier functions for each strained structure using Wannier90 [125][126][127].Spin-orbit coupling is included by first performing the DFT calculation without it and then adding a local SOC λ SOC = 200 meV to account for the correlation-induced enhancement over the DFT value.
Functional renormalization group.-In order to account for strong local electronic correlations in this multiorbital system, we consider the Hubbard-Kanamori interaction Hamiltonian [88] where U is the intra-orbital on-site Coulomb repulsion, while J dd (J ss ) is the density-density (spin-flip and pair hoping) part of the Hund's coupling.In the rotationally invariant formulation where O(3) symmetry is satisfied, J dd = J ss = J.Although usually equal, the distinct effect of J dd and J ss when they are treated independently is discussed in Appendix B. Our calculations are performed in the regime of U/W = 0.3 − 0.5, with W being the total bandwidth and U being the Hubbard interaction.We restrict ourselves to a situation with an on-site Hubbard-Kanamori interaction, as with this type of interaction, we reproduce the dominant spin-spin suszeptibilities peak position in comparision to the interacting spin-spin susceptibility measured in neutron scattering experiments [102], see App.B. Any interaction resulting in an analog peak structure should give the same physical results, as we argue in App.E. The strong electronic correlations emerging from the Hubbard-Kanamori interactions are incorporated to the non-interacting downfolded systems using the functional renormalization group (FRG) [74,128].FRG is technically an exact method to calculate the effective action functional of a given quantum action.It does so by introducing a scale-dependent cutoff (here we use a sharp energy cutoff) in the non-interacting propagator of the system.By taking derivatives with respect to this cut-off, one generates an infinite hierarchy of flow equations.In practice, this hierarchy must be truncated to become numerically tractable, making the method pertubatively motivated.
In this work, we employ the standard level-2 truncation, neglecting all three and more particle vertices.The validity of this approximation in the weak-tointermediate coupling regime can be motivated by a power-counting argument to prove the RG-irrelevance of higher order terms [74].Furthermore, we neglect the frequency dependence of the interaction, again motivated by the power counting argument, and the self-energy.This approach was applied to various systems including SRO [46,48,[76][77][78][89][90][91][92][93][94][95][96][97][98][99][100] and can be viewed as an diagrammatically unbiased extension of the random phase approximation.
In practice, we solve the flow equations from an energy scale much larger than the bandwidth and then integrate towards lower energies until we hit a divergence in one of the three diagrammatic channels labelled the particle-particle (PP), particle-hole (PH) and crossed particle-hole (PH) channels.A divergence is associated to a phase transition as the corresponding susceptibility also diverges.Information of the ordering type can be extracted from the susceptibilities as well as linearized gap equations [129].
We employ the truncated unity approximation which allows us to reduce the memory required computationally [130][131][132].
For the FRG simulations, we use the TU 2 FRG code [132].For convergence, we include all form-factors up to a distance of 8.2 Å, which amounts to a total number of 75 basis functions per orbital in the unit cell, i.e. all lattice harmonics up to the fifth (sixth) are included (depending on whether the constant is counted as harmonic or not).We checked for convergence by increasing the number of form-factors included near the phase transition between the magnetic and the superconducting phase for a few data points in the phase diagram.The simulations are performed on a 36 × 36 momentum-mesh in the x − y plane for the vertex function.The loop integration is performed using a FFT approach and an additional refinement of 45 × 45 is employed to achieve higher energy resolution.The results of the integration do not differ upon changing the resolution of the loop integration.
By using an enhanced value of the SOC, the effects of local interactions on SOC are already included on the single-particle level of the calculations.We do not suffer from double counting at that level since we neglect the flow of the self-energy.As a consequence however, our calculation does not take into account the renormalized effective mass of quasi-particles.
Acknowledgments.-JBP thanks Prof. E. Pavarini, Friedrich Krien, Lennart Klebl and Jacob Beyer for helpful discussion.The authors gratefully acknowledge the computing time granted through JARA on the supercomputer JURECA [133]  clearly overestimates the correlations at the X-point and on the connection line between X and M , similarly to what is obtained using the random phase approximation [67,68].Beyond that, FRG reproduces roughly the same shape and structure of the susceptibility as DMFT, but cannot reproduce its shifting of the peaks [55,108].As discussed in App.E, correcting for this overestimation in the effective interaction might influence the critical energy scales, but is not expected to drastically alter the hierarchy of the different order parameters.Most critically, we show analytically that the leading SDW peak observed in experiments [102] will also lead to an attractive interaction in the singlet channel, see Appendix E.
where l is an orbital and i, j are sites.We first examine the magnetic channel susceptibility which is given in the two-particle basis |l 1 ⟩ ⊗ |l 2 ⟩.The components of the bare particle-hole susceptibility are given by where G is the one-particle Green's function and K, K ′ and Q are four-momenta [68,135].From this expression, we find that the bare particle-hole susceptibility is diagonal in orbital space.Presuming the two orbitals to be degenerate, we write In RPA, the irreducible vertex in the particle-hole channel is the anti-symmetric static and local Coulomb tensor Γ PH = Λ PH .We spin-diagonalize it and take the magnetic channel such that Λ m PH ≡ Λ ↑↑↑↑ PH − Λ ↑↑↓↓ PH .It describes the irreducible set of particle-hole interactions that generates spin-fluctuations.As shown in the Appendix C of Ref. 68, it can be found, up to a minus sign due to conventions, as The parameters are the same as those characterizing the Hamiltonian of Eq. ( 1), in particular with the pairhopping and spin-flip terms characterized by J ph and J sf respectively.Consequently, the eigenvalues of Λ m PH χ 0 PH which quantify the amount of spin-fluctuations are In other words, a pair on l 1 can constructively or destructively interfere with a pair on l 2 via J ph , which enhances the intra-orbital components of the susceptibility.
On the other hand, J sf enhances the inter-orbital components.On the other hand, J dd acts in the same way as spin flips, affecting only the inter-orbital components of the interacting susceptibility.However it has twice the magnitude.Which terms play what role also crucially depends on the signs of χ 1 and χ 2 ; The corresponding numerical experiments are visualized in Fig. 8.They quantify the effects of the two different couplings and help to gauge their importance for the phase diagram.extended s-wave or d x 2 −y 2 (or p x and p y ).It merely gives a numerical prefactor to be put into the gap equation.Secondly, that prefactor is just dependent on the value of C(r) for r on the nearest neighbors; i.e. as long as the Fourier transform of C is attractive on the nearest neighbor, we generate superconductivity whose symmetry is determined by the particle-particle loops fermionic argument only.C(r) = dq e iqr C(q) (E5) For now, lets assume that C(q) (again we focus on the d xy orbital) is consisting of δ-peaks at q 2 or q 1 /q 3 with unit weight, for simplicity we furthermore go back to a square lattice unit cell.This allows us to calculate the Fourier transformation C(r) analytically.We focus on C(r x = 1, r y = 0), since all other terms can be understood via symmetries.We find C(r x = 1, r y = 0; q 2 ) = −1 and C(r x = 1, r y = 0; q 1 ) = e i2π/3 , meaning that Γ P (q 2 ) = −1 and Γ P (q 1 ) = −0.5 on the nearest neighbor shell.Thus the overestimation of correlations at X amplifies tendency towards superconductivity, but not towards one specific type.In more general terms, any transfer momentum of the form q i = (q, q) will generate attraction in the even channel as long as rq i > π/2.Below this threshold, the Fourier transformation of C(q) will become positive and thus generate attractive interactions in the odd-channel.Thus also the predicted leading momentum transfer found in DMFT should generate attraction for a singlet state.
To make again a connection to the model at hand, we analyse the attraction induced by an spin-fluctuation vertex for the full 3D model from RPA-flow.We present the attraction values generated by the leading eigenvalues of the interaction at each q in the xy plane, see Fig. 10.

Appendix F: Comparison to 3D results
To ensure that the results from our quasi 2D simulation are applicable to the realistic material, we check for consistency of our results at a few selected points between 2D and 3D calculations.In full 3D calculations more order parameters are allowed to be non zero [63,70,71,73], therefore besides the convergence check of our calculations this is also a consistency check of that the quasi-2D model indeed contains the leading instabilities of the 3D model.For this we simulate a 20 3 grid for the vertex function with an additional refinement of 21 3 for the momentum integration in the loop.First, we estimate the most relevant form-factor contributions.For this, we calculate C o1,o3 (q) and apply what has been discussed in section E, revealing that the most attractive contributions are the x and y direction nearest bond neighbor bonds as in our 2D simulations and their higher harmonics, see Fig. 10.
Thus, we include the first and second harmonic of this contribution in our calculation -i.e.all neighbors within a distance of 7.8 Å.The q = 0 SCOP results for U = 1.5eV,J = 0.15U and ϵ xx = 0.0%, 0.8% are summarized in Fig. 11.
Our results presented, indicate that the order parameter is indeed purely 2D, thus our 2D simulations capture the main aspects of the real material.Further the critical scale changes slightly from 3.15 (3.69) meV in 3D to 2.48(2.80)meV in 2D with 0.0(0.8)%uniaxial strain applied.Two remarks have to be made: The nearly layered structure gives rise to sub-leading pair-density wave contributions with the same fermionic momentum dependence as presented in Fig. 11 and finite transfer momentum of the form q = (0, 0, q z ).These are however subleading to the q z = 0 SCOP.The existence of these competing pair-density waves can be understood from a simple 2D toy model -assume a two band system which is non-dispersive in one direction: this immediately implies that both the particle-hole and particle-particle loop have to be independent of k y .Thus the spin-spin vertex is of the form and the linearized gap equation also being independent of k y .This implies that in the simple picture, all gaps of the form ∆ o1,o2 (k x , q y ) have equal eigenvalue λ independent of q y explaining why we observe these pair-density waves in the nearly layered system.

FIG. 1 .
FIG.1.Partial density of states and Fermi surface for the three t2g orbitals of the unstrained (left) and the optimally strained (right) systems in the upper anf lowe panels respectively.The optimal stain of 0.8% corresponds to the system being closest to the Lifschitz transition.Here, a (b) is the lattice parameter in the x (y) direction, with a = b in the ϵxx = 0 case.The black dots indicates the position of the vHs of the dxy orbital.The three dominant spin-density wave ordering vectors q1, q2 and q3 are highlighted in black, gray and pink, respectively.The first Brillouin zone is marked by a black square.We mark the Γ and Z point by red crosses.Furthermore, we labelled the α, β and γ sheets on the FS for ϵxx = 0.In the unstrained case, the partial DOS for the dxz and dyz orbitals is identical, which is emphasized by the color mixing.

FIG. 2 .
FIG. 2. Phase diagrams in the U -J parameter space for the unstrained (a) and the optimally strained (b) systems.The background color indicates the critical scale Λc, proportional to the ordering energy scale, with the corresponding color bar on the right.The phases are either a B1g superconductor (Ag under uniaxial strain), two different spin-density waves (SDWs) or a Fermi-liquid (FL).The q-vectors associated to the SDWs are shown in Fig. 1.

FIG. 3 .
FIG. 3. Example of superconducting order parameters found in Fig. 2, calculated from a linearized gap equation in orbital space.They correspond to U = 1.1 eV and J = 0.22 U .The 3 × 3 matrices represent the spin-orbital states of the paired electrons in the inter-pseudospin channel since the intra-pseudospin terms are vanishing.Each panel shows the momentum structure.We present in panel (a) and (b) the unstrained case with ϵxx = 0 and the optimally strained case with ϵ vHs xx respectively.The SCOP transforms like the one-component B1g (Ag) irreducible representation of the D 4h (D 2h) group which requires some components to be purely imaginary due to the transformation behavior of the spin-orbitals under rotational symmetries[68].Since the gap function is only known up to a prefactor, we rescale it from −1 to 1, the sign and value is encoded in the colorbar.The first Brillouin-zone is marked by the grey square.
FIG. 4.Effect of strain ϵxx on the critical scale Λc for different values of U and J (upper left) and values of U and J for each line given in the upper right.DOS in the dxy orbital depending on ϵ (lower left) and theoretically predicted enhancement of Λc due to uniaxial strain as a function of its value for ϵxx = 0 (lower right).Each dotted line in the upper left plot corresponds to one U − J combination given in the upper right.There is a clear correlation between Λc(ϵxx = 0) and the ratio of increase in Tc which can be seen in both the upper left and lower right panels.The experimental data points in the lower right plot are extracted from[33].A proposed mechanism that explains this enhancement is detailed in the text.

FIG. 10 .
FIG.10.Most attractive values of the projection of a spin-spin vertex to the pairing channel in the x-y plane.For this, we calculate an RPA-like flow for the Cchannel up to a scale of Λ = 1e −3 .We plot here the real part of Γ P (b, b) calculated from the dominant eigenvalue of the C channel at each k-point.To incorporate the two leading contributions, we require a from-factor cutoff of 7.8 Å in the full FRG calculation, visually indicated by a circle with radius 8.2 Å to emphasize that the largest attraction value is included.

FIG. 12 .
FIG.11.2D cuts through the leading eigenvector of the pairing vertex from TUFRG in kxky (upper panels) and kxkz (lower panels) in spin-orbital-momentum representation for the unstrained (left) and 0.8% strained (right) case.We find identical results to our 2D calculation.
at Forschungszentrum Jülich.JBP and DMK are supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under RTG 1995, within the Priority Program SPP 2244 "2DMP" -443273985 and under Germany's Excellence Strategy -Cluster of Excellence Matter and Light for Quantum Computing (ML4Q) EXC 2004/1 -390534769.The Flatiron Institute is a division of the Simons Foundation.Competing Interest.-The Authors declare no Competing Financial or Non-Financial Interests.