Superconducting order of $\mathrm{Sr}_2\mathrm{RuO}_4$ from a three-dimensional microscopic model

We compute and compare even- and odd-parity superconducting order parameters of strontium ruthenate ($\mathrm{Sr}_2\mathrm{RuO}_4$) in the limit of weak interactions, resulting from a fully microscopic three-dimensional model including spin-orbit coupling. We find that odd-parity helical and even-parity $d$-wave order are favored for smaller and larger values of the Hund's coupling parameter $J$, respectively. Both orders are found compatible with specific heat data and the recently-reported nuclear magnetic resonance (NMR) Knight shift drop [A. Pustogow et al. Nature 574, 72 (2019)]. The chiral $p$-wave order, numerically very competitive with helical order, sharply conflicts with the NMR experiment.


I. INTRODUCTION
Superconductivity was discovered in the layered perovskite strontium ruthenate, Sr 2 RuO 4 (SRO), about 25 years ago 1 . Muon spin relaxation and Kerr effect experiments indicated time-reversal symmetry breaking (TRSB) in the superconducting phase 2, 3 . The accompanied absence of a drop in the spin susceptibility 4 pointed toward a chiral p-wave order parameter [5][6][7] , which would make SRO an electronic analog of the A-phase of 3 He 8, 9 . In addition to the general interest in instances of unconventional superconductivity, intrinsic chiral p-wave superconductors are of particular importance owing to the possibility of enabling topological quantum computation with non-Abelian anyons 10-12 . However, a series of key experiments conflict with the above interpretation. The linear temperature dependence of specific heat 13 at low temperature implies nodes or deep minima in the gap 14 . Recent thermal Hall conductivity measurements further suggest vertical (out-ofplane) line nodes 15,16 . Uniaxial strain experiments see no indications of a T c -cusp, as expected for chiral p-wave order [17][18][19] . A very recent in-plane field NMR experiment measured a significant spin susceptibility drop 20 , contradicting the original measurements, which in the absence of spin-orbit coupling (SOC) would exclude all models featuring vectorial order parameters (triplet) pointing out of the basal plane 21 . This has reignited a longstanding debate, possibly making the case of helical or even-parity order parameters plausible 19,[22][23][24] . We note, however, that strong SOC [25][26][27] in a multi-orbital system complicates the analysis of the magnetic susceptibility compared to the single orbital case 9 .
Most studies so far have used a two-dimensional (2D) model, taking advantage of the quasi-2D nature of the dispersion relation of the relevant bands. However, the small corrugation of the cylindrical Fermi surfaces is deceptive, and actually hides a non-trivial k z dependence of the orbital content of the bands due to SOC, whose effect was shown to be three-dimensional (3D) 25 . A full 3D calculation is therefore warranted in order to study superconductivity in SRO 15,23,25,[28][29][30][31] , and especially to study the effect of SOC on the recent Knight shift experiments. Further, a 3D calculation is also required to study the possibility of horizontal line nodes, which have been proposed as a way to reconcile a nodal superconducting gap with TRSB 32 .
We propose an effective three-band three-dimensional model with on-site interaction, and we calculate the superconducting order parameter in the weak-coupling limit. As a function of the ratio of the Hund's coupling J to the Hubbard interaction strength U we find a transition at J/U ≈ 0.15 from an odd-parity helical phase with accidental (near-)nodes to an even-parity phase with symmetry-imposed vertical line nodes, both of which are compatible with several key experiments, but are not compatible with the observation of TRSB. This will be commented on in the Conclusions.
The paper is organized as follows. In Section II we describe the effective three-dimensional tight-binding model considered in this work. In Section III we outline the weak-coupling approach used to numerically determine the superconducting order parameter. The main results are presented in Section IV, where we provide detailed calculations of -and compare results with -two key experimental probes: the specific heat (Section IV B) and the NMR Knight shift (Section IV C). We conclude and discuss future prospects in Section V.

