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

TL;DR

This study uses a detailed three-dimensional microscopic model to analyze the superconducting order parameters of Sr2RuO4, revealing the favored pairing symmetries depending on Hund's coupling and their consistency with experimental data.

Contribution

It provides a comprehensive microscopic analysis of Sr2RuO4's superconducting states, highlighting the influence of Hund's coupling on pairing symmetry and reconciling theoretical orders with experimental observations.

Findings

Odd-parity helical order favored at smaller Hund's coupling.

Even-parity d-wave order favored at larger Hund's coupling.

Chiral p-wave order conflicts with NMR data.

Abstract

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.

Tables

Table 1. Table 1: Tight-binding parameters for ε 1 ​ D subscript 𝜀 1 D \varepsilon_{\mathrm{1D}} in Eq. ( 3 ) obtained with Monte Carlo sampling.

Parameter	$t_1$	$t_2$	$t_3$	$t_4$	$t_5$	$t_6$	$t_7$	$t_8$	${\mu }_{1D}$
Value [meV]	$257.8$	$27.8$	$-35.5$	$-22.4$	$-4.7$	$-2.4$	$3.2$	$54.5$	$286.9$

Table 2. Table 2: Tight-binding parameters for Eq. ( 4 ), ( 5 ), ( 6 ), and ( 7 ) obtained with Monte Carlo sampling.

Parameter	${\overline{t}}_1$	${\overline{t}}_2$	${\overline{t}}_3$	${\overline{t}}_4$	${\overline{t}}_5$	${\mu }_{2D}$	$\eta$	$t_{\mathrm{int},1}$	$t_{\mathrm{int},2}$	$t_{\mathrm{int},3}$
Value [meV]	$356.8$	$126.3$	$17.0$	$22.3$	$-4.1$	$351.9$	$59.2$	$2.0$	$-15.5$	$-5.4$

Table 3. Table 3: Irreducible representations of the tetragonal point group D 4 ​ h subscript 𝐷 4 ℎ D_{4h}   Sigrist and Ueda ( 1991 ) . Even-parity representations (subscript g 𝑔 g ) are described by a scalar ( d 0 subscript 𝑑 0 d_{0} ) order parameter, while odd-parity (subscript u 𝑢 u ) order parameters are described by a vector ( 𝒅 𝒅 \boldsymbol{d} ); see Eq. ( 14 ). In the second column one should associate f j subscript 𝑓 𝑗 f_{j} with any function that transforms like sin ⁡ ( k j ) subscript 𝑘 𝑗 \sin{k_{j}} under the point group operations, and f j 2 superscript subscript 𝑓 𝑗 2 f_{j}^{2} with a function that transforms like cos ⁡ ( k j ) subscript 𝑘 𝑗 \cos{k_{j}} for j = x , y , z 𝑗 𝑥 𝑦 𝑧 j=x,y,z . Representations E u subscript 𝐸 𝑢 E_{u} and E g subscript 𝐸 𝑔 E_{g} are two-dimensional and can favor TRSB combinations as indicated.

Representation	Order parameter
$A_{1g}$	$d_0\left(𝒌\right)=f_x^2+f_y^2$
$A_{2g}$	$d_0\left(𝒌\right)=f_x{f}_y\left(f_x^2-f_y^2\right)$
$B_{1g}$	$d_0\left(𝒌\right)=f_x^2-f_y^2$
$B_{2g}$	$d_0\left(𝒌\right)=f_x{f}_y$
$E_g$	$d_0\left(𝒌\right)=f_z\left(f_y\pm i{f}_x\right)$
$A_{1u}$	$𝒅\left(𝒌\right)=f_x\hat{x}+f_y\hat{y}$
$A_{2u}$	$𝒅\left(𝒌\right)=f_y\hat{x}-f_x\hat{y}$
$B_{1u}$	$𝒅\left(𝒌\right)=f_x\hat{x}-f_y\hat{y}$
$B_{2u}$	$𝒅\left(𝒌\right)=f_y\hat{x}+f_x\hat{y}$
$E_u$	$𝒅\left(𝒌\right)=\left(f_x\pm i{f}_y\right)\hat{z}$

