Superconducting gap anisotropy and topological singularities due to lattice translational symmetry and their thermodynamic signatures

TL;DR

This paper investigates how lattice translational symmetry affects the superconducting gap structure and thermodynamic signatures in chiral superconductors, revealing challenges in interpreting experimental data and supporting a chiral $f$-wave pairing in Sr$_2$RuO$_4$.

Contribution

It provides a symmetry-based analysis of gap anisotropy due to lattice translations and compares thermodynamic signatures of chiral and nodal superconducting states.

Findings

Strong anisotropy in chiral states can mimic nodal behavior in thermodynamic measurements.

Available measurements are consistent with a chiral $f$-wave order parameter.

Lattice translation effects complicate the identification of the superconducting order parameter.

Abstract

Symmetry arguments based on the point group of a system and thermodynamic measurements are often combined to identify the order parameter in unconventional superconductors. However, lattice translations, which can induce additional momenta with vanishing order parameter in the Brillouin zone, are neglected, especially in gap functions otherwise expected to be constant, such as in chiral superconductors. After a general analysis of the symmetry conditions for vanishing gap functions, we study the case of chiral $p$- and chiral $f$-wave pairing on a square lattice, a situation relevant for Sr$_2$RuO$_4$. Specifically, we calculate the impurity-induced density of states, specific heat, superfluid density and thermal conductivity employing a self-consistent T-matrix calculation and compare our results to the case of a nodal ($d$-wave) order parameter. While there is a clear distinction between a fully gapped chiral state and a nodal state, the strongly anisotropic case is almost indistinguishable from the nodal case. Our findings illustrate the difficulty of interpreting thermodynamic measurements. In particular, we find that the available measurements are consistent with a chiral ($f$-wave) order parameter. Our results help to reconcile the thermodynamic measurements with the overall picture of chiral spin-triplet superconductivity in Sr$_2$RuO$_4$.

Equations

$$g:{\Delta }_{\mathbf{k}}\mapsto {\hat{U}}_g{\Delta }_{U_g^{-1}\mathbf{k}}\stackrel{!}{=}{e}^{i{\varphi }_g}{\Delta }_{\mathbf{k}},$$

$$U_g^{-1}\mathbf{k}+\mathbf{G}=\mathbf{k}$$

$$\mathcal{H}=\sum _{i,j,s}{t}_{ij}{a}_{i,s}^{†}{a}_{j,s}^{}+V_{\text{pair}},$$

$${\xi }_{\mathbf{k}}=-2t_1(\cos k_x+\cos k_y)-4t_2\cos k_x\cos k_y-\mu ,$$

$$V_{\text{pair}}=\sum _{\begin{array}{c}i,j\\ s_1\dots s_4\end{array}}{V}_{i,j}{\sigma }_{s_1{s}_2}^x{\sigma }_{s_3{s}_4}^x{a}_{i,s_1}^{†}{a}_{j,s_2}^{†}{a}_{j,s_3}^{}{a}_{i,s_4}^{},$$

$$V_{\text{pair}}=\sum _{\begin{array}{c}\mathbf{k},{\mathbf{k}}^{\prime }\\ s_1{s}_2\end{array}}{V}_{\mathbf{k}{\mathbf{k}}^{\prime }}{c}_{\mathbf{k},s_1}^{†}{c}_{-\mathbf{k},-s_1}^{†}{c}_{-{\mathbf{k}}^{\prime },-s_2}^{}{c}_{{\mathbf{k}}^{\prime },s_2}^{},$$

$$V_{\mathbf{k}{\mathbf{k}}^{\prime }}=\sum _{a=x,y}\left[V_1{\Phi }_{1a}(\mathbf{k}){\Phi }_{1a}({\mathbf{k}}^{\prime })+V_2{\Phi }_{2a}(\mathbf{k}){\Phi }_{2a}({\mathbf{k}}^{\prime })\right],$$

$${\Phi }_{1x}(\mathbf{k})$$

$${\Phi }_{2x}(\mathbf{k})$$

$$\begin{gathered} {\Delta }_{\mathbf{k}}={\Delta }_p(\sin k_x\pm i\sin k_y) \\ +{\Delta }_f(\sin k_x\cos k_y\pm i\sin k_y\cos k_x), \end{gathered}$$

$${\Delta }_{\mathbf{k}}=-T\sum _n\sum _{{\mathbf{k}}^{\prime }}{V}_{\mathbf{k}{\mathbf{k}}^{\prime }}\frac{{\Delta }_{{\mathbf{k}}^{\prime }}}{{\omega }_n^2+{\xi }_{{\mathbf{k}}^{\prime }}^2+{\Delta }_{{\mathbf{k}}^{\prime }}^2},$$

$$\left(\begin{array}{c}{\Delta }_p\\ {\Delta }_f\end{array}\right)=\sum _{\mathbf{k}}{\mathcal{C}}_{\mathbf{k}}\left(\begin{array}{cc}{V}_1& V_2\cos k_y\\ V_1\cos k_y& V_2{\cos }^2{k}_y\end{array}\right)\left(\begin{array}{c}{\Delta }_p\\ {\Delta }_f\end{array}\right).$$

$${\mathcal{C}}_{\mathbf{k}}=-T\sum _n\frac{{\sin }^2{k}_x}{{\tilde{\omega }}_n^2+{\xi }_{\mathbf{k}}^2+|{\Delta }_{\mathbf{k}}{|}^2}.$$

$$N_C=\frac{1}{2\pi }{\oint }_{C}d\mathbf{k}\cdot {\nabla }_{\mathbf{k}}{\theta }_{\mathbf{k}},$$