II. EFFECTIVE THREE-DIMENSIONAL MODEL
In SRO, three bands cross the Fermi energy and form quasi-two-dimensional Fermi surfaces, commonly denoted α, β, and γ 27,33-35 . Sheets α and β are formed mostly by the 4d xz and 4d yz ruthenium (Ru) orbitals, whereas the γ sheet stems mostly from the 4d xy Ru orbital. We construct a tight-binding model for the three active bands, based on the three Ru t 2g orbitals: (2) and ψ s (k) = [c As (k), c Bs (k), c Cs (k)] T . We here used the shorthand notation A = xz, B = yz, C = xy,s = −s, and s denotes spin (s = +1 is understood as ↑ and s = −1 means ↓). The annihilation operator for an electron with wavevector k and spin s on Ru orbital 4d a is denoted by c as (k). The matrix elements ε ab (k) account for intraand inter-orbital hopping, both in-and out-of-plane, and η sets the SOC amplitude. A priori we retain terms up to three sites apart in-plane and leading order terms, including inter-orbital terms, out-of-plane: To relate the above terms to Eq.
All terms in this model respect the crystal symmetries; they preserve inversion and time-reversal symmetry 23 . With the above conventions the first Brillouin zone is here defined as BZ = [−π, π] 2 × [−2π, 2π].

A. Band structure fit with Monte Carlo sampling
Within the 19-dimensional parameter space we seek the set {t} that globally minimizes the quantity where µ = α, β, γ is the band index, ξ µ (u µ C ) is the band energy (orbital content C = xy, as given by the eigenvector components of H s=+1 ) of the model in Eq. (2). Similarly,ξ µ (ũ µ C ) is the band energy (orbital content) of the model of Ref. 26, which is based on a tight-binding fit from spin-resolved angle-resolved photoemission spectroscopy (ARPES) data. Finally, w k (w k ) are chosen energy (orbital) weights. For the q µ 's we choose the three in-plane directions θ = 0, π/6, π/4 for the three k z cuts 0, π, 2π.
The hopping amplitudes were obtained by fitting the dispersion and orbital content of the 17-band model of Ref. 26   Specifically, we draw a set {t} for each MC cycle and accept it if it makes D({t}) smaller than the previously drawn set. Otherwise, it is retained as the new optimal set with probability exp −D({t}) 1/2 /T , where T is an artificially introduced 'temperature' that we gradually lower. For the momentum path we fit the band structure at the fitting points marked with crosses in Fig. 1. The points (Z, Γ, M , X, A, R) are weighted four times as much as the majority of the points, and points close to the Fermi energy are weighted four times as much as the remaining points. The above-mentioned points are defined in the primitive tetragonal unit cell as Γ = (000), Z = 00 1 2 , R = 1 2 0 1 2 , X = 1 2 00 , M = 1 2 1 2 0 , A = 1 2 1 2 1 2 . The orbital weightsw q were fixed to be comparatively smaller than the energy weights w k .
The optimal tight-binding parameters are summarized in Tables I and II. The expected uncertainties for the tight-binding parameters are O(10 meV). Since the data from the 17-band model in Ref. 26 incorporated terms with an accuracy threshold of 10 meV, our effective model can not hope to accurately describe terms of lower energy than this threshold. Thus, tight-binding parameters smaller than this threshold, having a negligible effect on e.g. the band structure, were neglected in the implementation of the method described in Section III.
The Fermi surfaces produced with the effective tightbinding model are shown in Fig 2; in Figs. 2 (a) -2 (c) we display the orbital content of the three Fermi surfaces, revealing that the β and γ sheets exhibit significant mixtures of the semi-2D orbtial (d xy ) and semi-1D Ru orbitals (d xz and d yz ), whereas the α sheet is dominated by the semi-1D orbitals.