Equations

$$H_0=\sum _{\mathbf{k},s}{\mathbf{\psi }}_s^{†}(\mathbf{k})\text{pazocal}{H}_s(\mathbf{k}){\mathbf{\psi }}_s(\mathbf{k})$$

$$\text{pazocal}{H}_s(\mathbf{k})=\left(\begin{array}{ccc}{\epsilon }_{AA}(\mathbf{k})& {\epsilon }_{AB}(\mathbf{k})-is\eta & {\epsilon }_{AC}(\mathbf{k})+i\eta \\ {\epsilon }_{BA}(\mathbf{k})+is\eta & {\epsilon }_{BB}(\mathbf{k})& {\epsilon }_{BC}(\mathbf{k})-s\eta \\ {\epsilon }_{CA}(\mathbf{k})-i\eta & {\epsilon }_{CB}(\mathbf{k})-s\eta & {\epsilon }_{CC}(\mathbf{k})\end{array}\right),$$

$${\epsilon }_{1D}(k_{\parallel },k_{\perp },k_z)$$

$$-2t_6\cos (3k_{\parallel })-2t_7\cos (2k_{\perp })-2t_8\cos (k_{\parallel }/2)\cos (k_{\perp }/2)\cos (k_z/2)-{\mu }_{1D},$$

$${\epsilon }_{2D}(\mathbf{k})$$

$$-4{\bar{t}}_4\left[\cos (2k_x)\cos (k_y)+\cos (2k_y)\cos (k_x)\right]-2{\bar{t}}_5\cos (k_z/2)\cos (k_x/2)\cos (k_y/2)-{\mu }_{2D},$$

$${\epsilon }_{AB}(\mathbf{k})$$

$${\epsilon }_{AC}(\mathbf{k})$$

$${\epsilon }_{BC}(\mathbf{k})$$

$$\{t\}=\big{\{}\{t_{i}\}_{i=1}^{8},\{\bar{t}_{i}\}_{i=1}^{5},\{t_{\mathrm{int},i}\}_{i=1}^{3},\mu_{\mathrm{1D}},\mu_{\mathrm{2D}},\eta\big{\}},$$

$$D$$

$$+$$

$$H_I$$

$$+\frac{J}{2}\sum _{i,a\ne b,s,s^{\prime }}{c}_{ias}^{†}{c}_{ib{s}^{\prime }}^{†}{c}_{ia{s}^{\prime }}^{}{c}_{ibs}^{}$$

$$+\frac{J^{\prime }}{2}\sum _{i,a\ne b,s\ne s^{\prime }}{c}_{ias}^{†}{c}_{ia{s}^{\prime }}^{†}{c}_{ib{s}^{\prime }}^{}{c}_{ibs}^{},$$

$$\sum _{\nu }{\int }_{S_{\nu }}\frac{\mathrm{d}{\mathbf{k}}_{\nu }}{|S_{\nu }|}g({\mathbf{k}}_{\mu },{\mathbf{k}}_{\nu })\phi ({\mathbf{k}}_{\nu })=\lambda \phi ({\mathbf{k}}_{\mu }),$$

$$g({\mathbf{k}}_{\mu },{\mathbf{k}}_{\nu })=\sqrt{\frac{{\rho }_{\mu }{\bar{v}}_{\mu }}{v_{\mu }({\mathbf{k}}_{\mu })}}\Gamma ({\mathbf{k}}_{\mu },{\mathbf{k}}_{\nu })\sqrt{\frac{{\rho }_{\nu }{\bar{v}}_{\nu }}{v_{\nu }({\mathbf{k}}_{\nu })}}.$$

$$\Delta ({\mathbf{k}}_{\mu })\sim \sqrt{\frac{v_{\mu }({\mathbf{k}}_{\mu })}{{\bar{v}}_{\mu }{\rho }_{\mu }}}\phi ({\mathbf{k}}_{\mu })$$