$$0$$

$${\mathbf{\eta }}_p$$

$${\mathbf{\eta }}_f$$

$$f$$

$$+\frac{b_2}{2}({\eta }_x^{\ast 2}{\eta }_y^2+{\eta }_x^2{\eta }_y^{\ast 2})+b_3|{\eta }_x{|}^2|{\eta }_y{|}^2,$$

$${\mathbf{\eta }}_p$$

$${\mathbf{\eta }}_f$$

$$f$$

$$+4\gamma {\eta }_p{\eta }_f\cos (\Delta \theta )+b_p{\eta }_p^4+b_f{\eta }_f^4$$

$$+4c_1{\eta }_p^3{\eta }_f\cos (\Delta \theta )+4c_2{\eta }_p^2{\eta }_f^2\left[1+\frac{1}{2}\cos (2\Delta \theta )\right]$$

$$+4c_3{\eta }_p{\eta }_f^3\cos (\Delta \theta ),$$

$${\mathcal{H}}_{\text{imp}}=\sum _{\mathbf{k},{\mathbf{k}}^{\prime },s}{U}_{\mathbf{k},{\mathbf{k}}^{\prime }}{c}_{\mathbf{k}s}^{†}{c}_{{\mathbf{k}}^{\prime }s}^{},$$

$$T_{\mathbf{k},{\mathbf{k}}^{\prime }}(i{\omega }_n)$$

$$T_{\mathbf{k},{\mathbf{k}}^{\prime }}(i{\omega }_n)=T(i{\omega }_n).$$

$$\Sigma (i{\omega }_n)=cT(i{\omega }_n)$$

$$i{\tilde{\omega }}_n=i{\omega }_n-\Sigma (i{\omega }_n).$$

$$G(\mathbf{k},i{\omega }_n)$$

$$F(\mathbf{k},i{\omega }_n)$$

$$F^{†}(\mathbf{k},i{\omega }_n)$$

$$\displaystyle\xi^{2}=\frac{\sum_{\mathbf{k}}\Big{|}\nabla_{k}\frac{\Delta_{\mathbf{k}}}{E_{\mathbf{k}}}\Big{|}^{2}}{\sum_{\mathbf{k}}\Big{|}\frac{\Delta_{\mathbf{k}}}{E_{\mathbf{k}}}\Big{|}^{2}}.$$

Full text

I Introduction

Generalizing unconventional superconductivity from the isotropic case to a lattice is usually done through the analysis of the generating point group. While the resulting classification is in terms of the irreducible representations of the point group, the nomenclature of the isotropic case, in other words in terms of angular momentum, is still generally used for simplicity and corresponds to the long-wavelength behavior of the order parameter. A result of such a symmetry analysis is the imposition of zeros in the gap function, resulting in a nodal structure of the excitation spectrum, with point or line nodes, depending on the symmetry properties under rotations or mirror operations. In two dimensions, for example, a $d$-wave order parameter has point nodes on either both axes or diagonals, while a $p$-wave order parameter only has point nodes on one axis. However, combining order parameters of a two-dimensional irreducible representation into a chiral state, these nodes can be lifted and a fully gapped state would be expected.

Experimentally, the nodal structure of a superconductor and the resulting low-energy excitations lead to distinct low-temperature thermodynamic responses, in particular power-law behavior in the temperature. These include the specific heat, superfluid density or London penetration depth, and thermal conductivity, all of which are usually employed to characterize the superconducting state and its symmetry.

However, the above described analysis neglects another set of symmetries of a lattice, namely translation symmetry. Translations in real space lead to translational symmetry also in momentum space and thus, the Brillouin zone (BZ). Furthermore, the translation symmetries in momentum space imprint further symmetry requirements onto the order parameter, which can result in additional zeros in the gap function as defined over the whole BZ. This is particularly important in systems, where no other nodes are expected, such as chiral superconductors.

Such considerations have come into focus in Sr2RuO4, where strong evidence for unconventional superconductivity has been gathered Maeno et al. (1994, 2001). Early on, the superconducting state was predicted to be of chiral $p$-wave symmetry Rice and Sigrist (1995), characterized by the $d$ vector $\mathbf{d}(\mathbf{k})=\boldsymbol{\hat{z}}\Delta(T)(\sin k_{x}\pm i\sin k_{y})$ Mackenzie and Maeno (2003); Maeno et al. (2011); Kallin (2012); Ishida et al. (1998); Luke et al. (1998, 2000); Duffy et al. (2000); Nelson et al. (2004); Xia et al. (2006); Ishida et al. (2015). Still, several discrepancies exist 111 Note that due to recent Knight shift measurements Pustogow et al. (2019), the assumed spin-triplet pairing symmetry of the superconducting state of Sr2RuO4 has been questioned once again., amongst them low-temperature measurements of thermodynamic quantities, such as specific heat NishiZaki et al. (1999, 2000); Deguchi et al. (2004) and the London penetration depth Bonalde et al. (2000), or thermal conductivity measurements Tanatar et al. (2001a, b); Izawa et al. (2001); Hassinger et al. (2017). These measurements all indicate nodes in the superconducting order parameter and thus seemingly contradict the proposal of an intrinsically nodeless $p$-wave state in Sr2RuO4. Note, however, that several theoretical works predict a substantial contribution of higher-angular-momentum spin-triplet states in Sr2RuO4, which manifest in strong gap anisotropies Raghu et al. (2010); Wang et al. (2013); Scaffidi et al. (2014).