III. WEAK-COUPLING THEORY
The projection of the Coulomb interaction on the onsite t 2g orbitals is given by where i is the lattice site, a is the orbital index, and n ias = c † ias c ias is the density operator. We assume that the phenomenological parameters satisfy U = U − 2J and J = J 37 . Since we will consider the weak-coupling limit in the following, this leaves the single parameter J/U characterizing the interaction.
We base our analysis on the weak-coupling scheme for repulsive Hubbard models, introduced and developed in Refs. 36, 38-50. We first diagonalize H 0 and index the eigenstates by a band index µ = α, β, γ and a pseudospin index, which we mostly keep implicit below. In this basis, the linearized gap equation reads where S ν is the Fermi surface of band ν, with |S ν | the corresponding Fermi surface area, and g the dimensionless matrix Here, Γ is the two-particle interaction vertex (see Ref. 36 for details) at leading (second) order, is the density of states, andv −1 µ = Sµ dk |Sµ| v µ (k) −1 . Approaching the weak-coupling limit U/t → 0 asymptotically, an eigenfunction ϕ of the Eq. 11 corresponding to a negative eigenvalue λ yields the superconducting gap below the critical temperature T c ∼ W e −1/|λ| , where W is the bare bandwidth.
Since we have chosen a pseudo-spin basis which is consistent with the tetragonal point group, each eigenvector ϕ belongs to one of its ten irreducible representations D 4h 46,51,52 . In Table III we list all ten irreducible representations of the tetragonal point group, with corresponding order parameter structures, using the standard decomposition of the order parameter into even-parity, d 0 (−k) = d 0 (k), and odd-parity, d(−k) = −d(k), components as 21 where s and s are pseudo-spin indices. Only the oddparity E u and even-parity E g representations are twodimensional and permit TRSB order without the need of fine tuning model parameters.
The Pauli principle assures that odd-parity (respectively, even parity) solutions correspond to pseudo-spin . Even-parity representations (subscript g) are described by a scalar (d0) order parameter, while odd-parity (subscript u) order parameters are described by a vector (d); see Eq. (14). In the second column one should associate fj with any function that transforms like sin kj under the point group operations, and f 2 j with a function that transforms like cos kj for j = x, y, z. Representations Eu and Eg are twodimensional and can favor TRSB combinations as indicated.

Representation
Order parameter A1g triplets (respectively, singlets). One should, however, keep in mind that a Zeeman field, as considered in Sec. IV C, couples to the physical spin, and not the pseudo-spin, which means that the behavior of the magnetic susceptibility cannot be deduced from the parity of the order parameter alone, and always requires a numerical calculation. The method outlined above is valid in the limit of weak interactions. One might reasonably object that this limit is not strictly satisfied in the case of SRO (and possibly other real materials) [53][54][55][56] . However, recent work (Ref. 57) reviewing a range of weak-and strong-coupling numerical methods applied to the twodimensional Hubbard model found that both the kdependency and the symmetry of the superconducting order exhibit surprisingly little variation between methods, suggesting a smooth transition from weak to strong coupling. Our aim is that the present study, using an advanced three-dimensional tight-binding model and a controlled approximation, provides helpful guidance for future studies.

IV. RESULTS AND DISCUSSION
For a repulsive interaction, there is no superconducting instability at first order 58 . At second order, which we limit ourselves to in the calculation presented in the main text, the vertex we calculate is described in depth in Refs. 36, 46, and 48. The leading eigenvalues in each irreducible representation are displayed as a function of J/U in Fig. 3. Whereas the even-parity orders show qualitatively different trends with J/U , the odd-parity states all show the same trend, and the splitting between them always remains small. The highest-T c state is the oddparity helical order A 1u for J/U < 0.15, and the evenparity d x 2 −y 2 order B 1g for J/U > 0. 15 The results presented here were obtained by discretization of the Fermi surfaces, using 3318 k-points for each of the three Fermi sheets shown in Figs. 2 (a) -2 (c) (20 points in the k z direction), thus solving the linearized gap Eq. (11) as a regular matrix eigenvalue problem. In the numerical implementation, the eigenvalues of Eq. (2) appear repeatedly in the two-particle vertex 36 of Eq. (12), and an effective diagonalization routine for the kinetic Hamiltonian constitutes the bulk of the numerical procedure 59 . Numerically, we computed sub-blocks of the gmatrix (Eq. (12)) simultaneously, using about 500 cores over the duration of a few weeks.