$$\Delta_{ss^{\prime}}(\boldsymbol{k})=\big{[}\left(d_{0}(\boldsymbol{k})\mathds{1}+\boldsymbol{d}(\boldsymbol{k})\cdot\boldsymbol{\sigma}\right)i\sigma_{y}\big{]}_{ss^{\prime}},$$

$$\frac{\Delta (0)}{k_B{T}_c}=\pi \exp (-\gamma -⟨\text{leavevmode}⟨\log |\Delta (\mathbf{k})|{⟩\text{leavevmode}⟩}_{\mathrm{FS}}),$$

$$⟨\text{leavevmode}⟨\log |\Delta (\mathbf{k})|{⟩\text{leavevmode}⟩}_{\mathrm{FS}}=\frac{{\sum }_{\mu }{\int }_{S_{\mu }}\mathrm{d}\mathbf{k}\frac{|{\Delta }_{\mu }(\mathbf{k}){|}^2}{v_{\mu }(\mathbf{k})}\log |{\Delta }_{\mu }(\mathbf{k})|}{{\sum }_{\nu }{\int }_{S_{\nu }}\mathrm{d}\mathbf{k}\frac{|{\Delta }_{\nu }(\mathbf{k}){|}^2}{v_{\nu }(\mathbf{k})}}.$$

$$\frac{C(T)}{T\gamma_{n}}=\frac{3}{2\pi^{2}(k_{B}T)^{3}}\int_{0}^{\infty}\mathrm{d}\xi\hskip 2.84526pt\Big{\langle}\frac{\xi^{2}+\Delta(T)^{2}\lvert\Delta_{\mu}(\boldsymbol{k})\rvert^{2}-\frac{T}{2}\frac{\partial\Delta(T)^{2}}{\partial T}\lvert\Delta_{\mu}(\boldsymbol{k})\rvert^{2}}{\cosh^{2}(\frac{E_{\mu}(\boldsymbol{k})}{2k_{B}T})}\Big{\rangle}_{\mathrm{FS}},$$

$$⟨A{⟩}_{\mathrm{FS}}=\frac{1}{{\sum }_{\nu }{\rho }_{\nu }}\sum _{\mu }\frac{{\rho }_{\mu }}{|S_{\mu }|}{\int }_{S_{\mu }}\mathrm{d}\mathbf{k}A.$$

$$H=H_0+H_S+H_L,$$

$$H_0$$

$$H_S$$

$$=-{\mu }_B\sum _{\mathbf{k},s,s^{\prime }}\sum _{a,b}{H}_b{\sigma }_{s{s}^{\prime }}^b{c}_{as}^{†}(\mathbf{k})c_{a{s}^{\prime }}(\mathbf{k}),$$

$$H_L$$

$$=-{\mu }_{B}i\sum _{\mathbf{k},s,s^{\prime }}\sum _{a,b,c}{H}_c{ϵ}_{abc}{c}_{as}^{†}(\mathbf{k})c_{b{s}^{\prime }}(\mathbf{k}).$$

$$H$$

$$\mathbf{\psi }(\mathbf{k})$$

$$c_{as}(\mathbf{k})=\sum _{\mu ,\sigma }[u_{as}^{\mu \sigma }(\mathbf{k}){]}^{\ast }{c}_{\mu \sigma }(\mathbf{k}).$$

$$\text{pazocal}{M}^{ii}$$

$$M_a^{s_1{s}_2}$$

$$M_a^{s_1{s}_2}={\mu }_B\sum _{\mathbf{k},\mu ,\sigma }{u}_{a{s}_1}^{\mu \sigma }(\mathbf{k}){\left[u_{a{s}_2}^{\mu \sigma }(\mathbf{k})\right]}^{\ast }f({\xi }_{\mu \sigma }(\mathbf{k})).$$

$$H_{\Delta }=\sum _{\mathbf{k}}\sum _{a_1,a_2,s,s^{\prime }}{\Delta }_{s_1{s}_2}^{a_1{a}_2}(\mathbf{k})c_{a_{1}s}^{†}(\mathbf{k})c_{a_2{s}^{\prime }}^{†}(-\mathbf{k})+\text{h.c.},$$