Motivated by the general observation of gap anisotropies unexpected by pure point group analysis and the specific case of Sr2RuO4, we analyze here the occurrence of gap anisotropy in chiral superconductors and its connection to the thermodynamic response in the presence of disorder. For this purpose, we employ a tight-binding description of a two-dimensional lattice and calculate microscopically the density of states, specific heat, superfluid density and thermal conductivity. The disorder consists of non-magnetic impurities and the respective scattering is assumed to be isotropic and unitary, which allows us to employ a self-consistent T-matrix approximation.

To analyze all translation-symmetry-imposed zeros, we include nearest-, as well as next-nearest-neighbor interactions. This leads to an almost isotropic ($p$-wave) and a strongly anisotropic ($f$-wave) pairing gap. For comparison, we further add a nodal spin-singlet order parameter ($d$-wave). We observe that the strongly anisotropic case is hardly distinguishable from the nodal case once disorder is added.

The present article is structured as follows. In Sec. II we first discuss the origin and properties of symmetry-imposed zeros in a chiral gap function. We elaborate on the resulting gap anisotropy for the lowest order spin-triplet states on the square lattice. We further analyze the movability of the zeros within the framework of Ginzburg-Landau (GL) theory. In Sec. III we add disorder to our system and examine how the superconducting state is affected. Furthermore, in Sec. IV, we study the influence of the gap anisotropy on several thermodynamic properties. Finally, we discuss our findings and conclude in Sec. V.

II Origin and properties of nodes in the chiral spin-triplet state

II.1 General Symmetry considerations

Throughout this paper, we use the following nomenclature: Defining the gap function over the whole BZ, we refer to zeros, where the gap function vanishes. The transformation properties of a gap function under the symmetry operations of the generating point group can enforce zeros along high-symmetry lines or at high-symmetry points. Nodes, on the other hand, denote zeros in the excitation spectrum, in other words zeros of the gap function coinciding with the normal state Fermi surface.

Given the point group $\mathcal{G}$ of the underlying lattice, the superconducting order parameter transforms under an element $g\in\mathcal{G}$ as

[TABLE]

where $U_{g}$ and $\hat{U}_{g}$ are representations of $g$ in k space and the space of the gap matrix, usually spin or orbital space. The right-hand side of Eq. (1) with $\phi_{g}\in[0,2\pi)$ follows from the fact that $\Delta_{\mathbf{k}}$ transforms as an irreducible representation of $\mathcal{G}$ 222 Note that in case of a higher-dimensional irreducible representation, the actual gap function breaks additional symmetries and thus always transforms as a one-dimensional irreducible representation of the remaining symmetry group.. For any non-trivial phase $\phi_{g}$, this transformation behavior defines zeros in the gap function for any ${U_{g}^{-1}\bf k}={\bf k}$. Examples include the vanishing of the gap function at the $\Gamma$ point for any angular-momentum channel $l\neq 0$ with $g$ a rotation or along the main axis or diagonals for a $d$-wave order parameter with $g$ a mirror operation.

A lattice, however, obbeys additional symmetries, namely translations. In addition to the ones discussed above, Eq. (1) then indicates further zeros in the gap function for any momentum satisfying

[TABLE]

with $\mathbf{G}$ a reciprocal lattice vector if $\phi_{g}\neq 0$. To name an example, this condition implies a vanishing order parameter on the BZ boundary for the case of $d_{xy}$-wave pairing on a square lattice. More interestingly, the condition implies points of zeros for chiral superconductors as we will discuss in the following using specific lattice realizations.

To illustrate the appearance of zeros in the gap function and the most important features in the topological nature of the chiral superconducting phase, we first introduce a tight-binding model on the square lattice. The pairing interaction is formulated on the lattice, too, with nearest- and next-nearest-neighbor pairing. While we use this model primarily to study the chiral spin-triplet order parameters, it can also be used to study a nodal spin-singlet gap function for comparison in Sec. IV. Due to the lattice formulation, the structure of the gap function is adapted to the Brillouin zone (BZ) with a corresponding nodal structure. Note that the same arguments apply to lattices with a three-fold rotation axis, as we will discuss at the end of this section.

II.2 Two-dimensional square-lattice model

In the following, we study a system described by a Hamiltonian including nearest-neighbor (NN) and next-nearest-neighbor (NNN) hopping,

[TABLE]

where $a^{\dagger}_{i,s}$ ($a^{\phantom{\dagger}}_{i,s}$) denotes the creation (annihilation) operator for an electron with spin $s=\uparrow,\downarrow$ on site $i=(x_{i},y_{i})$. The dispersion in $k$ space is then given by

[TABLE]

where we have introduced the chemical potential $\mu$ and the first BZ is given by $k_{x},k_{y}\in[-\pi,\pi]$ taking the lattice constant $a=1$. In the following, we set the NN-hopping $t_{1}=1$ and the NNN-hopping $t_{2}=0.3$. We restrict the pairing potential $V_{\text{pair}}$ to the spin-triplet channel with (in-plane) equal-spin pairing,

[TABLE]

where $\sigma^{x}$ is a Pauli matrix, $V_{ij}=V_{1}$ for NN-pairing, and $V_{ij}=V_{2}$ for NNN-pairing. In momentum space, we obtain

[TABLE]

where

[TABLE]

with basis funtions

[TABLE]

These represent the simplest odd-parity pairing states on a square lattice with $\Phi_{ia}(-\mathbf{k})=-\Phi_{ia}(\mathbf{k})$, $i=1,2$.

Within standard mean-field theory, the pairing interaction results in the quasiparticle gap function

[TABLE]