A. Gap structure
We focus first on the odd-parity helical phase realized for J/U < 0.15. The magnitude of the helical order parameter is displayed in Figs. 4 (a) -4 (c). Deep vertical minima, as first predicted in two dimensions 49,60 , are present on all three bands. Notably, the gap on the β band has min θ |∆ β (θ, k z ≈ π)| 0.02∆ 0 , where θ is the in-plane azimuthal angle and ∆ 0 the maximal gap, at locations in agreement with previous predictions 16,49,61 . The character of the minima also appear likely to agree with thermal Hall measurements 15,16 . For larger values of the Hund's coupling an even-parity B 1g phase is realized. The order parameter, shown in Figs. 4 (d) -4 (f), has symmetry-imposed vertical line nodes on all bands and additionally a suppressed gap in large regions of the β and γ sheets. We note that in two dimensions the B 1g phase does not appear in the weak-coupling limit for the physical range of J/U , but it does appear with the random phase approximation for finite values of U 22,49,60 . In Figs. 5 and 6 we display k z cuts of the helical and d-wave order parameters at J/U = 0.06 and J/U = 0.20, respectively. Note in particular the (near-)nodes on β in Fig. 5 (b).
In three dimensions, the possibility of an E g order parameter with a horizontal line node at k z = 0 emerges 32 . Interest in this state has been fueled by recent specific heat measurements 62 combined with the possibility of explaining both TRSB and a nodal gap. However, this sector turns out to be strongly disfavored in our weakcoupling limit: At (e.g.) J/U = 0.20, the best candidate has λ Eg /λ B1g ≈ 0.03 and thus does not come close to competing with the semi-two-dimensional order parameters found.

B. Specific heat
The eigenvector calculated at weak coupling, ϕ µ (k), is related to the order parameter via Eq. (13). We assume further that the order parameter factorizes as ∆(T )∆ µ (k), with max k ∆ µ (k) = 1. The temperature dependency of ∆(T ) is assumed to be that of a conventional BCS superconductor 63 . The generalized BCS relation 64 , fixing the overall size of the gap, when given the experimental value of T c ≈ 1.48 K, is where γ ≈ 0.5772 is Euler's constant, and where we introduced the average Here, θ is the in-plane polar angle, defined with vertex at (0, 0, kz/2) for β and γ, and vertex at (π, π, kz/2) for α.
multiband superconductor the specific heat per temper-ature per normal state value, γ n , can be expressed as 64 where E µ (k) = ξ 2 + ∆(T ) 2 |∆ µ (k)| 2 , and where the average here is evaluated as In Fig. 7 we show the specific heat as computed from Eq. 17 compared to experimental values from Ref. 13. For the helical order parameter (at J/U = 0.06) obtained at weak coupling, the minima are deep enough to practically behave as accidental vertical nodes, and to provide a decent agreement with specific heat data, as shown in Fig. 7 (a). The B 1g order parameter (at J/U = 0.20) also provides a fairly good match with specific heat data, as revealed by Fig. 7 (b). In contrast to previous studies 65 , there is no fine tuning of model parameters in these cal-culations (the specific heat dependency on J/U is weak). Crucially, and as pointed out in previous work, a gap comparable in size on all three bands seems necessary to reproduce the features of the experimental data 28,49,65 .