$$\Delta_{s_{1}s_{2}}^{a_{1}a_{2}}(\boldsymbol{k})=\big{[}\left(d_{0}^{a_{1}a_{2}}(\boldsymbol{k})\mathds{1}+\boldsymbol{d}^{a_{1}a_{2}}(\boldsymbol{k})\cdot\boldsymbol{\sigma}\right)i\sigma_{y}\big{]}_{s_{1}s_{2}}.$$

$$c_{as}(\mathbf{k})=\sum _{\mu ,\sigma }\left(u_{as}^{\mu \sigma }(\mathbf{k})c_{\mu \sigma }(\mathbf{k})+v_{as}^{\mu \sigma }(-\mathbf{k})c_{\mu \sigma }^{†}(-\mathbf{k})\right).$$

Full text

I Introduction

Superconductivity was discovered in the layered perovskite strontium ruthenate, Sr2RuO4 (SRO), about $25$ years ago Maeno et al. (1994). Muon spin relaxation and Kerr effect experiments indicated time-reversal symmetry breaking (TRSB) in the superconducting phase Luke et al. (1998); Xia et al. (2006). The accompanied absence of a drop in the spin susceptibility Ishida et al. (1998) pointed toward a chiral $p$-wave order parameter Mackenzie and Maeno (2003); Maeno et al. (2012); Kallin and Berlinsky (2016), which would make SRO an electronic analog of the A-phase of ${}^{3}\mathrm{He}$ Leggett (1975); Rice and Sigrist (1995). 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 Read and Green (2000); Ivanov (2001); Nayak et al. (2008).

However, a series of key experiments conflict with the above interpretation. The linear temperature dependence of specific heat NishiZaki et al. (2000) at low temperature implies nodes or deep minima in the gap Deguchi et al. (2004). Recent thermal Hall conductivity measurements further suggest vertical (out-of-plane) line nodes Hassinger et al. (2017); Dodaro et al. (2018). Uniaxial strain experiments see no indications of a $T_{c}$-cusp, as expected for chiral $p$-wave order Hicks et al. (2014); Steppke et al. (2017); Li et al. (2019). A very recent in-plane field NMR experiment measured a significant spin susceptibility drop Pustogow et al. (2019), 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 Balian and Werthamer (1963). This has reignited a longstanding debate, possibly making the case of helical or even-parity order parameters plausible Rømer et al. (2019a); Ramires and Sigrist (2019); Li et al. (2019); Tanaka et al. (2009). We note, however, that strong SOC Haverkort et al. (2008); Veenstra et al. (2014); Tamai et al. (2019) in a multi-orbital system complicates the analysis of the magnetic susceptibility compared to the single orbital case Rice and Sigrist (1995).

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) Haverkort et al. (2008). A full 3D calculation is therefore warranted in order to study superconductivity in SRO Haverkort et al. (2008); Mackenzie et al. (2017); Hassinger et al. (2017); Røising et al. (2018); Ramires and Sigrist (2019); Huang et al. (2019); Gingras et al. (2018), 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 Žutić and Mazin (2005).

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\approx 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.2) and the NMR Knight shift (Section IV.3). 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 $\alpha$, $\beta$, and $\gamma$ Damascelli et al. (2000); Bergemann et al. (2000, 2003); Tamai et al. (2019). Sheets $\alpha$ and $\beta$ are formed mostly by the $4d_{xz}$ and $4d_{yz}$ ruthenium (Ru) orbitals, whereas the $\gamma$ 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:

[TABLE]

where

[TABLE]

and $\boldsymbol{\psi}_{s}(\boldsymbol{k})=[c_{As}(\boldsymbol{k}),\hskip 2.84526ptc_{Bs}(\boldsymbol{k}),\hskip 2.84526ptc_{C\bar{s}}(\boldsymbol{k})]^{T}$. We here used the shorthand notation $A=xz$, $B=yz$, $C=xy$, $\bar{s}=-s$, and $s$ denotes spin ($s=+1$ is understood as $\uparrow$ and $s=-1$ means $\downarrow$). The annihilation operator for an electron with wavevector $\mathbf{k}$ and spin $s$ on Ru orbital $4d_{a}$ is denoted by $c_{as}\left(\mathbf{k}\right)$. The matrix elements $\varepsilon_{ab}(\boldsymbol{k})$ account for intra- and inter-orbital hopping, both in- and out-of-plane, and $\eta$ 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:

[TABLE]

To relate the above terms to Eq. (2) we set $\varepsilon_{AA}(\boldsymbol{k})=\varepsilon_{\mathrm{1D}}(k_{x},k_{y},k_{z})$, $\varepsilon_{BB}(\boldsymbol{k})=\varepsilon_{\mathrm{1D}}(k_{y},k_{x},k_{z})$, and $\varepsilon_{CC}(\boldsymbol{k})=\varepsilon_{\mathrm{2D}}(k_{x},k_{y},k_{z})$. All terms in this model respect the crystal symmetries; they preserve inversion and time-reversal symmetry Ramires and Sigrist (2019). With the above conventions the first Brillouin zone is here defined as $\text{BZ}=[-\pi,\pi]^{2}\times[-2\pi,2\pi]$.

II.1 Band structure fit with Monte Carlo sampling

Within the $19$-dimensional parameter space

[TABLE]

we seek the set $\{t\}$ that globally minimizes the quantity

[TABLE]

where $\mu=\alpha,\beta,\gamma$ is the band index, $\xi_{\mu}$ ($u_{C}^{\mu}$) is the band energy (orbital content $C=xy$, as given by the eigenvector components of $\pazocal{H}_{s=+1}$) of the model in Eq. (2). Similarly, $\tilde{\xi}_{\mu}$ ($\tilde{u}^{\mu}_{C}$) is the band energy (orbital content) of the model of Ref. Veenstra et al., 2014, which is based on a tight-binding fit from spin-resolved angle-resolved photoemission spectroscopy (ARPES) data. Finally, $w_{\boldsymbol{k}}$ ($\tilde{w}_{\boldsymbol{k}}$) are chosen energy (orbital) weights. For the $\boldsymbol{q}_{\mu}$’s we choose the three in-plane directions $\theta=0,\hskip 2.84526pt\pi/6,\hskip 2.84526pt\pi/4$ for the three $k_{z}$ cuts $0,\hskip 2.84526pt\pi,\hskip 2.84526pt2\pi$.

The hopping amplitudes were obtained by fitting the dispersion and orbital content of the 17-band model of Ref. Veenstra et al., 2014 using Monte Carlo (MC) sampling, see Fig. 1. 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$, $\Gamma$, $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 $\Gamma=\left(000\right)$, $Z=\left(00\frac{1}{2}\right)$, $R=\left(\frac{1}{2}0\frac{1}{2}\right)$, $X=\left(\frac{1}{2}00\right)$, $M=\left(\frac{1}{2}\frac{1}{2}0\right)$, $A=\left(\frac{1}{2}\frac{1}{2}\frac{1}{2}\right)$. The orbital weights $\tilde{w}_{\boldsymbol{q}}$ were fixed to be comparatively smaller than the energy weights $w_{\boldsymbol{k}}$.

The optimal tight-binding parameters are summarized in Tables LABEL:tab:HoppingParameters1 and LABEL:tab:HoppingParameters2. The expected uncertainties for the tight-binding parameters are $\pazocal{O}(10~{}\text{meV})$. Since the data from the 17-band model in Ref. Veenstra et al., 2014 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 tight-binding 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 $\beta$ and $\gamma$ sheets exhibit significant mixtures of the semi-2D orbtial ($d_{xy}$) and semi-1D Ru orbitals ($d_{xz}$ and $d_{yz}$), whereas the $\alpha$ sheet is dominated by the semi-1D orbitals. Figures 2 (d) and (e) compare two slices of the Fermi surfaces produced with the effective model with the model considered in Ref. Veenstra et al., 2014, and Fig. 2 (f) shows the full band structure of the latter model for completeness.

III Weak-coupling theory