which is a time-reversal-symmetry breaking phase with chiral character. For convenience, we refer to the NN part as $p$-wave and the NNN part as $f$-wave component. The coefficients $\Delta_{p,f}$ are obtained through the solution of the self-consistency equation Mineev et al. (1999),

[TABLE]

where $V_{\mathbf{k}\mathbf{k^{\prime}}}$ is the pairing potential in momentum space for the spin-triplet channel, $\omega_{n}=(2n+1)\pi k_{B}T$ are the fermionic Matsubara frequencies and $T$ denotes the temperature. Using the full gap $\Delta_{\mathbf{k}}$, Eq. (10), we can trace out the two basis functions $\Phi_{ix}({\mathbf{k}})$ to find two coupled equations, one for $\Delta_{p}$ and one for $\Delta_{f}$,

[TABLE]

The factor $\mathcal{C}_{\mathbf{k}}$ is given by

[TABLE]

For a given strength of the pairing potentials for $p$- and $f$-wave pairs, $V_{1}$ and $V_{2}$, we can self-consistently tune the values of $\Delta_{p}$ and $\Delta_{f}$.

II.3 Singularities and topological properties of the gap function

The NN- and NNN-pairing states in Eq. (8) and (9) have to satisfy several symmetry properties due to underlying lattice symmetries. In App. A, we use the previous example of spin-triplet pairing on a square lattice to show in detail how the singularities of the gap function follow from general symmetry considerations alone. In this section however, we start the discussion directly with the quasiparticle gap function [Eq. (10)].

To each of the zeros in $\Delta_{\mathbf{k}}$ we can attribute a topological charge corresponding to the winding number

[TABLE]

where $\theta_{\mathbf{k}}=\arg[\Delta_{\mathbf{k}}]$ and $C$ denotes a closed path in positive direction around the zero without enclosing any other zero. Considering $p$- and $f$-wave order parameters, all zeros possess a charge of either $+1$ or $-1$. The relations (32) and (33) connect zeros of the same charge. The zeros and their charge are shown in Fig. 1. Note that the charges change sign under time reversal operation.

The charges defined above are connected to the topology of the superconducting phase. For chiral superconducting states in two dimensions, the topological invariant is the Chern number, which is defined through Eq. (14) by taking the (normal-state) Fermi surface as the path $C$. The Chern number is then the sum over all charges encircled by the Fermi surface.

Figures 2 a) and b) show the gap function magnitude for the chiral $p$-wave state and the phase $\theta_{\mathbf{k}}$ for the state with positive chirality along the normal-state Fermi surface for different band fillings. While the gap function has no nodes on the Fermi surface, it shows strong minima for fillings, where the Fermi surface comes close to the van Hove points at $(\pm\pi,0)$ and $(0,\pm\pi)$, here for $\mu=0.47$. Lowering the band filling leads to a more and more isotropic gap. For all band fillings shown the phase adds up to a total winding of $2\pi$ and the Chern number is +1 for our choice of chirality. Interestingly, the angle dependence of the phase is non-monotonic for band fillings for which the Fermi surface crosses the lines connecting the van Hove points.

Figures 2 c) and d) display the gap anisotropy and phase for the chiral $f$-wave state. Again, we find a pronounced dip in the gap magnitude when the Fermi surface approaches the van Hove points ($\mu=0.47$). In addition, however, even more striking near nodes appear when the Fermi surface approaches the gap zeros at $(\pi/2,\pi/2)$. Further, we see how the phase now changes depending on whether the FS encompasses these points or not. If it does, the Chern number is $-3$, as one would expect from a chiral $f$-wave state in the isotropic case, while it is $+1$ otherwise.

Having both components $\Delta_{f}$ and $\Delta_{p}$, the zeros corresponding to Eq. (35) can be shifted, while the zeros of Eq. (31) remain fixed at the BZ boundary. The former zeros are located at the intersection of four lines satisfying the equations

[TABLE]

such that the position along the $[11]$ direction depends on the ratio of $\Delta_{p}$ to $\Delta_{f}$. Note, however, that the winding numbers associated with the zeros remain. For $\Delta_{p}=\Delta_{f}$, the gap function develops zeros along the BZ boundary or main axes allowing for the winding numbers to change. Depending on the Fermi surface topology, this is accompanied with a topological transition in the superconducting state.

Before continuing with the analysis of the topological transition, we want to emphasize that the given discussion of the singularities and their topological properties can be applied to other systems as well. As an example, in App. B, we comment on the situation of lattices with a three-fold rotation axis.

II.4 Topological transition

In the following we will be interested in the behavior of superconductors when the zeros described by Eq.(15,16) for $|\Delta_{p}|<|\Delta_{f}|$ lie close to the Fermi surface and give rise to a near nodal spectrum. Therefore, it is important to know whether a near nodal situation is stable, if we give the positions of the zeros some flexibility combining $p$- and $f$-wave components. Within a weak-coupling picture one might expect that the order parameter would change in a way as to avoid the zero near the Fermi surface, since it diminishes the condensation energy.

Here, we analyze this question using a generalized Ginzburg-Landau theory for the two order parameters, whose coefficients are derived from the microscopic weak-coupling theory assuming a variation of the chemical potential to shift the Fermi surface as could be done by doping or applying uniform pressure. The zero passing through the Fermi surface corresponds to a topological phase transition changing the Chern number.

The order parameter of the superconducting state for both $p$- and $f$-wave have two complex components and belong to the $\Gamma_{5}^{-}$ representation of the tetragonal point group,

[TABLE]

The Ginzburg-Landau free energy density for a simple two-component order parameter belonging to $\Gamma_{5}^{-}$ is given by Sigrist and Ueda (1991)

[TABLE]