C. Magnetic susceptibility
Since both leading order parameters exhibit nodal behavior in practice, we turn to a different probe of the superconducting order: the spin susceptibility, as measured by the NMR Knight shift 20,66 .
We consider the normal state Hamiltonian with Zeeman terms when an external magnetic field is applied, where H 0 is here the normal state Hamiltonian in the absence of a magnetic field (Eq. (1)), H S (H L ) denotes coupling between the magnetic field H and the spin S (orbital L) degrees of freedom (cf. Ref. 67), Here, ε ab (k) are the single-particle orbital terms, σ a is the a'th Pauli matrix, µ B is the Bohr magneton (the vacuum permeability was fixed to µ 0 = 1), abc is the Levi-Civita symbol, and a, b, c run over orbitals xz, yz, xy (A, B, C, respectively). The matrix U (k) is defined to diagonalize the matrix H 0 (k) + H S + H L , where The (six) Fermi surfaces are defined by ξ µσ (k) = 0 for µ = α, β, γ and σ = +, −. In the absence of a magnetic field the two energies ξ µ± (k) become degenerate due to restored time-reversal symmetry. The matrix U (k) and the energies ξ µσ (k) contain all information needed to calculate the magnetization in the normal state, as defined below. The transformation between orbitals/spins (a, s) and bands/pseudo-spins (µ, σ) is given by the components of U (k): We define the magnetization (with the magnetic field and the measured response along direction i = x, y, z) as Using Eq. (25) and c † µσ1 (k)c νσ2 (k) = δ µν δ σ1σ2 f (ξ µσ1 (k)), where f is the Fermi function, the matrix elements of Eq. (27) are given by ). (28) In the superconducting phase we add superconducting terms at orbital level ∆ a1a2 s1s2 (k) = (d a1a2 The electron operators in orbital and spin basis are now expressed as linear combinations of their particle (u's) and hole (v's) constituents, In turn, this leads to the magnetization matrix elements Note that due to particle-hole symmetry the two terms of Eq. (32) contribute equally. To relate the orbital gaps of Eq. (30) to the order parameters obtained at weak coupling, we make use of the transformation where ∆ µ σ1σ2 (k) is the order parameter in band and pseudospin basis, and where u µσ as (k) are eigenvector components of H 0 (k) (i.e. in the absence of a magnetic field, crucially with the same gauge choice as in the weakcoupling calculation).
We use the spin susceptibility normalized by its normal state value as a proxy for the Knight shift in the superconducting state, where the subscripts refer to the normal state value, ∆ indicates a small increment in the external magnetic field in the linear response regime, and χ ii is a diagonal element of the magnetic susceptibility tensor. In the numerical evaluation of the momentum integral of Eq. (32) we associate for any given k the weak-coupling order parameter solution from the Fermi surface point closest to k, convoluted with a Gaussian damping factor set by the distance from the Fermi surface. In practice, the orbital coupling in Eq. (22) did not change the resulting Knight shift by any significant amount and was consequentially not included in the numerical results presented here. With SOC, the normal state spin susceptibility has a 'bulk' interband contribution that is not related to the Fermi surface, and that is therefore not affected by superconductivity. Further, as mentioned before, Cooper pairs only form well-defined singlets (respectively, triplets) in the pseudo-spin basis, but not in the physical spin basis. This leads to similar values of K x for both order parameters at zero temperature: K x (T = 0) = 0.45 for the d-wave order at J/U = 0.20 and K x (T = 0) = 0.59 for the helical order at J/U = 0.06. These numbers are in rough agreement with a recent NMR experiment 20 , which indicates a drop of around K x (T = 20 mK) ≈ 0.5. The only order parameter clearly seen to conflict with the experimental value is the chiral p-wave, which shows almost no drop (K x (T = 0) = 0.99). Note that, for textbook order parameters without SOC, there would have been a sharp contrast between d-wave (K x (T = 0) = 0) and helical (K x (T = 0) = 0.5) 9 .
Recently, the experimental results of Ref. 20 were reproduced by Ref. 66, where the temperature dependence of the Knight shift was also measured. Values of K x (T ) and K z (T ) for representative order parameters are presented in Fig. 8. With the the current experimental data it appears difficult to sharply distinguish a helical from a d-wave order parameter. However, the low-temperature values of K x arguably appear in best agreement with dwave order. Crucially, it would be desirable, although perhaps technically difficult, to have the NMR experiment repeated for out-of-plane fields as this would yield a crisp way to distinguish helical from d-wave order, as Fig. 8 (b) shows.

V. CONCLUSIONS
Both the d-wave and helical orders found in this calculation have vertical (near-)nodes, and seem compati-ble with specific heat data and recent Knight shift measurements 20 . However, despite fairly strong SOC, a chiral order parameter appears incompatible with the observed Knight shift drop. Further microscopic multiband Knight shift calculations would help in quantifying this, and an out-of-plane NMR experiment would help in definitely distinguishing the helical states from even-parity order parameters. While the OPs exhibit a substantial k z dependence on the β band, they remain overall fairly two-dimensional. We do not see any microscopic evidence of a favored E g gap with symmetry-imposed horizontal line nodes, at least in the weak-coupling limit.
An important outstanding aspect requiring further assessment, both theoretically and experimentally, is how to unify evidence of TRSB with either helical or evenparity order 68 . Should it turn out that TRSB is spurious, or unrelated to superconductivity, the scenario of a single-component d-wave order parameter would become a natural contender. Another possibility would be the formation of a two-component order parameter which couples different irreducible representations with accidentally close critical temperatures 22,30,69 . The neardegeneracy of the various odd-parity orders, found here and in previous work 49,60 , could potentially provide evidence for this scenario. We note that the four helical states in general are non-degenerate when J/U = 0, η = 0, and ε AB = 0 49 . However, the splitting of these states is a relatively small effect, as Fig. 3 shows. The (accidental) formation of a TRSB combination of helical states would be non-unitary.