The projection of the Coulomb interaction on the on-site $t_{2g}$ orbitals is given by

[TABLE]

where $i$ is the lattice site, $a$ is the orbital index, and $n_{ias}=c_{ias}^{\dagger}c_{ias}$ is the density operator. We assume that the phenomenological parameters satisfy $U^{\prime}=U-2J$ and $J^{\prime}=J$ Dagotto et al. (2001). 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. Kohn and Luttinger, 1965; Baranov and Kagan, 1992; Kagan and Chubukov, 1989; Chubukov and Lu, 1992; Baranov et al., 1992; Chubukov, 1993; Fukazawa and Yamada, 2002; Hlubina, 1999; Raghu et al., 2010a, b; Cho et al., 2013; Scaffidi et al., 2014; Šimkovic et al., 2016; Scaffidi, 2017. We first diagonalize $H_{0}$ and index the eigenstates by a band index $\mu=\alpha,\beta,\gamma$ and a pseudo-spin index, which we mostly keep implicit below. In this basis, the linearized gap equation reads

[TABLE]

where $S_{\nu}$ is the Fermi surface of band $\nu$, with $|S_{\nu}|$ the corresponding Fermi surface area, and $g$ the dimensionless matrix

[TABLE]

Here, $\Gamma$ is the two-particle interaction vertex (see Ref. Scaffidi, 2017 for details) at leading (second) order, $\rho_{\mu}=\lvert S_{\mu}\rvert/[\bar{v}_{\mu}(2\pi)^{3}]$ is the density of states, and $\bar{v}_{\mu}^{-1}=\int_{S_{\mu}}\frac{\mathrm{d}\boldsymbol{k}}{\lvert S_{\mu}\rvert}v_{\mu}(\boldsymbol{k})^{-1}$. Approaching the weak-coupling limit $U/t\to 0$ asymptotically, an eigenfunction $\varphi$ of the Eq. 11 corresponding to a negative eigenvalue $\lambda$ yields the superconducting gap

[TABLE]

below the critical temperature $T_{c}\sim W\ e^{-1/\lvert\lambda\rvert}$, where $W$ is the bare bandwidth.

Since we have chosen a pseudo-spin basis which is consistent with the tetragonal point group, each eigenvector $\varphi$ belongs to one of its ten irreducible representations $D_{4h}$ Sigrist and Ueda (1991); Raghu et al. (2010a); Annett (1990). In Table 3 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}(-\boldsymbol{k})=d_{0}(\boldsymbol{k})$, and odd-parity, $\boldsymbol{d}(-\boldsymbol{k})=-\boldsymbol{d}(\boldsymbol{k})$, components as Balian and Werthamer (1963)

[TABLE]

where $s$ and $s^{\prime}$ are pseudo-spin indices. Only the odd-parity $E_{u}$ and even-parity $E_{g}$ representations are two-dimensional 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 triplets (respectively, singlets). One should, however, keep in mind that a Zeeman field, as considered in Sec. IV.3, 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) Vaugier et al. (2012); Mravlje et al. (2011); Huo et al. (2013); Behrmann et al. (2012). However, recent work (Ref. Rømer et al., 2019b) reviewing a range of weak- and strong-coupling numerical methods applied to the two-dimensional Hubbard model found that both the $\boldsymbol{k}$-dependency 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 111This can be proven by noticing that the (even-parity) vertex at first order can be recast as a Gramian matrix.. 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. Scaffidi, 2017; Raghu et al., 2010a; Cho et al., 2013. 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 odd-parity helical order $A_{1u}$ for $J/U<0.15$, and the even-parity $d_{x^{2}-y^{2}}$ order $B_{1g}$ for $J/U>0.15$. Both types of order parameter are shown in Figs. 4, with $k_{z}$ cuts displayed in Figs. 5 and 6.