where we assume homogeneity and ignore gradient terms. Since both order parameters, $\boldsymbol{\eta}_{p}$ and $\boldsymbol{\eta}_{f}$ belong to $\Gamma_{5}^{-}$, they couple in a complex way. To simplify the problem, we assume for both order parameters a chiral state and make the ansatz

[TABLE]

leading to the free energy expansion,

[TABLE]

where $\Delta\theta=\theta_{p}-\theta_{f}$ and $b_{p,f}=4b_{1p,f}-b_{2p,f}+b_{3p,f}$. The coefficients as a function of the chemical potential $\mu$ are obtained through the expansion of the self-consistent gap equation [Eq. (11)].

Following the equilibrium state for varying $\mu$ we see in Fig. 3 that no anomalous behavior appears when the zero crosses the Fermi energy. The critical $\mu$ of the topological transition is indicated by the dashed line, where the Chern number changes from 1 to - 3. Although both components change and shift the zeros slightly closer towards the $\Gamma$ point, they vary smoothly. This indicates that the zeros can be located arbitrarily close to the Fermi surface.

III Disorder and gap anisotropy

Unconventional Cooper pairing is very susceptible to scattering even on non-magnetic impurities. Before analyzing the thermodynamic properties of the individual pairing possibilities, we thus consider how disorder affects the two spin-triplet pairing states of $p$- and $f$-wave symmetry, introduced above, and combinations thereof.

III.1 Theoretical framework

We introduce disorder as identical impurity scatterers with concentration $c$. A single impurity at the origin is described by

[TABLE]

where $U_{\mathbf{k},\mathbf{k^{\prime}}}$ is the scattering matrix element. Since we consider strong scattering potentials, whose energy scale exceeds the band width, we employ a T-matrix approach, which takes multiple scatterings at the same impurity into account Hirschfeld et al. (1986, 1988); Balatsky et al. (2006). The T-matrix is defined through the equation

[TABLE]

where $G(\mathbf{k},i\omega_{n})$ is the (normal) electron Green’s function. Note that we have omitted terms containing $\sum_{\mathbf{k}}F(\mathbf{k},i\omega_{n})$, where $F(\mathbf{k},i\omega_{n})$ is the anomalous Green’s function, since these terms vanish for non-$s$-wave gap functions Hirschfeld et al. (1988); Keller et al. (1988); Nomura (2005). Consequently, there are no off-diagonal entries in the T-matrix expansion, which would renormalize the gap directly. In the following, we consider isotropic, $s$-wave scattering for a point-like impurity potential, removing the momentum dependence from the scattering matrix elements $U_{\mathbf{k},\mathbf{k^{\prime}}}=U$. Thus, the T-matrix is a scalar in momentum space as well,

[TABLE]

Restricting our investigation to small impurity concentrations $c$, we can safely neglect impurity interference effects coming from multiple scatterings on different impurities. Hence, the self-energy is proportional to $c$,

[TABLE]

and can be absorbed into a renormalized Matsubara frequency

[TABLE]

The normal and anomalous Green’s functions and the self-energy are then related through Gorkov’s equations Mineev et al. (1999); Nomura (2005). Assuming a flat density of states over an energy scale larger than the gap size, we obtain the following simplified form for the normal and anomalous Greens functions,

[TABLE]

where Eq.(25) has to be determined self-consistently. Note that the assumption of a flat density of states is generically not justified, but still affects results for the thermodynamic response considered below in a negligible way Monien et al. (1987); Nomura (2005). Using the renormalized Matsubara frequencies in the self-consistent gap equation [Eq. (11)], we can examine the influence of disorder on the various spin-triplet pairing states.

III.2 Influence on gap anisotropy

If we consider the two spin-triplet pairing states separately, disorder treated within the above formalism merely leads to a gradual shape-conserving suppression of the gap magnitude ($\Delta_{p}$ or $\Delta_{f}$) with increasing $c$. However, we can study the influence of disorder on the gap anisotropy by studying a combination of $p$- and $f$-wave state. For this purpose, we treat simultaneously NN- and NNN-pairing interactions self-consistently, which results in a mutual interaction between the two pairing channels. For the chemical potential, we choose $\mu=0.1$, such that the Fermi surface passes rather close to the $f$-wave zeros at $\mathbf{k}=(\pm\pi/2,\pm\pi/2)$.

In Fig. 4, we show the relative change of the ratio $\Delta_{p}/\Delta_{f}$ for two sets of pairing interactions. Both the temperature and the impurity concentration change this ratio. For a dominant $p$-wave component, we find a decrease of the ratio with temperature from $\Delta_{p}/\Delta_{f}\approx 2.6$ at $T_{c}$ to $1.9$ at $T_{0}=T_{c}/50$ in the clean limit. Interestingly, the ratio decreases with disorder at $T_{c}$ while this trend is reversed at low temperature. At low temperature, the change can be understood from the difference in the coherence length of the two pairing states, whereby the one with the shorter coherence length is at an advantage. The coherence length can be defined as

[TABLE]

for both states separately. For dominant $p$-wave pairing, a ratio of $\xi_{p}/\xi_{f}\approx 0.24$ at $T_{0}$ thus implies that the $f$-wave component is more fragile. At $T_{c}$, however, the coherence length of the (total) order parameter diverges. Being closer to its pure $T_{c}$, the dominant channel has a longer coherence length and is thus more fragile. This behavior is corroborated by the opposite case observed for a stronger $f$-wave pairing interaction, where at $T_{c}$ disorder favors the weaker state while at $T_{0}$ the $p$-wave state is more quickly suppressed by impurities with $\xi_{p}/\xi_{f}\approx 1.25$. Note, finally, that all these findings apply also for combination with a relative phase of $\pi$ between the two components. In other words, the influence of impurity scattering is generally not a reduction in anisotropy, but rather a suppression of the component with larger coherence length.

