Intrinsic Mechanism for Magneto-Thermal Conductivity Oscillations in Spin-Orbit-Coupled Nodal Superconductors

We describe a mechanism by which the longitudinal thermal conductivity $\kappa_{xx}$, measured in an in-plane magnetic field, oscillates as a function of field angle in layered nodal superconductors. These oscillations occur when the spin-orbit splitting at the nodes is larger than the nodal scattering rate, and are complementary to vortex-induced oscillations identified previously. In sufficiently anisotropic materials, the spin-orbit mechanism may be dominant. As a particular application, we focus on the cuprate high-temperature superconductor YBa$_2$Cu$_3$O$_{6+x}$. This material belongs to the class of Rashba bilayers, in which individual CuO$_2$ layers lack inversion symmetry although the crystal itself is globally centrosymmetric. We show that spin-orbit coupling endows $\kappa_{xx}/T$ with a characteristic dependence on magnetic field angle that should be easily detected experimentally, and argue that for underdoped samples the spin-orbit contribution is larger than the vortex contribution. A key advantage of the magneto-thermal conductivity is that it is a bulk probe of spin-orbit physics, and therefore not sensitive to inversion breaking at surfaces.


I. INTRODUCTION
Nodal superconductors are characterized by an energy gap that vanishes at point or line "nodes" on the Fermi surface. Low-energy quasiparticle excitations exist in the neighborhood of the nodes, and these excitations are reflected in characteristic power laws in the temperature dependence of various thermodynamic quantities. 1 These power laws can distinguish different node types (i.e. line versus point), but contain incomplete information about the symmetry of the superconducting state. Magnetothermal conductivity measurements are particularly useful in this regard, as the magnetic field dependence contains information about the k-space structure of the gap nodes. 1 A typical experiment involves the measurement of the thermal conductivity tensor κ ij in a magnetic field, which is swept through polar and/or azimuthal angles. In a nodal superconductor, the longitudinal thermal conductivity, for example κ xx , will oscillate as a function of the relative orientation of the field and the gap nodes. 2,3 The oscillation pattern has a nontrivial dependence on field strength and temperature, but with appropriate modeling will reveal the symmetry of the superconducting state. 3,4 Experimentally, this technique has been used to study the gap symmetry for a variety of unconventional superconductors, including organic 5 and heavy Fermion materials, 1,6,7 optimally doped YBa 2 Cu 3 O 7−δ , [8][9][10] and Sr 2 RuO 4 . 11 There are two established mechanisms underlying these oscillations, both of which are associated with the vortex lattice formed by the magnetic field. At low temperatures, circulating vortex currents "Doppler shift" the quasiparticle energies by an amount v s (r) · k, where k is the quasiparticle wavevector and v s (r) the superfluid velocity in the neighborhood of r. 2,12 The Doppler shift induces a nonzero density of states at each node that depends on the angle between the Fermi wavevector k F at that node and v s (r). The total induced density of states, and consequently the thermal conductivity, therefore changes with the orientation of the vortex lattice or, equivalently, the field angle. At high temperatures, a second mechanism takes over, namely anisotropic quasiparticle scattering by vortices becomes the dominant source of field-angle dependence. 3 In this work, we discuss a third mechanism that is especially relevant to layered superconductors in the highly anisotropic (quasi-two-dimensional) limit. In this geometry, an in-plane magnetic field generates only weak circulating vortex currents because of the small quasiparticle mobility along the interlayer direction. Under such circumstances, we show that one may still observe pronounced field-angle oscillations in the presence of spinorbit coupling (SOC). As a particular application of this mechanism, we focus on so-called "hidden spin-orbit" superconductors.
In the past few years, a number of materials have been discovered that exhibit signatures of SOC despite being both centrosymmetric and time-reversal symmetric. [13][14][15][16][17][18] Naive considerations would suggest this is not possible: by Kramers' theorem, materials that satisfy both inversion and time-reversal symmetry must have degenerate energies E k↑ and E k↓ , which seems to eliminate the possibility of k-space spin textures. Such spin textures, which are a hallmark of SOC, have nonetheless been observed. Key to this is that the Kramers-degenerate states are spatially separated, which leads to spin textures that are localized in space. [19][20][21] Rashba bilayers form a prominent subclass of hidden spin-orbit materials. In these materials, the unit cell con-tains pairs of conducting layers; while the unit cells are centrosymmetric, the individual layers are not. Rather, the layers are "inversion pairs", meaning that they transform into one another under inversion. 20,21 The individual layers thus exhibit some combination of Rashba and Dresselhaus SOC, with the Rashba contribution typically being larger in layered materials; 20 however, global inversion symmetry guarantees that spin textures in one layer are compensated for by opposite textures in the other layer. It is thus essential that the coupling between the layers be weak, as the k-space spin textures will be quenched when the two layers are strongly hybridized.
Rashba bilayers have been investigated as possible topological insulators, 22 as semimetals with electrically tunable Dirac cones, 23 and as model systems with nontrivial superconducting 24-27 and nematic 28,29 phases. Furthermore, many high-temperature superconductors, including YBa 2 Cu 3 O 6+x (YBCO 6+x ) and Bi 2 Sr 2 CaCu 2 O 8+x (Bi2212), satisfy the structural requirements to be Rashba bilayers; however, the relevance of this fact to cuprate physics is not established and hinges on the size of the effect.
Rashba-like spin polarization patterns have been directly measured in Bi2212 via spin-polarized angle-resolved photoemission spectroscopy (ARPES) experiments. 30 While important, these observations require independent confirmation because ARPES is a surface probe, and therefore sensitive to inversion symmetry breaking at surfaces. Indirect evidence for SOC has been obtained from the magnetic breakdown energy scale that one infers from many quantum oscillation experiments in YBCO 6+x . 31,32 However, it remains open whether the observed splitting is due to spin-orbit physics, 33,34 or to interlayer coupling. 35 A recent microscopic model for YBCO 6+x 36 suggests that the Fermi surface is spin-split by ∼ 10-20 meV due to a Rashba-like SOC; however, this is an upper bound as interlayer coupling may quench spin-orbit physics.
In this work, we show that the nodal structure of the d-wave superconducting gap allows for an elegant and straightforward observation of spin-orbit coupling through the longitudinal thermal conductivity in a transverse magnetic field. This effect will be present whether or not the superconductor is quasi-two-dimensional (quasi-2D), but must be disentangled from Doppler shift contributions if circulating vortex currents are not negligible. Importantly, this is a bulk measurement that is insensitive to inversion symmetry breaking at sample surfaces, and is complementary to other recent proposals: Kaladzhyan et al. showed that the dominant Friedel oscillation wavevectors associated with impurity scattering (as measured by scanning tunneling spectroscopy) reflect the spin-splitting of the Fermi surface and can be used to obtain the spin-orbit coupling constant, 37 while Raines et al. discussed the practicality of spin-Hall and Edelstein effects as probes of SOC. 38 We provide a simple description of the effect in Sec. II. While the mechanism has some similarities to the Doppler shift mechanism described in Ref. 2, there is the essential difference that a finite density of states is induced by Zeeman coupling to the quasiparticles rather than by circulating currents. Thermal conductivity calculations are described in Sec. III A, although details are left to the appendices, and results of these calculations are given in Sec. III B. We address the important question of how to distinguish spin-orbit and Doppler-shift contributions to the magneto-thermal conductivity in Sec. IV. In that same section, we make an estimate that suggests that magneto-thermal oscillations in YBa 2 Cu 3 O 6.5 are dominated by spin-orbit effects.