The results presented here were obtained by discretization of the Fermi surfaces, using $3318$ $\boldsymbol{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 Scaffidi (2017) of Eq. (12), and an effective diagonalization routine for the kinetic Hamiltonian constitutes the bulk of the numerical procedure Kopp (2008). Numerically, we computed sub-blocks of the $g$-matrix (Eq. (12)) simultaneously, using about $500$ cores over the duration of a few weeks.

IV.1 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 Scaffidi et al. (2014); Zhang et al. (2018), are present on all three bands. Notably, the gap on the $\beta$ band has $\min_{\theta}\lvert\Delta_{\beta}(\theta,k_{z}\approx\pi)\rvert\lesssim 0.02\Delta_{0}$, where $\theta$ is the in-plane azimuthal angle and $\Delta_{0}$ the maximal gap, at locations in agreement with previous predictions Firmo et al. (2013); Scaffidi et al. (2014); Dodaro et al. (2018). The character of the minima also appear likely to agree with thermal Hall measurements Hassinger et al. (2017); Dodaro et al. (2018). 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 $\beta$ and $\gamma$ 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$ Scaffidi et al. (2014); Zhang et al. (2018); Rømer et al. (2019a). 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 $\beta$ 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 Žutić and Mazin (2005). Interest in this state has been fueled by recent specific heat measurements Kittaka et al. (2018) combined with the possibility of explaining both TRSB and a nodal gap. However, this sector turns out to be strongly disfavored in our weak-coupling limit: At (e.g.) $J/U=0.20$, the best candidate has $\lambda_{E_{g}}/\lambda_{B_{1g}}\approx 0.03$ and thus does not come close to competing with the semi-two-dimensional order parameters found.

IV.2 Specific heat

The eigenvector calculated at weak coupling, $\varphi_{\mu}(\boldsymbol{k})$, is related to the order parameter via Eq. (13). We assume further that the order parameter factorizes as $\Delta(T)\Delta_{\mu}(\boldsymbol{k})$, with $\max_{\boldsymbol{k}}\Delta_{\mu}(\boldsymbol{k})=1$. The temperature dependency of $\Delta(T)$ is assumed to be that of a conventional BCS superconductor Tinkham (1975). The generalized BCS relation Sigrist (2005), fixing the overall size of the gap, when given the experimental value of $T_{c}\approx 1.48$ K, is

[TABLE]

where $\gamma\approx 0.5772$ is Euler’s constant, and where we introduced the average

[TABLE]

For a uniform gap, the average in Eq. (16) is 0, and the BCS relation $\Delta(0)=1.764~{}k_{B}T_{c}$ is recovered Tinkham (1975). For a multiband superconductor the specific heat per temperature per normal state value, $\gamma_{n}$, can be expressed as Sigrist (2005)

[TABLE]

where $E_{\mu}(\boldsymbol{k})=\sqrt{\xi^{2}+\Delta(T)^{2}\lvert\Delta_{\mu}(\boldsymbol{k})\rvert^{2}}$, and where the average here is evaluated as

[TABLE]

In Fig. 7 we show the specific heat as computed from Eq. 17 compared to experimental values from Ref. NishiZaki et al., 2000. 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 Nishizaki et al. (1998), there is no fine tuning of model parameters in these calculations (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 Nishizaki et al. (1998); Mackenzie et al. (2017); Scaffidi et al. (2014).

IV.3 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 Pustogow et al. (2019); Ishida et al. (2019).

We consider the normal state Hamiltonian with Zeeman terms when an external magnetic field is applied,

[TABLE]

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 $\boldsymbol{H}$ and the spin $\boldsymbol{S}$ (orbital $\boldsymbol{L}$) degrees of freedom (cf. Ref. Ramires and Sigrist, 2017),

[TABLE]

Here, $\varepsilon_{ab}(\boldsymbol{k})$ are the single-particle orbital terms, $\sigma^{a}$ is the $a$’th Pauli matrix, $\mu_{B}$ is the Bohr magneton (the vacuum permeability was fixed to $\mu_{0}=1$), $\epsilon_{abc}$ is the Levi-Civita symbol, and $a,b,c$ run over orbitals $xz,yz,xy$ ($A,B,C$, respectively).

The matrix $U(\boldsymbol{k})$ is defined to diagonalize the matrix $\pazocal{H}_{0}(\boldsymbol{k})+\pazocal{H}_{S}+\pazocal{H}_{L}$, where

[TABLE]

The (six) Fermi surfaces are defined by $\xi_{\mu\sigma}(\boldsymbol{k})=0$ for $\mu=\alpha,\beta,\gamma$ and $\sigma=+,-$. In the absence of a magnetic field the two energies $\xi_{\mu\pm}(\boldsymbol{k})$ become degenerate due to restored time-reversal symmetry. The matrix $U(\boldsymbol{k})$ and the energies $\xi_{\mu\sigma}(\boldsymbol{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 ($\mu$, $\sigma$) is given by the components of $U(\boldsymbol{k})$:

[TABLE]

We define the magnetization (with the magnetic field and the measured response along direction $i=x,y,z$) as

[TABLE]

Using Eq. (25) and $\langle c_{\mu\sigma_{1}}^{\dagger}(\boldsymbol{k})c_{\nu\sigma_{2}}(\boldsymbol{k})\rangle=\delta_{\mu\nu}\delta_{\sigma_{1}\sigma_{2}}f(\xi_{\mu\sigma_{1}}(\boldsymbol{k}))$, where $f$ is the Fermi function, the matrix elements of Eq. (27) are given by

[TABLE]

In the superconducting phase we add superconducting terms at orbital level

[TABLE]

[TABLE]

The electron operators in orbital and spin basis are now expressed as linear combinations of their particle ($u$’s) and hole ($v$’s) constituents,

[TABLE]

In turn, this leads to the magnetization matrix elements

[TABLE]

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

[TABLE]

where $\Delta_{\sigma_{1}\sigma_{2}}^{\mu}(\boldsymbol{k})$ is the order parameter in band and pseudospin basis, and where $u_{as}^{\mu\sigma}(\boldsymbol{k})$ are eigenvector components of $\pazocal{H}_{0}(\boldsymbol{k})$ (i.e. in the absence of a magnetic field, crucially with the same gauge choice as in the weak-coupling calculation).

We use the spin susceptibility normalized by its normal state value as a proxy for the Knight shift in the superconducting state,

[TABLE]

where the subscripts refer to the normal state value, $\Delta$ indicates a small increment in the external magnetic field in the linear response regime, and $\chi^{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 $\boldsymbol{k}$ the weak-coupling order parameter solution from the Fermi surface point closest to $\boldsymbol{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 Pustogow et al. (2019), which indicates a drop of around $K^{x}(T=20~{}\text{mK})\approx 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$) Rice and Sigrist (1995).

Recently, the experimental results of Ref. Pustogow et al., 2019 were reproduced by Ref. Ishida et al., 2019, 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 $d$-wave 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 compatible with specific heat data and recent Knight shift measurements Pustogow et al. (2019). 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 $\beta$ 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 even-parity order Kashiwaya et al. (2019). 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 Huang and Yao (2018); Huang et al. (2019); Rømer et al. (2019a). The near-degeneracy of the various odd-parity orders, found here and in previous work Scaffidi et al. (2014); Zhang et al. (2018), could potentially provide evidence for this scenario. We note that the four helical states in general are non-degenerate when $J/U\neq 0$, $\eta\neq 0$, and $\varepsilon_{AB}\neq 0$ Scaffidi et al. (2014). 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.

Acknowledgments

Helpful conversations with Fabian Jerzembeck, Clifford Hicks, Steven Kivelson, Yoshiteru Maeno, Daniel Agterberg, Assa Auerbach, Stephen Blundell, Mats Horsdal, Andrew Mackenzie, Peter Hirschfeld, and Catherine Kallin are acknowledged. We thank Kenji Ishida for providing the experimental data from Ref. Ishida et al., 2019. T.S. acknowledges support from the Emergent Phenomena in Quantum Systems initiative of the Gordon and Betty Moore Foundation. H.S.R.  and G.F.L.  are both supported by the Aker Scholarship. F.F. acknowledges support from the Astor Junior Research Fellowship of New College, Oxford. S.H.S. is supported by EPSRC Grant No. EP/N01930X/1.