IV Thermodynamic properties

We now turn to the (low-temperature) thermodynamic properties of the various states discussed above. For comparison to a truly nodal state, we additionally study a $d_{x^{2}-y^{2}}$-wave pairing state. While the clean-limit behavior at $T\ll T_{c}$ is entirely dominated by the presence or absence of nodes, we show in the following how disorder complicates the picture. Finally, as the real situation is generically not a pure $p$-wave or pure $f$-wave order parameter, we study the influence of the gap anisotropy given by changing the ratio $\Delta_{p}/\Delta_{f}$.

IV.1 Pure pairing channels

First, we study the three pairing states independently including disorder by means of the Greens function formalism introduced above and linear response theory Luttinger and Ward (1960); Langer (1962); Ambegaokar and Tewordt (1964); Luttinger (1964); Schrieffer (1971); Keller et al. (1988); Lee and Wen (1997). For this purpose, we calculate the specific heat, the superfluid density, and the thermal conductivity. A short overview over their derivations is provided in App. C. Our results are shown in Figs. 5 (a)-(c).

Figure 5 (a) shows the specific heat divided by the temperature, $C/T$. In the clean limit, we observe that both $f$- and $d$-wave exhibit a linear temperature dependence down to low temperatures. While the linear dependence goes all the way to $T=0$ for the $d$-wave case, the $f$-wave state with no true node has a downwards bending at low temperatures. In the same temperature range, the $p$-wave state clearly shows the exponential suppression of the specific heat at low temperature. For all three states, there is no residual value of $C/T$ at zero temperature due to the vanishing density of states in the clean system. This changes with disorder, which leads to a non-vanishing $C/T$ at $T=0$ all three states. Importantly, impurities lead to an even stronger qualitative resemblance of the $f$-wave and the $d$-wave state. Note that for the calculations the same chemical potential was used for all three pairing states.

Turning to the superfluid density, Fig. 5 (b), both spin-triplet states exhibit an exponential saturation of the superfluid density at low temperatures due to their fully gapped excitation spectrum. The nodal $d$-wave case, on the other hand, shows the expected linear temperature dependence in the clean limit. Once disorder is included, however, the $f$-wave and $d$-wave states can hardly be distinguished anymore. Both curves have approximatively the same slope and only a different residual value at zero temperature.

Finally, the thermal conductivity, shown in Fig. 5 (c), yields a similar picture; In the clean limit, both the nodal spin-singlet and the very anisotropic $f$-wave state look much alike on an overall level but have a clear qualitative difference at very low temperature, where the former vanishes linearly in temperature and the latter has an exponential suppressing indicating a tiny gap. The suppression is naturally more pronounced for the large gap $p$-wave state. Much like the above discussed specific heat and superfluid density, disorder renders the $d$- and $f$-wave states practically indistinguishable with only different residual values at $T=0$. Again, the $p$-wave state shows a much stronger suppression and lower residual conductivity.

The residual value of the thermal conductance, in particular its saturation for finite impurity density, is often used as a characterizing feature of symmetry-imposed nodes, such as in the case of a $d$-wave state. Indeed, we see in Fig. 6, that the residual value for the system with near nodes always vanishes in the limit of low impurity concentrations. However, for deep enough nodes, the concentration necessary to see this trend can be very low.

We can shed some light into the general trends shown in Figs. 5 and 6 by looking at the angle-resolved density of states. Figure 7 shows the radially-integrated density of states for a moderate impurity concentration of $c/c_{c}=0.18$ restricted to a quarter of the BZ, $\phi\in[0,\pi/2]$. For all three states, there is a clear, slightly softened gap with excitations above it. The $d$-wave and $f$-wave states have very similar density of states due to their nodal or near-nodal structure, with low-energy states induced by the impurity scattering well visible. The dominant feature appears around the nodal direction ($\phi=\pi/4$) with a peak in the density of low-energy states whose weight has been shifted from the higher-energy region. The region of the nodes has the highest Fermi velocity, making states in this region especially important for heat transport. It is the relocation of states that is crucially affecting the low-temperature behavior and the two pairing states’ resemblance observed in the thermodynamic properties. In the case of $p$-wave pairing, the impurities only give rise to subgap states at low energy, with a broad accumulation, where the Fermi surface approaches the van Hove points ($\phi=0,\pi/2$). The low-temperature behavior is thus still close to that of a fully gapped system with only little residual density of states at $T=0$. Finally, note again how Fig. 6 shows that the ”nodal” behavior of the $f$-wave state depends crucially on how close the gap zero is to the Fermi surface.

IV.2 Tuning the spin-triplet anisotropy strength

While we have restricted our discussion of thermodynamic properties to pure $p$- or $f$-wave states so far, we expect a superposition of the two in the general case. The effect of impurities in this case is then twofold: First, it changes the ratio of the two channels, thus changing the anisotropy. Second, the impurities lead to low-energy density of states, which depends on the anisotropy of the gap in the excitation spectrum. Therefore, we finish our discussion of thermodynamic properties with the case of different superpositions. In the following, we assume a relative phase of $\pi$ between the two components, such that the zeros of the gap are moved from $\mathbf{k}=(\pm\pi/2,\pm\pi/2)$ towards the BZ center and are always enclosed by the Fermi surface ($\mu=0.1$). Under these conditions the increase of $|\Delta_{p}|$ relative to $|\Delta_{f}|$ leads to a more isotropic gap, since the gap zero moves away from the FS.