II. ORIGIN OF THE DENSITY OF STATES OSCILLATIONS
For a singlet superconductor, the BCS Hamiltonian takes the form where h.c. is the hermitian conjugate, ǫ k is the normalstate dispersion measured relative to the chemical potential, and χ k is the superconducting order parameter. The pairing term is typically simplified by making the permutation −c † k↓ c † −k↑ = c † −k↑ c † k↓ and recognizing that χ k = χ −k . However, Eq. (1) can be extended easily to include SOC, and is therefore left as-is for this discussion.
Written in this form, Eq. (1) generates four flavors of BCS quasiparticle, described by the creation operators where the coherence factors are and E k = ǫ 2 k + χ 2 k is the usual BCS quasiparticle excitation energy. The operators γ † kσs defined by Eqs. (2)-(5) are labeled by their spin-state σ and band index s = ±; we have followed the convention that the quasiparticle spectrum has two branches, with energies E kσ± = ±E k , corresponding to the quasiparticle operators γ † kσ± . The branches are independent of σ, and are thus doubly degenerate with the upper (lower) branches empty (fully occupied) at zero temperature.
These branches are pictured in Fig. 1(a) along two cuts through the Brillouin zone for the case of a nodal d x 2 −y 2 superconductor. The locations of the cuts are indicated (f) Spin-momentum locking on the normal-state Fermi surface determines the spin polarization of quasiparticle states near the gap nodes. Arrows correspond to the spin polarizations of the bands in (e). For node n, θn is defined as the smallest angle between the nodal spin polarization and the horizontal axis for the positive helicity band. (g) In an in-plane magnetic Zeeman field, the quasiparticle branches are shifted by −aµB cos(θn + φ), where φ is the angle between the field and the horizontal axis. (h) As a result, the sizes of the Bogoliubov Fermi surfaces are node-dependent. Note that the SOC (α = 40 meV) and the field strength (µB = 30 meV) are artifically inflated for clarity. Other parameters are as described in the main text.
in Fig. 1(b), which shows also the normal-state Fermi surface and gap nodes, i.e. points on the Fermi surface where χ k vanishes such that E k = 0. In Fig. 1(a), the excitation branches disperse linearly near the node, and are gapped away from the node.
The quasiparticles defined by Eqs. (2)-(5) have welldefined spins, such that the z-component of the spin operator isŜ A magnetic Zeeman field adds a term −gµ B BŜ z to the electronic Hamiltonian, with g the electronic g-factor and µ B the Bohr magneton. BecauseŜ z is diagonal in the quasiparticle operators γ kσs , this additional term leaves the quasiparticles intact, but rigidly shifts their dispersions by an amount −µB (σ =↑) or +µB (σ =↓), where µ = 1 2 gµ B is the electron dipole moment. These shifts are independent of the field direction, as the spin quantization axis is arbitrary in the absence of SOC, and is the same for each of the nodes. The resultant superconducting bands are shown in Fig. 1(c) along the same cuts as in Fig. 1(a). The band shifts inflate the nodal points to form so-called Bogoliubov Fermi surfaces 39 that separate occupied and empty quasiparticle states. These are shown as green ellipses in Fig. 1(d); they are the same for all nodes and are independent of field direction. 40 The situation changes when the SOC is nonzero. Here, the size of the induced Bogoliubov Fermi surfaces varies from node to node and depends on the field angle. For a generic 2D d x 2 −y 2 superconductor with Rashba SOC, the Hamiltonian isĤ = k where C k = ( c k↑ , c k↓ , c † −k↑ , c † −k↓ ) T , the prime indicates the summation is over a reduced Brillouin zone Here, h k is the Hamiltonian for the normal state, τ and τ 0 are the Pauli spin matrices, g k is the Rashba spin-orbit term, and B is the in-plane magnetic field. As shown in Fig. 1(f), SOC splits the normal-state Fermi surface into spin-polarized bands, with the electron spins locked to their momentum. For the d x 2 −y 2 superconductor, the nodes shown in Fig. 1(a) are shifted by displacements ±δk [Figs. 1(e) and (f)], where δk ∼ α/v F . This doubles the number of nodes in each quadrant of the Brillouin zone, but the dispersion near each of the shifted nodes has the same structure as when SOC is absent [ Fig. 1(e)]. (Note that although calculations are performed in a reduced Brillouin zone, we continue to show the full zone for illustrative purposes.) We remark that the superconducting order parameter ∆ k must develop a triplet component in response to the SOC, and that this may alter the structure of the gap nodes. For simplicity, we make the assumption that the triplet component is small and can be neglected, which is certainly the case in YBCO 6+x (see Appendix A).
In most SOC materials, the coupling constant α is orders of magnitude larger than µB ∼ 1 meV, and in this limit the physics of the nodal dispersion is easily understood. Crucially, the SOC selects a preferred polarization axis near each of the gap nodes, and this axis is largely unchanged by the Zeeman field in the limit µB ≪ α. In each quadrant (labeled n = 1, . . . , 4) one may locally rotate the spin-quantization axis such that the Hamiltonian for each band has a BCS-like form (Appendix B). The quasiparticle creation operators are then similar to Eqs. (2)-(5) near the gap nodes, but with "up" and "down" spin directions aligned with the red and blue arrows, respectively, in Fig. 1(f). For quadrant n, we denote the angle between the "up" direction and the k x axis by θ n .
In the limit µB ≪ α, the principal effect of the Zeeman field is to shift the nodal dispersions by an amount −aµB cos(θ n + φ), where φ is the angle between the magnetic field and the k x axis, and a = ± is the helicity of the band (positive helicity indicates that the spin winds clockwise around the center of the Brillouin zone). This has several consequences. First, the dispersions near the two spin-split nodes are shifted in opposite directions because their helicities are opposite [ Fig. 1 The sizes of the Fermi surface pockets [ Fig. 2(a)] are directly related to the induced density of states ∆ρ(ε F ) at the Fermi energy. Analytic expressions for ∆ρ(ε F ) may be obtained (Appendix C), and their field-angledependence, shown for each of the nodal regions in Fig. 2(b) is a signature of Rashba SOC. The two regions oscillate out of phase with each other, and the total induced density of states obtained from their sum has minima at field angles φ = (m ± 1 4 )π. This is reflected in the linear specific heat coefficient, shown in Fig. 2(c). The size of the oscillations depends on both B and on the single-particle scattering rate γ. For the benchmark case of YBa 2 Cu 3 O 6.5 shown in Fig. 2(c), the induced specific heat coefficient ∆γ is comparable to typical measured values for an out-of-plane magnetic field (i.e. for the vortex phase). 41 Significantly, the predicted magneto-thermal oscillations are a factor of 4 larger than the expected vortex contributions for an in-plane field. In Fig. 2(c), the vortex contribution is calculated using an estimate from Ref. 2 and is discussed in detail in Sec. IV.
To obtain quantitative estimates for YBa 2 Cu 3 O 6.5 , we have used a dispersion ǫ k that was obtained from a tightbinding fit to ARPES measurements on YBCO 6+x , 42 and a realistic gap function χ dk = χ d (cos k x −cos k y )/ √ 2 with χ d = 50 meV. These choices give a nodal Fermi velocity v F = 1.2 eVÅ, in close agreement with experiments, 43 and a superconducting nodal group velocity v 2 = |∇ k ∆ k | node = 0.17 eVÅ. Unless specified otherwise, the Rashba coupling constant is taken to be α = 10 meV throughout this work, which gives a spin-splitting at the gap nodes of 25 meV. 36 Although our explanation of the density of states oscillations is based on the weak-field limit, the effect is general provided the spin-splitting of the nodes is greater than the quasiparticle scattering rate γ. The effect is present regardless of the dimensionality of the system and will dominate over the vortex contribution in highly anisotropic materials; however, even in three dimensional materials the SOC and vortex effects can be comparable.
We finish this section with a comment that is specific to Rashba bilayers. The model explored in this section describes a single CuO 2 layer, and there are two issues that might limit its applicability to the bilayer. First, in a Rashba bilayer, the sign of α is opposite in each layer, and it is a concern that the contributions from each layer might cancel. However, from Fig. 2(a), it is clear that reversing the direction of each nodal polarization will have no effect on the induced density of states. Indeed, as we show explicitly below, the field-angle dependence of the thermal conductivity is an even function of α. Second, one must keep in mind that hybridization of the two layers via a hopping matrix element t ⊥ will quench the spin polarization. As discussed elsewhere, 33,34,36 the spin polarization at the gap nodes is of order α/ t 2 ⊥ + α 2 . The analysis contained in this work assumes that t ⊥ ≪ α, which is supported by an apparent collapse of bilayer splitting at the gap nodes in underdoped YBCO 6+x . 44 An experimental failure to measure the predicted thermal conductivity oscillations likely implies that the limit t ⊥ ≪ α does not apply.

A. Calculations
In this section, we discuss calculations of the longitudinal thermal conductivity in the presence of an inplane magnetic Zeeman field. As in the previous section, we assume that the triplet contribution to the superconducting order parameter can be neglected. For YBCO 6+x , we have checked numerically that neglect of the triplet components has no observable effect on the calculated longitudinal thermal conductivity. This is essentially different, then, from the intrinsic thermal Hall effect in a perpendicular Zeeman field, which depends crucially on the triplet component; together with SOC, a perpendicular Zeeman field creates a gapful mixed-parity topological superconductor 45-47 whose finite Chern number determines the T -linear part of the thermal Hall conductivity. 48 Following Refs. 49 and 50, the thermal current opera-tor is whereĊ k,i indicates a time derivative of C k,i , and where the velocity matrix is In this expression, v k = ∇ k h k and v ∆,k = ∇ k ∆ k , with h k and ∆ k given by Eqs. (9) and (10), respectively. Equation (13) does not include corrections due to circulating thermal currents that appear when time-reversal symmetry is broken, 51 as these do not contribute to the longitudinal thermal conductivity. From the Kubo formula, the longitudinal thermal conductivity satisfies where f (x) is the Fermi function, d is the mean interlayer distance (for the bilayer case of YBCO 6+x , it is the c-axis lattice constant divided by two), and where is the dimensionless thermal conductivity kernel. In Eq. (16), N k is the number of k-points in the reduced Brillouin zone, a 0 is the lattice constant, and A k (x) is the spectral function obtained from the Hamiltonian H k , with γ the quasiparticle scattering rate. At low T , Eq. (15) simplifies to In the limit of strong SOC, |α| ≫ µB, γ, it is further possible to obtain an analytic result for Π xx (x) (see Appendix B), with E(y) = y + 1 y tan −1 y. B. Results Figure 3 shows the longitudinal thermal conductivity in an in-plane magnetic field, as a function of the angle between the field and the x-axis. We focus initially on the large-SOC limit, with α = 10 meV, since this is the regime that we expect to be relevant to YBCO 6+x . The essential point of this figure is that at sufficiently low temperatures, κ xx /T exhibits clear and pronounced oscillations. The size of the oscillations is approximately proportional to µB/γ, and grows by an order of magnitude between the intermediate scattering [ Fig. 3(a) and (b)] and clean [ Fig. 3(c)] limits.
The oscillations change their qualitative character as a function of temperature, and in the clean limit [ Fig. 3(c)], the crossover between low-and high-temperature patterns occurs at k B T ∼ µB. This crossover is shifted to somewhat higher temperatures when γ ∼ µB [ Fig. 3(a)  and (b)]. The shape of the oscillations depends also on the ratio µB/γ, taking an approximately sinusoidal form when µB ∼ γ, and deviating strongly from it when µB ≫ γ.
The figure also shows that the approximate expression (19) for the transport kernel Π xx (x) works well when α is large. Inspection of Eq. (19) reveals that, in the zerotemperature limit, κ xx /T is a function of µB/γ, and is independent of α. This is apparent in Fig. 4(a), which shows that the oscillation amplitude at T = 0 saturates at an approximately constant value when |α| ≫ γ. We also note that the amplitude is a symmetric function of α, confirming our earlier assertion that the thermal conductivity oscillations due to the two layers making up a Rashba bilayer are additive. Figure 4(b) shows the dependence of κ xx /T on magnetic field strength. When µB → 0, the conductivity kernel becomes independent of the scattering rate, 50 yielding lim T →0 κ xx /T = 0.062 mW/K 2 ·cm. This value is close to that measured by Sutherland et al. in YBCO 6.54 . 52 When µB is not zero, the field dependence reflects the structure of the function E(µB/γ) in Eq. (19). For small argument, the field dependence is quadratic, with E(µB/γ) ≈ 1+ 2 3 (µB/γ) 2 , while for large arguments it is linear with E(µB/γ) ≈ πµB/2γ.
The value of γ is thus central to the observability of the thermal conductivity oscillations, which are suppressed when γ > µB. Angle-resolved photoemission spectroscopy (ARPES) experiments have placed an upper bound of γ = 12 meV on the nodal scattering rate for the bilayer cuprate superconductor Bi2212, 53 consistent with the strong inhomogeneity observed in that material in tunneling experiments. 54,55 Spin-orbit effects might thus be hard to observe in Bi2212. Conversely, microwave conductivity measurements have found very small transport scattering rates γ tr < ∼ 0.1 meV in YBCO 6.50 at low temperatures. 56 The reasonable assumption γ tr ∼ γ places YBCO 6.50 in the clean limit, where thermal conductivity oscillations should be easily observable. Eq. (19). Equation (19) is nominally valid for |α| ≫ µB but provides a good quantitative fit to the data even for α = µB.
Results are for T = 0, µB = 3 meV, and γ = 0.1 meV. Figure 5 shows the range of behavior that can be expected for κ xx /T in the clean limit for different values of α spanning |α| ≪ µB to |α| ≫ µB. For comparison, the analytic result for large SOC is also plotted, and is quantitatively similar to the numerical data for all values of α. Remarkably, there is very little to distinguish the different thermal conductivity data sets, even though α changes by a factor of 20 across the figure. Thus pronounced oscillations in the thermal conductivity should be observable so long as |α| > γ. By fitting Eq. (19) to experimental measurements of κ xx /T at low T and for µB < ∼ |α|, it is possible to extract both the ratio v F /v 2 and the scattering rate γ.

IV. DISCUSSION AND CONCLUSIONS
In this work, we have shown that spin-orbit coupling generates a characteristic field-angle dependence of the longitudinal thermal conductivity in nodal superconductors. These oscillations reflect how the density of states induced by a magnetic Zeeman field depends on the angle between the field and the spin-polarization at the gap nodes. Although we have focused on the case of a d x 2 −y 2 superconductor with Rashba SOC, the mechanisms described in this work will be present for any nodal superconductor in which spin-orbit physics leads to spinmomentum locking at the gap nodes. At low temperature (T < ∼ µB), the oscillation pattern is a function of the angle between the nodal spin axis and the magnetic field, and therefore depends both on the location of the nodes and the structure of the SOC. When the SOC is known, the magneto-thermal oscillations can be used to determine the locations of the nodes; conversely, when the node positions are known, the oscillations yield a fingerprint of the SOC. This analysis is simplest if there is a clean separation between spin-orbit and vortex contributions to the magneto-thermal oscillations. Indeed, the field-angle oscillations are qualitatively similar in both cases, and the question of how one may distinguish them is important. For d-wave superconductors, it is common to write the field-angle dependence of the thermal conductivity as a series, κ xx = κ 0,xx + κ 2 cos 2φ + κ 4 cos 4φ + . . .
The terms κ 2 and κ 4 are the amplitudes of the twofoldand fourfold-symmetric contributions to κ xx , respectively, and φ is the field angle as before. The fourfold term, κ 4 , is a direct consequence of the fourfold symmetry of the excitation spectrum. It is generally attributed to the symmetry of the gap function, but may also reflect the underlying band structure. 4 κ 4 may be positive or negative, with κ xx having its minimum at φ = π/4 in the first case, and its maximum at π/4 in the second. For oscillations due to the vortex lattice, the general trend at low fields is that κ xx is positive at low T , but then goes negative as T is increased. 3 Fig. 3 shows that the trend is similar here, with the crossover happening at k B T ≈ µB in the clean limit.
The T -dependence of κ 4 therefore does not allow us to isolate the source of the oscillations. However, the fieldstrength dependence is qualitively different for the two mechanisms, and does allow one to distinguish between them. In the clean limit (µB ≫ γ) and at low T , κ 4 grows as √ B for vortex-driven oscillations, 2 but as B for SOC-driven oscillations.
Another distinction between SOC-and vortex-driven magneto-thermal conductivity oscillations lies in κ 2 . In previous work, κ 2 was found to arise from the scattering of thermal currents by the vortices, and the twofold anisotropy reflects the difference between driving currents parallel to and perpendicular to the vortices. While in many materials this term is dominant at elevated temperatures, it is suppressed in quasi-2D materials where the circulating currents are small. Formally, κ 2 = 0 in our calculations, and any nonzero value of κ 2 would therefore signal a nonzero vortex contribution. A weak κ 2 in conjunction with a significant κ 4 is a strong hint that SOC is the dominant factor in observed magnetothermal conductivity oscillations.
One consequence of κ 2 vanishing is that, for the SOCdriven oscillations, the longitudinal thermal conductivity should be the same for both [ For purely SOC-driven oscillations the intrinsic contribution to κ xy /T in the limit T → 0 vanishes for a d-wave superconductor in an in-plane Zeeman field, 48,51 and the leading-order intrinsic contribution to κ xy must therefore be O(T 2 ). Since the leading-order contributions to κ xx and κ yy are O(T ), lim T →0 κ l reduces to κ xx and 2κ xx κ yy /(κ xx + κ yy ) for the [100] and [110] directions, respectively. These are equal for tetragonal superconductors.
For the case of YBCO 6+x , we note that the material itself is orthorhombic due to the presence of onedimensional CuO chains. This twofold anisotropy is tied to the crystal lattice, rather than the vortex lattice, and will show up primarily as a difference in the zeroth order term in the expansions of κ xx and κ yy , so that κ 0,xx = κ 0,yy .
We can estimate the relative importance of the vortex and SOC contributions for YBCO 6.5 by focusing on the specific heat. At this doping level, YBCO 6+x is strongly anisotropic, with the c-axis conductivity a factor of 10 3 smaller than the in-plane conductivity. The density of states induced by circulating vortex currents can be obtained from the clean-limit (γ → 0) approximation given in Ref. 2, with ρ 0 the normal state density of states and where Φ 0 = π /e is the superconducting flux quantum. In YBCO 6.5 , the ratio of in-plane and out-of-plane penetration depths λ c /λ ab ≈ 35, 57,58 and the fitted dispersion ǫ k gives the normal state density of states at the Fermi level ρ 0 = 3.35 eV −1 . The resultant shift in the specific heat constant, ∆γ 0 is shown in Fig. 1(d) and is considerably smaller than that obtained from SOC. It is reasonable to expect a similar disparity for the thermal conductivity, so that experimental observations of magnetothermal conductivity oscillations for in-plane fields, along with a small value for κ 2 , would be consistent with significant spin-orbit coupling. We finish with a few caveats. First, we note that a detailed quantitative description of any particular material will depend on the details of the scattering rate. In particular, we have taken a simple model in which γ depends on neither wavevector or energy. While such an assumption is sensible for many materials, it can be problematic in nodal superconductors where the scattering rate can have a nontrivial energy dependence that is determined by the strength of the scattering potential. In cuprates, for example, quantitatively accurate models typically require an admixture of Born and unitary scatterers. [59][60][61] Each of these has a characteristic energy dependence that will modify the temperature dependence of the oscillations. Rigorous modeling of experiments will thus require a realistic disorder model.
Second, the question of how to extract the size of the spin-orbit coupling constant α remains open. While the existence of SOC in a nodal superconductor is easy to establish via the thermal conductivity oscillations, the amplitude of the spin-orbit coupling constant is harder to determine. In clean materials, the size of the oscillations is nearly independent of α [ Fig. 4(a)] and there is no clear crossover in behavior as a function of magnetic field strength between µB ≪ α and µB ≫ α. Rather, the amplitude of the oscillations is an indication of whether the spin splitting of the Fermi surfaces is large or small relative to the scattering rate.
Third, we have neglected in this work variations of the chemical potential with the field angle, which is another potential contribution to the field-angle dependence of κ xx . In our analytical calculations, these variations would manifest themselves as modulations of both the Fermi velocity and the anomalous velocity, v 2 . We have checked numerically that, at least for the model used in this work, chemical potential modulations have a negligible effect on the thermal conductivity.
In summary, we have demonstrated the existence of a novel intrinsic mechanism for oscillations of the longitudinal magneto-thermal conductivity as a function of magnetic field angle. Zeeman coupling to nodal quasiparticles inflates the nodes into Bogoliubov Fermi surfaces, whose sizes depend on the field angle, the structure of the SOC, and the location of the gap nodes. The magnitude of the induced specific heat and thermal conductivity both depend on the sizes of the induced Fermi surfaces, and the angle-dependence of the thermal properties provides a tool to explore the SOC (if the locations of the gap nodes are known) or the nodal structure (if the SOC is known). As a probe of SOC, this technique has the advantage of being a bulk measurement, and is therefore insensitive to inversion-symmetry breaking at surfaces.
The structure of the magneto-thermal oscillations is qualitatively similar to what was found earlier for vortexinduced oscillations, but can be distinguished by details of the field and angle dependence. Interestingly, we find that the pattern of the oscillations inverts at high temperature, similar to what was reported earlier for vortexinduced oscillations. We consider the superconducting state of a single CuO 2 layer with Rashba SOC. The superconducting state is predominantly singlet, with a d-wave symmetry, but additional triplet pieces are induced by both the SOC and the in-plane Zeeman field. The goal of this section is to evaluate the size of the triplet components.
Taking the basis the Hamiltonian has the formĤ = k ′ C † k H k C k , where the primed sum is restricted to half of the Brillouin zone, and and whereg k = α(sin k y + i sin k x ) − µBe −iφ and φ is the angle between the x-axis and the in-plane magnetic field. We take the simplest form of pairing interaction appropriate for the cuprates, Under the assumption that the singlet order parameter is d-wave, we write χ dk = χ d η dk with Similarly, the triplet components have odd spatial parity, and can therefore be written with (a = x, y) We choose parameters such that the singlet component of the order parameter is χ d = 50 meV, which is comparable to the antinodal gap in underdoped YBCO 6+x . SOC induces a triplet component of the form ±d xk + id yk that goes along with the singlet piece. 62 For our model parameters, self-consistent calculations find that the triplet component is approximately 1% of the dominant singlet component, and is nearly independent of magnetic field strength and direction. In addition, the Zeeman field induces a second triplet component d zk that depends on field angle. This component is found to be three orders of magnitude smaller than the singlet component.

Appendix B: Thermal Conductivity: Large SOC Limit
In this appendix, we derive an analytic approximation for Π xx (x) that is valid in the limit α ≫ γ, µB. To evaluate Eq. (16), we transform both the Hamiltonian and the quasiparticle velocity operators to the helical basis. When the SOC is large, we can neglect the mixing of bands of different helicities, either by impurities or by the magnetic field. This simplification allows us to derive an explicit expression for Π xx (x).
a. The Hamiltonian in the helical basis. As a first step, we transform the Hamiltonian, Eq. (A2), to the helical basis via the unitary transformation, with e iθ k = g k /|g k | and g k = α(sin k y + i sin k x ). In zero field, U k diagonalizes the nonsuperconducting Hamiltonian h k , with U † k h k U k = diag(ξ k+ , ξ k− ), and ξ k± = ǫ k ± |g k |.
In the superconducting state and with nonzero in-plane magnetic field, we obtain the transformed Hamiltonian, The diagonal block is The off-diagonal block has a similar structure to Eq. (A4); however, it simplifies considerably if the the triplet components of the order parameter can be neglected. Then with ∆ ± k = −e ±iθ k χ dk . When B = 0, the two helical bands are not mixed by singlet superconductivity under the restriction of zero-momentum pairing, and the Hamiltonian H ξ k thus describes two independent superconducting bands, each of which has a BCS-like structure.
b. The velocity matrix in the helical basis. Next, we transform the matrix defined by Eq. (14) via The top left block transforms as while the top right block is where v ξ kα = ∇ k ξ kα , v a ∆,k = −e iaθ k ∇ k χ dk , and the SOC-related anomalous velocity is v an k = i(∇ k θ k )g k . c. Re-ordering the basis. To move forward, it is useful to rearrange the Hamiltonian and velocity into blocks with given helicity, The velocity operator becomes where V 0 k contains the diagonal blocks of V ξ k and V an k contains the off-diagonal blocks. The Hamiltonian is where a = ± labels the band helicity, and The matrix h ′ k mixes bands with different helicities and may be dropped when µB ≪ |α|.
d. Spectral function. To zeroth order in h ′ k , H ξ k has the energy eigenvalues where a = ± is the helicity index and s = ± indicates whether the quasiparticle branch is upward or downward dispersing. The Green's function is where a and b represent the different helicities, andz ± = z ± µB cos(θ k + φ) with z a complex frequency. To obtain the spectral function for (real) frequency x, we take In this expression, the δ functions are understood as Lorentzians, δ(x) = π −1 γ/(x 2 + γ 2 ). e. Thermal conductivity. We evaluate the transport kernel Π(x) from Eq. (16), using both the leading-order spectral function, Eq. (B18), and the leading-order velocity operator V 0 k defined in Eq. (B11), which neglects interband mixing due to the anomalous velocity term v an k . As noted before, this calculation is zeroth order in interband mixing, but not zeroth order in the magnetic field.
We take the Brillouin zone to be 0 < k x < π, −π < k y < π, and sum over four nodal regions [ Fig. 6(b)]; two of these correspond to a sum over helicty index a; the other two will be indicated by a sum over nodal index n = 1, 2. We use rotated momenta k 1 and k 2 near each node, with the understanding that they are rotated by 90 • between nodes in regions 1 and 2, and that the zero of the coordinate system is at the nodal point for each (n, a). The advantage of this definition is that we can write Brillouin zone, which consists of the region 0 < kx < π, −π < ky < π. Nodal points, at which the superconducting d-wave order parameter vanishes, are indicated by dots. Arrows located at each dot indicate the directions of the normal and superconducting group velocities. Local coordinate axes k1 and k2 are attached to each of the nodes, and are parallel to the velocities. Each of the four nodes is labeled by its quadrant (n = 1, 2) and its helicity (a = ±).
for each node, with the approximation that θ k can be treated as constant in the neighborhood of each node. Recalling that tan θ k = sin k x / sin k y , we obtain θ 1 = π/4 for node 1 and θ 2 = 3π/4 for node 2. Similarly, we will assume that the quasiparticle velocities depend only on the nodal index n, and are constant in the vicinity of each node.
Under these transformations, where a 0 is the lattice constant. The kernel for the thermal conductivity is therefore where we have made the linearized nodal approximation, ξ = E cos ζ, ∆ a = Ee iaθn sin ζ and ζ ∈ [0, 2π]. The trace in these equations is over particle-hole channels associated with superconductivity, so the matrices are 2 × 2.
The velocity matrix v x n,aa is the 2 × 2 matrix obtained from the top-left (a = +) or bottom-right (a = −) block of V 0 k . The superscript x refers to the component of the quasiparticle velocity along the x-direction. Noting that v In linearized nodal coordinates, the spectral function in Eq. (B22) is Separating the thermal conductivity kernel into contributions from the two helicity bands, Π xx (x) = a Π a (x), we obtain for the + helicity band Using the Lorentzian form for the delta-function, the energy integrals easily give Performing the matrix multiplications, taking the traces, and integrating over the angle ζ gives where A nearly identical calculation gives the contribution Π − (x) to the thermal conductivity kernel from the − helicity bands. Π − (x) has teh same form as Eq. (B28), but with µB → −µB. The total kernel is then The thermal conductivity follows from In the limit T → 0, this expression simplifies to µB cos( 3π 4 + φ) γ .