In order to illustrate the effect of disorder on the superposition, we investigate the heat conductance at very low temperature and follow the residual value, $\kappa_{0}/T$ for $T\to 0$. In Fig. 8, we display $\kappa_{0}/T$ for different ratios $\Delta_{p}/\Delta_{f}$ at $c=0$ and compare it also with the nodal $d$-wave gap. We observe that for the $d$-wave state, $\kappa_{0}/T$ varies little with $c$. On the other hand, for all considered ratios $\Delta_{p}/\Delta_{f}$ the values of $\kappa_{0}/T$ are suppressed to zero for $c\to 0$ and approach each other with increasing $c$. While these residual heat conductance behave similar to the $d$-wave state for larger $c/c_{c}$, they always lie lower. As Fig. 6 shows the residual value increases the closer the gap zeros come to the FS, when we compare the $f$-wave states for different chemical potentials. From both Figs. 8 and 6 we see that the characteristic feature of downwards bending of $\kappa_{0}/T$ for $c$ going towards zero is gradually disappearing when the FS is removed from the gap zero. This behavior is characteristic for near nodal states.

V Discussion and Conclusion

Going beyond a pure point group analysis of gap-function zeros for chiral superconducting state, we have found a rich structure of gap singularities and associated winding numbers in the Brillouin zone. Depending on the exact form of the order parameter and the position of the Fermi surface, this structure influences the topological invariant of the system, the Chern number, and can lead to a pronounced anisotropy in the gap function as defined along the normal state Fermi surface. This anisotropy in the excitation spectrum has clear signatures in the thermodynamic properties of the system, exemplified by the specific heat, superfluid density, and thermal transport in Fig. 5. For a Fermi surface close to the gap-function zeros, the anisotropic chiral $f$-wave state becomes almost indistinguishable from a (nodal) $d$-wave state for even small disorder.

A generic order parameter is not a pure chiral $p$-wave or $f$-wave state, but a superposition of the two (and higher angular momentum channels). This raises two important issues that we have addressed in this work: First, is a gap-function zero close to the Fermi surface stable, or will the system find a superposition of $p$-wave and $f$-wave states to avoid this situation? Employing a Ginzburg-Landau analysis of this situation, we have indeed found that the zeros are generally not avoided. Instead, we find a topological transition when moving the Fermi surface through the zeros, which is of second order, see Fig. 3.

A second important issue concerns the influence of disorder on the gap anisotropy. Given the basic notion that disorder leads to a more isotropic order parameter, it is crucial to study the evolution of anisotropy in the presence of disorder for a superposition of two unconventional states. Interestingly, we find that the anisotropy can both increase and decrease when the impurity concentration is increased. We explain our unusual observation by first noting that both spin-triplet states are unconventional in nature and break time-reversal symmetry. Therefore, both gap functions are suppressed by disorder Anderson (1959). For substantially different gap contributions of the two states, the respective coherence lengths ($\xi\propto 1/|\Delta|^{2}$) differ as well and with it the two components’ susceptibility to disorder. The resulting stronger suppression of the pairing channel with longer coherence length can both increase and decrease the gap-function anisotropy.

To conclude, our work highlights the difficulty in interpreting thermodynamic measurements in case of chiral superconducting order parameters. Taken together with Wang et al. (2019); Dodaro et al. (2018), our findings show that thermodynamic measurements to distinguish a chiral $f$-wave state from a $d$-wave state are not necessarily conclusive even when considering different impurity concentrations. Such an analysis was done recently for Sr2RuO4 Hassinger et al. (2017) and motivated this work. However, there is a strong dependence of the anisotropy of an $f$-wave state on the Fermi surface position, in contrast to a $d$-wave state with point-group symmetry imposed nodes. Thus, Fermi surface engineering, as for example suggested by Hsu et al. (2016), combined with thermodynamic measurements could provide a path forward.

Acknowledgements

We are very grateful to Jose L. Lado, Jonathan Ruhman, and Catherine Kallin for many useful discussions. This work was financially supported by the Swiss National Science Foundation (SNSF) through Division II (No. 184739).

Appendix A Singularities of gap functions due to lattice translational symmetries

Using inversion in Eq. (1), we generally find for an odd-parity order parameter ($\phi_{\rm I}=\pi$)

[TABLE]

Together with Eq. (2), we identify

[TABLE]

and all equivalent points in the BZ as zeros of $\Delta_{\mathbf{k}}$.

Since the square lattice is bipartite, the NN (NNN) pairing states couple electrons on different (the same) sublattices. This leads to a distinguishing property between inter- and intra-sublattice pairing states, namely

[TABLE]

and

[TABLE]

with $\mathbf{Q}=(\pi,\pi)$. Note that this is a ‘reciprocal lattice vector’ considering a folded BZ. For intra-sublattice states, we can use the relations

[TABLE]

to identify, through $\mathbf{k}=-\mathbf{k}+\mathbf{Q}$, another set of points

[TABLE]

with $\Delta_{\mathbf{k}}^{\rm intra}=0$. $\Delta_{\mathbf{k}}^{\rm inter}$ has no further zeros. For the $p$-wave (inter-sublattice) and $f$-wave (intra-sublattice) pairing state, Eqs. (31) and (35) provide a complete list of zeros within the BZ. Note that for general odd-parity pairing states, more zeros can appear whose position require a more complex analysis. Obviously, all such zeros require that both components of the pairing function $\Phi_{x}(\mathbf{k})$ and $\Phi_{y}(\mathbf{k})$ vanish.

Appendix B Lattices with three-fold rotation axis

In a hexagonal lattice with full $D_{6h}$ symmetry Nandkishore et al. (2012); Kiesel et al. (2012), there are four two-dimensional representations, two even and two odd under inversion, allowing for both spin-singlet and spin-triplet chiral states. However, only one of each can be realized in the two-dimensional case. Using the primitive lattice vectors $\mathbf{T}_{1}=(1,0)$, and $\mathbf{T}_{2,3}=-(1/2,\mp\sqrt{3}/2)$, we can write them as

[TABLE]

with $\omega=\exp(2\pi i/3)$. These gap functions correspond, first, to a chiral $d$-wave and second, to a chiral $p$-wave order parameter. Due to the three-fold rotation, which maps $K$ ($K^{\prime}$) to itself, but adds a non-trivial phase to the gap function, the order parameter has to vanish at the $K$ and $K^{\prime}$ points for both order parameters. Additionally, the chiral $p$-wave order parameter has zeros on the $M$ points due to inversion symmetry, see Fig. 9.

While these two order parameters can not mix in the hexagonal lattice and the zeros at the high-symmetry points are fixed, the situation is different when the lattice symmetry is reduced to $D_{3h}$ by making the two sublattices in the hexagonal lattice distinct, a situation that is for example realized in the transition-metal dichalcogenides (TMDs) Ribak et al. (2019) or in a single layer of SrPtAs Nishikubo et al. (2011); Goryo et al. (2012); Fischer and Goryo (2015). Lacking inversion symmetry, the system allows for combinations of the spin-singlet and spin-triplet order parameters. In particular, angular momentum is only defined modulo 3, such that the $l=1$ pairing state can mix with the one with $l=-2$. Starting from the chiral $d$-wave situation with $l=-2$, the zero at the $\Gamma$ point splits up into one zero with winding number $1$ and three with winding number $-1$. These latter zeros move towards the $K$ points, see Fig. 9 333In this case, there are two non-degenerate bands with different relative phase in the order parameters, such that the other gap function has zeros moving towards the $K^{\prime}$ points. Depending on the Fermi surface topology, this can lead to several topological transitions until the zeros with $-1$ reach the $M$ points for a pure chiral $p$-wave state.

Appendix C Derivation of thermodynamic quantities

Here, we take a closer look at the derivation of the density of states, specific heat, superfluid density and thermal conductivity. We employ the Green’s function formalism and linear response theory Mineev et al. (1999); Langer (1962); Ambegaokar and Tewordt (1964); Luttinger (1964); Schrieffer (1971); Nomura (2005).

The density of states is in our framework Mineev et al. (1999) given by

[TABLE]

with $G^{\text{ret}}(\mathbf{k},E)$ being the retarded Green’s function, obtained by analytically continuing the frequencies to the real-axis, $i\omega_{n}\rightarrow E+i0$. In a dirty system, $E$ is replaced by $\tilde{E}=E-\Sigma^{\text{ret}}(E)$. The explicit expression of the Green’s functions is shown below.

Moving on to the specific heat, we start with the generalization Keller et al. (1988) of the ground-state energy formula by Luttinger and Ward Luttinger and Ward (1960) for an interacting electron system. The grand potential is given by

[TABLE]

with $i\tilde{\omega}_{n}=i\omega_{n}-\Sigma(i\omega_{n})$ and

[TABLE]

where $\Sigma(i\omega_{n})=cT(i\omega_{n})=\sum_{\nu}\Sigma_{\nu}(i\omega_{n})$. To make sure that the sum over $n$ converges rapidly, we consider the difference between superconducting and normal state, $\Omega_{s}-\Omega_{n}$.

Once the self-energy has been determined self-consistently, we calculate the specific heat difference through

[TABLE]

The specific heat in the normal state can now be added to the numerical solution of this equation to obtain $C_{s}/T$. Note that we have dropped the subscript in the main text.

Now let us take a look at the superfluid density $n_{s}$. In the clean system it is per definition Lee and Wen (1997) equal to

[TABLE]

Thus, it is enough to calculate $n_{n}$ by taking the real part of the current-current correlation function Schrieffer (1971),

[TABLE]

whose definition is

[TABLE]

It incorporates the electrical-flux operator,

[TABLE]

The $\mu$-component of the velocity, $v_{\mathbf{k}\mu}$, obeys the relation

[TABLE]

The final result then reads

[TABLE]

where, according to Gor’kov’s equations Mineev et al. (1999), the Green’s functions take the form

[TABLE]

Note that one has to use electron-hole symmetry, in other words a constant density of states, to obtain this specific structure of the Green’s functions. Even if the system may not have this symmetry, we assume it here for the sake of simplicity. Note that the results are not affected by this assumption.

Lastly, we consider the thermal conductivity $\kappa$. Again, linear response theory is applied. Our starting point is

[TABLE]

where $K^{\text{ret}}_{\mu\nu}(\mathbf{q},\Omega)$ is the retarded thermal-flux correlation function. Without the analytical continuation of the bosonic Matsubara frequencies, $i\Omega_{m}\rightarrow\Omega+i0$, and Fourier transformation, it is given by

[TABLE]

with $j_{\mu}^{\text{Th}}(\mathbf{q},\tau)$ the thermal-flux operator,

[TABLE]

After applying a standard mean-field approximation to Eq. (52) and reformulating everything in terms of Green’s functions, we make use of the spectral representation of the latter, perform the Matsubara frequency summation and subsequently take the limits. In the end Eq. (51) yields

[TABLE]

The retarded Green’s functions, which appear in this equation, are defined as

[TABLE]