# Spin susceptibility of superconductors with strong spin-orbit coupling

K. V. Samokhin111E-mail<EMAIL_ADDRESS>Department of Physics,
Brock University, St. Catharines, Ontario L2S 3A1, Canada

###### Abstract

We show that in some trigonal and hexagonal crystals the Zeeman coupling of
band electrons with an external magnetic field is strongly anisotropic and
necessarily vanishes along the main symmetry axis. This leads to qualitative
changes in the temperature dependence of the electron spin susceptibility in
the superconducting state. In particular, the power-law exponents at low
temperatures due to the contribution of the nodal quasiparticles are
considerably modified compared to their textbook values.

##  I Introduction

Measuring the spin magnetic response is an important characterization tool of
superconductors, both conventional and unconventional.SU-review ; TheBook  The
usual assumption is that the Zeeman interaction of the conduction electrons
with an external magnetic field has essentially the same form as for bare
electrons, i.e., is independent of the wave vector. Then, the temperature
behaviour of the spin susceptibility depends only on the superconducting gap
structure, but not on the character of the electron Bloch bands. If the
superconducting pairing occurs in the spin-singlet channel, then the
susceptibility tensor $\chi_{ij}(T)=\chi(T)\delta_{ij}$ is entirely determined
by the thermally excited quasiparticles and in a fully gapped superconducting
state one would see an exponentially decreasing $\chi(T)$ at low temperatures.
In contrast, if the gap has zeros on the Fermi surface, then $\chi(T)\propto
T^{2}$ for isolated first-order point nodes, or $\chi(T)\propto T$ for line
nodes or second-order point nodes. In the triplet case, the susceptibility
acquires a nontrivial tensor structure depending on the mutual orientation of
the order parameter $\bm{d}$ and the external magnetic field $\bm{H}$. For
$\bm{d}\parallel\bm{H}$ the spin response is determined only by the
excitations and the susceptibility has the same temperature dependence as in
the singlet case, whereas for $\bm{d}\perp\bm{H}$ both the Cooper pairs and
the excitations contribute to the spin response and the susceptibility is
unchanged compared to the normal state.

The electron-lattice spin-orbit (SO) coupling in many superconductors of
current interest, in particular, the heavy-fermion compounds, is strong
compared with the Zeeman energy and the energy scales associated with
superconductivity, and one can call into question the validity of the
assumption that the effective Zeeman interaction in the Bloch bands has the
same simple form $\hat{H}_{Z}=\mu_{B}\bm{H}\hat{\bm{\sigma}}$ as for bare
electrons (here the spin magnetic moment is equal to
$-\mu_{B}\hat{\bm{\sigma}}$, $\mu_{B}$ is the Bohr magneton, the electron
charge is $-e$, the Landé factor is set to $2$, and $\hat{\bm{\sigma}}$ are
the Pauli matrices). Even in the presence of the SO coupling, the Bloch bands
in a centrosymmetric crystal are twofold degenerate.Kittel-book  Focusing on
just one band, one can write the effective Zeeman Hamiltonian at the wave
vector $\bm{k}$ in the general form as
$\hat{H}_{Z}=-\hat{\bm{m}}(\bm{k})\bm{H}$, where the Hermitian $2\times 2$
matrix $\hat{\bm{m}}$ has the meaning of the magnetic moment of the band
electrons.CB60 ; Roth66 ; LL-9

While the effects of the Zeeman coupling anisotropy due to the SO coupling
have been extensively studied in the context of semiconductors, see, e.g.,
Ref. Winkler-book, , its consequences for the spin response of superconductors
have received relatively little attention in the literature. These
consequences are expected to be profound, especially in the light of the
recent developments in the field,WW09 ; j32 ; Fis13 ; RS16 ; NHI16 ; Sam19 ;
Sam19-PRB  which demonstrated that the SO coupling can significantly affect
the symmetry of the Bloch bands and invalidate the pseudospin pictureAnd84 ;
UR85  commonly used in the theory of unconventional superconductivity.

The pseudospin approach is based on the assumption that the twofold degenerate
Bloch states transform under the crystal point group operations in the same
way as the pure spin states $\lvert\uparrow\rangle$ and
$\lvert\downarrow\rangle$ at each wave vector, including the center of the
Brillouin zone (the $\Gamma$ point). However, the latter cannot always be
true, because, mathematically, not all double-valued irreducible
representations of the crystal point group (or, more accurately,
corepresentations of the magnetic point group) are equivalent to the spin-1/2
representation. We will show how a “non-pseudospin” character of the Bloch
states leads to a strongly anisotropic Zeeman interaction in some bands in
trigonal and hexagonal crystals.

Our main goal is to calculate the spin susceptibility of superconductors with
an anisotropic Zeeman coupling, whose form is constrained only by the symmetry
of the system. Impurities are neglected, as are the effects of the interaction
of the electron charges with the magnetic field. In Secs. II and III, we
develop the general symmetry-based theory of the spin response of the Bloch
bands in crystals with the SO coupling and identify the cases in which the
Zeeman interaction is qualitatively different from the usual isotropic
expression. In Sec. IV, the spin response in the superconducting state is
studied and the spin susceptibility of a general multiband superconductor is
derived. In Sec. V, we discuss applications to singlet and triplet
superconductors of hexagonal symmetry. Some technical details are relegated to
the Appendices. Throughout the paper we use the units in which $\hbar=1$ and
$k_{B}=1$, neglecting, in particular, the difference between the quasiparticle
momentum and the wave vector.

##  II Zeeman coupling in the band representation

We consider a centrosymmetric nonmagnetic crystal with $N$ twofold degenerate
bands crossing the Fermi level. The Bloch states are degenerate at each wave
vector $\bm{k}$ due to the conjugation symmetry ${\cal C}=KI$ (Ref. Kittel-
book, ): the states $|\bm{k},n,1\rangle$ and $|\bm{k},n,2\rangle\equiv{\cal
C}|\bm{k},n,1\rangle$ are orthogonal and have the same energy. Here $K$ is the
time reversal (TR) operation, $I$ is the space inversion, $n=1,...,N$ labels
the bands, and the index $s=1,2$, which distinguishes two degenerate states
within the same band, is called the conjugation index. The conjugation
operation is antiunitary and, since the Bloch states are spin-$1/2$ spinors,
we have ${\cal C}|\bm{k},n,2\rangle=-|\bm{k},n,1\rangle$, i.e., ${\cal
C}^{2}=-1$. The conjugation degeneracy, which is also called the Kramers
degeneracy, can be lifted by an external magnetic field. Note that the
assumption of inversion symmetry is important: in a non-centrosymmetric
crystal the band degeneracy is lifted by the SO coupling even at zero magnetic
field, with wide-reaching implications for superconductivity.NCSC-book

The Zeeman coupling of electrons with the magnetic field $\bm{H}$ is described
by the Hamiltonian $\hat{H}_{Z}=-\bm{H}\hat{\bm{m}}$, where
$(\hat{m}_{x},\hat{m}_{y},\hat{m}_{z})=-\mu_{B}(\hat{\sigma}_{1},\hat{\sigma}_{2},\hat{\sigma}_{3})$
is the operator of the spin magnetic moment. In the band representation, the
Zeeman Hamiltonian is diagonal in the wave vector, but not in the band and
conjugation indices:

$\langle\bm{k},n,s|\hat{H}_{Z}|\bm{k},n^{\prime},s^{\prime}\rangle=-\bm{H}\bm{m}_{nn^{\prime},ss^{\prime}}(\bm{k}).$
(1)

The matrix elements of the spin magnetic moment in the Bloch basis can be
represented as linear combinations of the unit matrix and the Pauli matrices
in the conjugacy space as follows:

$\displaystyle m_{nn^{\prime},ss^{\prime},j}(\bm{k})$ $\displaystyle\equiv$
$\displaystyle\langle\bm{k},n,s|\hat{m}_{j}|\bm{k},n^{\prime},s^{\prime}\rangle$
(2) $\displaystyle=$ $\displaystyle
iA_{nn^{\prime},j}(\bm{k})\delta_{ss^{\prime}}+\bm{B}_{nn^{\prime},j}(\bm{k})\bm{\sigma}_{ss^{\prime}}.$

The intraband and interband components correspond to $n=n^{\prime}$ and $n\neq
n^{\prime}$, respectively.

In principle, the functions $A$ and $\bm{B}$ in Eq. (2) can be calculated if a
treatable model of the band structure is available. Here we adopt a different
approach: without resorting to any particular microscopic model, we regard $A$
and $\bm{B}$ as phenomenological parameters, whose form is only constrained by
the symmetries of the system. In particular, from the requirement that
$\hat{\bm{m}}$ is a Hermitian operator, we have
$\langle\psi|\hat{\bm{m}}|\psi^{\prime}\rangle=\langle\psi^{\prime}|\hat{\bm{m}}|\psi\rangle^{*}$,
therefore,

$A_{nn^{\prime},i}(\bm{k})=-A_{n^{\prime}n,i}^{*}(\bm{k}),\quad\bm{B}_{nn^{\prime},i}(\bm{k})=\bm{B}_{n^{\prime}n,i}^{*}(\bm{k}).$
(3)

Additional constraints are obtained by analyzing the response of the spin
magnetic moment to TR and the point group operations, which in turn depends on
the transformation properties of the Bloch states.

From the group-theoretical point of view, the spinor Bloch states
$|\bm{k},n,1\rangle$ and $|\bm{k},n,2\rangle$ form the basis of an irreducible
double-valued two-dimensional (2D) corepresentation (corep) of the magnetic
point group of the wave vector $\bm{k}$, see Appendix A. Four-dimensional
coreps, which are possible in crystals of cubic symmetry, are not considered
here. If the corep at the $\Gamma$ point is described by $2\times 2$ matrices
$\hat{\cal D}_{n}(g)$, where $g$ is an element of the crystal point group
$\mathbb{G}$, then one can construct the Bloch bases at $\bm{k}\neq\bm{0}$
using the following expression:Sam19 ; Sam19-PRB

$g|\bm{k},n,s\rangle=\sum_{s^{\prime}}|g\bm{k},n,s^{\prime}\rangle{\cal
D}_{n,s^{\prime}s}(g),$ (4)

which defines the basis at the wave vector $g\bm{k}$ given the basis at
$\bm{k}$.

In a “pseudospin” band, the $\Gamma$-point corep is equivalent to the
spin-$1/2$ corep and the conjugate Bloch states affected by the SO coupling
respond to the point group operations in exactly the same way as the pure spin
states. Therefore, one can put $\hat{\cal D}_{n}(g)=\hat{D}^{(1/2)}(R)$ for
both proper rotations ($g=R$) and improper rotations ($g=IR$). Here

$\hat{D}^{(1/2)}(R)=e^{-i\theta(\bm{n}\hat{\bm{\sigma}})/2}$ (5)

is the spin-1/2 representation of a counterclockwise rotation $R$ through an
angle $\theta$ about an axis $\bm{n}$. Thus, in a pseudospin band Eq. (4)
reproduces the well-known Ueda-Rice formula.UR85

In general, however, the Bloch states at the $\Gamma$ point correspond to a
corep which is not equivalent to the spin-$1/2$ corep. Therefore, $\hat{\cal
D}_{n}(g)\neq\hat{D}^{(1/2)}(R)$ for some $g$ and the Ueda-Rice prescription
does not work. In this case, the band is called a “non-pseudospin” band. The
corep matrices in twofold degenerate non-pseudospin bands are shown in Table
1, see Appendix A for details.

Table 1: The $\Gamma$-point corepresentation matrices in twofold degenerate non-pseudospin bands; $R$ denotes the rotational generators of the point group $\mathbb{G}$ and the spin rotation matrix $\hat{D}^{(1/2)}(R)$ is given by Eq. (5). The notations for the point group elements are the same as in Ref. LL-3, . $\mathbb{G}$ | corep | $\hat{\cal D}(R)$
---|---|---
$\mathbf{C}_{4h}$ | $(\Gamma_{7},\Gamma_{8})$ | $\hat{\cal D}(C_{4z})=-\hat{D}^{(1/2)}(C_{4z})$
$\mathbf{D}_{4h}$ | $\Gamma_{7}$ |  $\hat{\cal D}(C_{4z})=-\hat{D}^{(1/2)}(C_{4z})$, $\hat{\cal D}(C_{2y})=\hat{D}^{(1/2)}(C_{2y})$
$\mathbf{C}_{3i}$ | $\Gamma_{6}$ | $\hat{\cal D}(C_{3z})=-\hat{\sigma}_{0}$
$\mathbf{D}_{3d}$ | $(\Gamma_{5},\Gamma_{6})$ |  $\hat{\cal D}(C_{3z})=-\hat{\sigma}_{0}$, $\hat{\cal D}(C_{2y})=\hat{D}^{(1/2)}(C_{2y})$
$\mathbf{C}_{6h}$ | $(\Gamma_{9},\Gamma_{10})$ | $\hat{\cal D}(C_{6z})=-\hat{D}^{(1/2)}(C_{6z})$
| $(\Gamma_{11},\Gamma_{12})$ | $\hat{\cal D}(C_{6z})=-\hat{D}^{(1/2)}(C_{2z})$
$\mathbf{D}_{6h}$ | $\Gamma_{8}$ |  $\hat{\cal D}(C_{6z})=-\hat{D}^{(1/2)}(C_{6z})$, $\hat{\cal D}(C_{2y})=\hat{D}^{(1/2)}(C_{2y})$
| $\Gamma_{9}$ |  $\hat{\cal D}(C_{6z})=-\hat{D}^{(1/2)}(C_{2z})$, $\hat{\cal D}(C_{2y})=\hat{D}^{(1/2)}(C_{2y})$
$\mathbf{O}_{h}$ | $\Gamma_{7}$ |  $\hat{\cal D}(C_{4z})=-\hat{D}^{(1/2)}(C_{4z})$, $\hat{\cal D}(C_{2y})=\hat{D}^{(1/2)}(C_{2y})$, $\hat{\cal D}(C_{3xyz})=\hat{D}^{(1/2)}(C_{3xyz})$

It follows from Eq. (4) that under the space inversion the conjugate Bloch
states transform as $I|\bm{k},n,s\rangle=p_{n}|-\bm{k},n,s\rangle$, where
$p_{n}=\pm$ denotes the parity of the band. The pseudovector quantity
$\hat{\bm{m}}$ must remain invariant under inversion,
$I\hat{\bm{m}}I^{-1}=\hat{\bm{m}}$, therefore

$\langle\bm{k},n,s|\hat{\bm{m}}|\bm{k},n^{\prime},s^{\prime}\rangle=p_{n}p_{n^{\prime}}\langle-\bm{k},n,s|\hat{\bm{m}}|-\bm{k},n^{\prime},s^{\prime}\rangle.$

Using the representation (2), we find that the parity of $A$ and $\bm{B}$ is
determined by the relative parity of the bands $p_{n}p_{n^{\prime}}$. For the
transformation under the TR operation $K={\cal C}I$, we have

$\begin{array}[]{l}K|\bm{k},n,1\rangle=p_{n}|-\bm{k},n,2\rangle,\vskip 6.0pt
plus 2.0pt minus 2.0pt\\\
K|\bm{k},n,2\rangle=-p_{n}|-\bm{k},n,1\rangle.\end{array}$ (6)

Being a magnetic moment operator, $\hat{\bm{m}}$ must change sign under TR:
$K\hat{\bm{m}}K^{-1}=-\hat{\bm{m}}$, therefore

$\displaystyle\langle\bm{k},n,s|\hat{\bm{m}}|\bm{k},n^{\prime},s^{\prime}\rangle=-\langle\bm{k},n,s|K^{\dagger}\hat{\bm{m}}K|\bm{k},n^{\prime},s^{\prime}\rangle$
$\displaystyle=-p_{n}p_{n^{\prime}}\sum_{s_{1}s_{2}}\sigma_{2,s^{\prime}s_{1}}\langle-\bm{k},n^{\prime},s_{1}|\hat{\bm{m}}|-\bm{k},n,s_{2}\rangle\sigma_{2,s_{2}s}.$

In the second line here we used Eq. (6) and the expression
$\langle\psi|K^{\dagger}|\psi^{\prime}\rangle=\langle\psi^{\prime}|K|\psi\rangle$,
which follows from the antiunitarity of $K$. In terms of $A$ and $\bm{B}$, we
obtain:
$A_{nn^{\prime},i}(\bm{k})=-p_{n}p_{n^{\prime}}A_{n^{\prime}n,i}(-\bm{k})$ and
$\bm{B}_{nn^{\prime},i}(\bm{k})=p_{n}p_{n^{\prime}}\bm{B}_{n^{\prime}n,i}(-\bm{k})$.
Putting all together, we find that the functions $A$ and $\bm{B}$ are real and
satisfy

$\begin{array}[]{c}A_{nn^{\prime},i}(\bm{k})=-A_{n^{\prime}n,i}(\bm{k}),\vskip
6.0pt plus 2.0pt minus 2.0pt\\\
A_{nn^{\prime},i}(-\bm{k})=p_{n}p_{n^{\prime}}A_{nn^{\prime},i}(\bm{k}),\end{array}$
(7)

and

$\begin{array}[]{c}\bm{B}_{nn^{\prime},i}(\bm{k})=\bm{B}_{n^{\prime}n,i}(\bm{k}),\vskip
6.0pt plus 2.0pt minus 2.0pt\\\
\bm{B}_{nn^{\prime},i}(-\bm{k})=p_{n}p_{n^{\prime}}\bm{B}_{nn^{\prime},i}(\bm{k}).\end{array}$
(8)

Additional constraints are imposed by the requirement of invariance under the
remaining elements of the crystal point group.

Under a point-group operation $g$ (e.g., a rotation), an arbitrary state of
the system is transformed into another state as
$|\psi\rangle\to|\tilde{\psi}\rangle=g|\psi\rangle$ and the magnetic field is
transformed into $g\bm{H}$. The matrix elements of the Zeeman Hamiltonian
$\hat{H}_{Z}(\bm{H})$ in the basis $|\psi\rangle$ are equal to the matrix
elements of the Hamiltonian $\hat{H}_{Z}(g\bm{H})$ in the transformed basis
$|\tilde{\psi}\rangle$. Therefore,

$\langle\bm{k},n,s|(\hat{\bm{m}}\cdot\bm{H})|\bm{k},n^{\prime},s^{\prime}\rangle=\langle\bm{k},n,s|g^{\dagger}(\hat{\bm{m}}\cdot
g\bm{H})g|\bm{k},n^{\prime},s^{\prime}\rangle,$

see, e.g., Ref. Bir-Pikus-book, . Using Eqs. (2) and (4), the invariance
condition under $g\in\mathbb{G}$ takes the following form:

$\displaystyle
iA_{nn^{\prime},j}(\bm{k})\hat{\sigma}_{0}+\bm{B}_{nn^{\prime},j}(\bm{k})\hat{\bm{\sigma}}$
$\displaystyle\quad=\sum_{k=1}^{3}R_{jk}(g)\bigl{[}iA_{nn^{\prime},k}(g^{-1}\bm{k})\hat{\tau}_{nn^{\prime},0}(g)$
$\displaystyle\quad+\bm{B}_{nn^{\prime},k}(g^{-1}\bm{k})\hat{\bm{\tau}}_{nn^{\prime}}(g)\bigr{]},$
(9)

where $\hat{R}$ is the $3\times 3$ rotation matrix and

$\displaystyle\hat{\tau}_{nn^{\prime},0}(g)=\hat{\cal D}_{n}(g)\hat{\cal
D}^{\dagger}_{n^{\prime}}(g),$
$\displaystyle\hat{\tau}_{nn^{\prime},\nu}(g)=\hat{\cal
D}_{n}(g)\hat{\sigma}_{\nu}\hat{\cal
D}^{\dagger}_{n^{\prime}}(g),\quad\nu=1,2,3.$

We used the fact that any element of the point group is either a proper
rotation $g=R$ or an improper rotation $g=IR$, and that the magnetic field is
a pseudovector and not affected by inversion, therefore $g\bm{H}=R(g)\bm{H}$.
The most general phenomenological expressions for the functions $A$ and
$\bm{B}$ in a given pair of bands $n$ and $n^{\prime}$ can be found by solving
Eq. (II) for all generators of $\mathbb{G}$, taking into account the
constraints (7) and (8).

##  III Intraband Zeeman coupling

We shall see below, in Sec. IV, that the intraband and interband components of
the Zeeman coupling produce qualitatively different contributions to the spin
magnetic response. Namely, the temperature dependence of the spin
susceptibility in the superconducting state is almost entirely determined by
the intraband coupling. Putting $n=n^{\prime}$ in Eq. (7), we find
$A_{nn,i}=0$. Introducing the notation $\bm{B}_{nn,i}=-\bm{\mu}_{n,i}$, the
spin magnetic moment in the $n$th band can be written as

$\hat{m}_{nn,i}(\bm{k})=-\bm{\mu}_{n,i}(\bm{k})\hat{\bm{\sigma}},$ (10)

so that the intraband Zeeman coupling takes the following form:

$\langle\bm{k},n,s|\hat{H}_{Z}|\bm{k},n,s^{\prime}\rangle=\sum_{i,\nu}\mu_{n,i\nu}(\bm{k})H_{i}\sigma_{\nu,ss^{\prime}},$
(11)

where $i=x,y,z$ and $\nu=1,2,3$. The coefficients $\mu_{n,i\nu}(\bm{k})$ form
a $3\times 3$ matrix, denoted as $\hat{\mu}_{n}(\bm{k})$, which generalizes
the Bohr magneton $\mu_{B}$ in the $n$th band. This matrix is real and even in
$\bm{k}$, according to Eq. (8), but not necessarily symmetric, as we shall see
below.

From the requirement of invariance under the point group operations, see Eq.
(II), we obtain:

$\hat{\mu}_{n}(\bm{k})=\hat{R}(g)\hat{\mu}_{n}(g^{-1}\bm{k})\hat{\cal
R}^{-1}_{n}(g),\quad g\in\mathbb{G},$ (12)

where the $3\times 3$ orthogonal matrix $\hat{\cal R}_{n}$ is defined by the
following expression:

$\hat{\cal D}_{n}^{\dagger}(g)\hat{\sigma}_{\mu}\hat{\cal
D}_{n}(g)=\sum_{\nu=1}^{3}{\cal R}_{n,\mu\nu}(g)\hat{\sigma}_{\nu}.$ (13)

This matrix depends on the symmetry of the $n$th band, namely the
$\Gamma$-point corep, but not on the band parity. We see that, while the
Zeeman coupling matrix $\hat{\mu}_{n}(\bm{k})$ is invariant under the point
group operations, in the sense given by Eq. (12), it does not transform like a
tensor, in general.

In all pseudospin bands, we have $\hat{\cal D}_{n}(g)=\hat{D}^{(1/2)}(R)$ and,
using the well-known identity

$\hat{D}^{(1/2),\dagger}(R)\hat{\sigma}_{\mu}\hat{D}^{(1/2)}(R)=\sum_{\nu}R_{\mu\nu}\hat{\sigma}_{\nu},$
(14)

we find that $\hat{\cal R}_{n}(g)=\hat{R}(g)$ for all elements of the point
group. Therefore, Eq. (12) becomes

$\hat{\mu}_{n}(\bm{k})=\hat{R}(g)\hat{\mu}_{n}(g^{-1}\bm{k})\hat{R}^{-1}(g).$
(15)

It is easy to see that the simplest solution of this last equation is given by

$\mu_{n,i\nu}(\bm{k})=\mu_{n}\delta_{i\nu}.$ (16)

Thus, in pseudospin bands one can use the “usual”, i.e., isotropic and
momentum-independent, expression for the intraband Zeeman coupling. Note that
$\mu_{n}$ includes the effects of the crystal field and the SO coupling and
therefore does not have to be equal to the Bohr magneton $\mu_{B}$. In
uniaxial crystals, one can use, if needed, a more general solution of Eq.
(15): $\mu_{n,x1}=\mu_{n,y2}=\mu_{n,\perp}$ and $\mu_{n,z3}=\mu_{n,z}$, which
distinguishes between the basal-plane and $z$-axis Zeeman couplings, but is
still $\bm{k}$-independent.

In non-pseudospin bands, using the $\Gamma$-point corep matrices from Table 1,
we obtain that in most cases $\hat{\cal R}_{n}(g)=\hat{R}(g)$ and the
expression (16) still works. However, in certain bands in trigonal and
hexagonal crystals we have $\hat{\cal R}_{n}(g)\neq\hat{R}(g)$ for some
elements of the point group, so that the invariance equation for
$\hat{\mu}_{n}(\bm{k})$ does not have the form (15). Namely, this happens in
the following four exceptional cases:

$\begin{array}[]{l}\Gamma_{6}\ \mathrm{of}\
\mathbf{C}_{3i},\quad(\Gamma_{5},\Gamma_{6})\ \mathrm{of}\
\mathbf{D}_{3d},\vskip 6.0pt plus 2.0pt minus 2.0pt\\\
(\Gamma_{11},\Gamma_{12})\ \mathrm{of}\ \mathbf{C}_{6h},\quad\Gamma_{9}\
\mathrm{of}\ \mathbf{D}_{6h}.\end{array}$ (17)

Note that these are the same bands in which the usual classification of
triplet superconducting states fails, see Ref. Sam19-PRB, . Below we show that
the intraband Zeeman coupling in the exceptional cases is very anisotropic and
that some components of $\hat{\mu}_{n}(\bm{k})$ necessarily vanish, for
symmetry reasons, along the main symmetry axis. The band index $n$ will be
dropped, for brevity.

###  III.1 $\Gamma_{6}$ of $\mathbf{C}_{3i}$

The group $\mathbf{C}_{3i}$ is generated by the rotations $C_{3z}$ about the
vertical ($z$) axis and by the inversion $I$. According to Table 1, the
$\Gamma$-point corep matrix is given by $\hat{\cal
D}_{\Gamma}(C_{3z})=-\hat{\sigma}_{0}$, therefore $\hat{\cal
R}(C_{3z})=\hat{\mathbb{1}}$ (the $3\times 3$ unit matrix) and the invariance
condition (12) takes the following form:

$\hat{\mu}(\bm{k})=\hat{R}(C_{3z})\hat{\mu}(C_{3z}^{-1}\bm{k}).$ (18)

Here the rotation matrix is given by

$\hat{R}(C_{3z})=\left(\begin{array}[]{ccc}-1/2&-\sqrt{3}/2&0\\\
\sqrt{3}/2&-1/2&0\\\ 0&0&1\end{array}\right)$

and the rotations act on the wave vector as $C_{3z}^{-1}k_{\pm}=e^{\mp 2\pi
i/3}k_{\pm}$ and $C_{3z}^{-1}k_{z}=k_{z}$, where $k_{\pm}=k_{x}\pm ik_{y}$.

The solution of Eq. (18) has a block-diagonal form:

$\hat{\mu}=\left(\begin{array}[]{ccc}\mu_{x1}&\mu_{x2}&0\\\
\mu_{y1}&\mu_{y2}&0\\\ 0&0&\mu_{z3}\end{array}\right).$ (19)

While the equation $\mu_{z3}(\bm{k})=\mu_{z3}(C_{3z}^{-1}\bm{k})$ is trivially
satisfied by a constant $\mu_{z3}(\bm{k})=\mu_{z}$, the solutions for the
other components are more involved. Introducing $q^{\pm}_{1}=\mu_{x1}\pm
i\mu_{y1}$ and $q^{\pm}_{2}=\mu_{x2}\pm i\mu_{y2}$, we obtain from Eq. (18):

$\begin{array}[]{l}q^{\pm}_{1}(\bm{k})=e^{\pm 2\pi
i/3}q^{\pm}_{1}(C_{3z}^{-1}\bm{k}),\vskip 6.0pt plus 2.0pt minus 2.0pt\\\
q^{\pm}_{2}(\bm{k})=e^{\pm 2\pi
i/3}q^{\pm}_{2}(C_{3z}^{-1}\bm{k}),\end{array}$ (20)

One can easily see that $q^{\pm}_{1,2}(0,0,k_{z})=0$, therefore

$\mu_{x1}=\mu_{x2}=\mu_{y1}=\mu_{y2}=0\qquad\mathrm{at}\ k_{x}=k_{y}=0.$ (21)

Thus, the basal-plane components of the Zeeman coupling matrix necessarily
depend on $\bm{k}$ and vanish along the threefold symmetry axis. Using the
lowest-order polynomial solutions of Eq. (20), we have

$\hat{\mu}(\bm{k})=\left(\begin{array}[]{ccc}(\alpha_{1}k_{+}+\alpha_{1}^{*}k_{-})k_{z}&(\alpha_{2}k_{+}+\alpha_{2}^{*}k_{-})k_{z}&0\\\
(-i\alpha_{1}k_{+}+i\alpha_{1}^{*}k_{-})k_{z}&(-i\alpha_{2}k_{+}+i\alpha_{2}^{*}k_{-})k_{z}&0\\\
0&0&\mu_{z}\end{array}\right),$ (22)

where $\alpha_{1,2}$ are complex constants and $\mu_{z}$ is a real constant.
Note that the vanishing of the basal-plane components of $\hat{\mu}$ at
$k_{z}=0$ is accidental, since one can consider a more general solution, e.g.,
$\mu_{x1}=(\alpha_{1}k_{+}+\alpha_{1}^{*}k_{-})k_{z}+(\alpha_{3}k_{+}^{2}+\alpha_{3}^{*}k_{-}^{2})$,
which has zeros only along the main axis, as dictated by symmetry.

###  III.2 $(\Gamma_{5},\Gamma_{6})$ of $\mathbf{D}_{3d}$

The group $\mathbf{D}_{3d}$ is generated by the rotations $C_{3z}$ and
$C_{2y}$, and also by $I$. From Table 1 we obtain: $\hat{\cal
R}(C_{3z})=\hat{\mathbb{1}}$ and $\hat{\cal R}(C_{2y})=\hat{R}(C_{2y})$. The
first of the two invariance conditions is the same as Eq. (18), while the
second one has the form

$\hat{\mu}(\bm{k})=\hat{R}(C_{2y})\hat{\mu}(C_{2y}^{-1}\bm{k})\hat{R}^{-1}(C_{2y}),$
(23)

where

$\hat{R}(C_{2y})=\left(\begin{array}[]{ccc}-1&0&0\\\ 0&1&0\\\
0&0&-1\end{array}\right)$

and $C_{2y}^{-1}(k_{+},k_{-},k_{z})=(-k_{-},-k_{+},-k_{z})$.

The solution for the Zeeman coupling matrix has a block-diagonal form (19),
with $\mu_{z3}(\bm{k})=\mu_{z}$. The condition (23) produces two additional
constraints: $q^{\pm}_{1}(\bm{k})=q^{\mp}_{1}(C_{2y}^{-1}\bm{k})$ and
$q^{\pm}_{2}(\bm{k})=-q^{\mp}_{2}(C_{2y}^{-1}\bm{k})$. Imposing them on Eq.
(22), we obtain:

$\hat{\mu}(\bm{k})=\left(\begin{array}[]{ccc}\beta_{1}k_{x}k_{z}&-\beta_{2}k_{y}k_{z}&0\\\
\beta_{1}k_{y}k_{z}&\beta_{2}k_{x}k_{z}&0\\\ 0&0&\mu_{z}\end{array}\right),$
(24)

where $\beta_{1,2}$ and $\mu_{z}$ are real constants. The zeros of the basal-
plane components of $\hat{\mu}$ at $k_{x}=k_{y}=0$ are required by symmetry,
whereas the zero at $k_{z}=0$ is accidental and can be removed by taking into
account additional terms in the solution, e.g.,
$\mu_{x1}=\beta_{1}k_{x}k_{z}+\beta_{3}(k_{x}^{2}-k_{y}^{2})$, with a real
$\beta_{3}$.

###  III.3 $(\Gamma_{11},\Gamma_{12})$ of $\mathbf{C}_{6h}$

The group $\mathbf{C}_{6h}$ is generated by $C_{6z}$ and $I$. From Table 1 we
obtain: $\hat{\cal R}(C_{6z})=\hat{R}(C_{2z})$, therefore the invariance
condition (12) takes the following form:

$\hat{\mu}(\bm{k})=\hat{R}(C_{6z})\hat{\mu}(C_{6z}^{-1}\bm{k})\hat{R}^{-1}(C_{2z}).$
(25)

Here

$\displaystyle\hat{R}(C_{6z})=\left(\begin{array}[]{ccc}1/2&-\sqrt{3}/2&0\\\
\sqrt{3}/2&1/2&0\\\ 0&0&1\end{array}\right),$
$\displaystyle\hat{R}(C_{2z})=\left(\begin{array}[]{ccc}-1&0&0\\\ 0&-1&0\\\
0&0&1\end{array}\right),$

while the action of rotations on the wave vector is given by
$C_{6z}^{-1}k_{\pm}=e^{\mp\pi i/3}k_{\pm}$ and $C_{6z}^{-1}k_{z}=k_{z}$.

The Zeeman coupling matrix has a block-diagonal form (19), with
$\mu_{z3}(\bm{k})=\mu_{z}$. The basal-plane components satisfy the equations
$q^{\pm}_{1}(\bm{k})=e^{\mp 2\pi i/3}q^{\pm}_{1}(C_{6z}^{-1}\bm{k})$,
$q^{\pm}_{2}(\bm{k})=e^{\mp 2\pi i/3}q^{\pm}_{2}(C_{6z}^{-1}\bm{k})$, and
necessarily vanish along the sixfold symmetry axis. For the lowest-order
polynomial solution we obtain:

$\hat{\mu}(\bm{k})=\left(\begin{array}[]{ccc}\alpha_{1}k_{+}^{2}+\alpha_{1}^{*}k_{-}^{2}&\alpha_{2}k_{+}^{2}+\alpha_{2}^{*}k_{-}^{2}&0\\\
i\alpha_{1}k_{+}^{2}-i\alpha_{1}^{*}k_{-}^{2}&i\alpha_{2}k_{+}^{2}-i\alpha_{2}^{*}k_{-}^{2}&0\\\
0&0&\mu_{z}\end{array}\right),$ (28)

where $\alpha_{1,2}$ are complex constants and $\mu_{z}$ is a real constant.

###  III.4 $\Gamma_{9}$ of $\mathbf{D}_{6h}$

The group $\mathbf{D}_{6h}$ is generated by the rotations $C_{6z}$ and
$C_{2y}$, and by $I$. It follows from Table 1 that $\hat{\cal
R}(C_{6z})=\hat{R}(C_{2z})$ and $\hat{\cal R}(C_{2y})=\hat{R}(C_{2y})$,
therefore the invariance conditions are given by Eqs. (25) and (23). As
before, the solution for the Zeeman coupling matrix has a block-diagonal form
(19), where $\mu_{z3}$ is $\bm{k}$-independent, but the basal-plane components
necessarily depend on $\bm{k}$.

Imposing the constraint (23) on Eq. (28), we obtain:

$\hat{\mu}(\bm{k})=\left(\begin{array}[]{ccc}\beta_{1}(k_{x}^{2}-k_{y}^{2})&2\beta_{2}k_{x}k_{y}&0\\\
-2\beta_{1}k_{x}k_{y}&\beta_{2}(k_{x}^{2}-k_{y}^{2})&0\\\
0&0&\mu_{z}\end{array}\right),$ (29)

where $\beta_{1,2}$ and $\mu_{z}$ are real constants. The vanishing of the
basal-plane components along the sixfold axis is required by symmetry.

###  III.5 Summary

We have shown that the intraband Zeeman coupling in the four exceptional non-
pseudospin bands (17) has a nontrivial matrix structure and momentum
dependence, see Eqs. (22), (24), (28), and (29). The components of the spin
magnetic moment operator in the band representation, Eq. (10), have the
following form:

$\displaystyle\hat{m}_{x,y}(\bm{k})=\left(\begin{array}[]{cc}0&\Phi_{x,y}(\bm{k})\\\
\Phi_{x,y}^{*}(\bm{k})&0\end{array}\right),$ (32)
$\displaystyle\hat{m}_{z}(\bm{k})=\left(\begin{array}[]{cc}-\mu_{z}&0\\\
0&\mu_{z}\end{array}\right).$ (35)

The basal-plane components $\hat{m}_{x,y}$ vanish along the main symmetry
axis, with the functions $\Phi_{x}$ and $\Phi_{y}$ being either linear in
$k_{x},k_{y}$ (in trigonal crystals), or quadratic in $k_{x},k_{y}$ (in
hexagonal crystals).

##  IV Spin susceptibility in the superconducting state

Let us now introduce the superconducting pairing of quasiparticles in the
Bloch bands constructed using Eq. (4). Suppose there are $N$ bands crossing
the Fermi level and participating in superconductivity, then the intraband
pairing affected by the Zeeman interaction with an external magnetic field can
be described by the following Bardeen-Cooper-Schrieffer (BCS) mean field
Hamiltonian:

$\hat{H}=\hat{H}_{0}+\hat{H}_{Z}+\hat{H}_{sc}.$ (36)

The three terms here are as follows:

$\hat{H}_{0}=\sum_{\bm{k},n,s}\xi_{n}(\bm{k})c^{\dagger}_{\bm{k},n,s}c_{\bm{k},n,s},$
(37)

where $\xi_{n}(\bm{k})=\xi_{n}(-\bm{k})$ are the band dispersions counted from
the chemical potential,

$\hat{H}_{Z}=-\sum_{\bm{k},nn^{\prime},ss^{\prime}}\bm{H}\bm{m}_{nn^{\prime},ss^{\prime}}(\bm{k})c^{\dagger}_{\bm{k},n,s}c_{\bm{k},n^{\prime},s^{\prime}},$
(38)

where the matrix elements of the spin magnetic moment are given by Eq. (2),
and

$\hat{H}_{sc}=\frac{1}{2}\sum_{\bm{k},n,ss^{\prime}}\left[\Delta_{n,ss^{\prime}}(\bm{k})c^{\dagger}_{\bm{k},n,s}\tilde{c}^{\dagger}_{\bm{k},n,s^{\prime}}+\mathrm{H.c.}\right]$
(39)

describes the pairing between time-reversed statesAnd59 ; Blount85  in the
same band, with

$\tilde{c}^{\dagger}_{\bm{k},n,s}=Kc^{\dagger}_{\bm{k},n,s}K^{-1}=p_{n}\sum_{s^{\prime}}c^{\dagger}_{-\bm{k},n,s^{\prime}}(-i\hat{\sigma}_{2})_{s^{\prime}s},$
(40)

according to Eq. (6). In the weak-coupling picture, the gap functions are
nonzero only inside the BCS shells near the Fermi surfaces:
$\hat{\Delta}_{n}(\bm{k})\propto\theta(\epsilon_{c}-|\xi_{n}(\bm{k})|)$, where
$\epsilon_{c}$ is the energy cutoff.

By analogy with the standard theory of unconventional superconductivity,SU-
review ; TheBook  one can consider separately the pairing in the singlet and
triplet channels, keeping in mind that the Kramers index $s$ is neither spin
nor pseudospin, in general. The matrix structure of the gap functions in the
Kramers space can be represented in the following form:

$\hat{\Delta}_{n}(\bm{k})=\psi_{n}(\bm{k})\hat{\sigma}_{0},\quad\psi_{n}(-\bm{k})=\psi_{n}(\bm{k}),$
(41)

in the singlet channel and

$\hat{\Delta}_{n}(\bm{k})=\bm{d}_{n}(\bm{k})\hat{\bm{\sigma}},\quad\bm{d}_{n}(-\bm{k})=-\bm{d}_{n}(\bm{k}),$
(42)

in the triplet channel. Note that, since we defined the gap function in Eq.
(39) as a measure of the pairing between the time-reversed states
$|\bm{k},n,s\rangle$ and $K|\bm{k},n,s^{\prime}\rangle$, instead of between
$|\bm{k},n,s\rangle$ and $|-\bm{k},n,s^{\prime}\rangle$, the expressions (41)
and (42) do not contain the factors $i\hat{\sigma}_{2}$, in contrast to the
convention used in Refs. SU-review,  and TheBook, .

According to Eq. (4), the transformation rule for the electron creation
operators in the Bloch states is given by
$gc^{\dagger}_{\bm{k},n,s}g^{-1}=\sum_{s^{\prime}}c^{\dagger}_{g\bm{k},n,s^{\prime}}{\cal
D}_{n,s^{\prime}s}(g)$. Using the fact that the TR operation $K$ commutes with
all elements of the point group and is antilinear, we obtain that the point-
group rotations and reflections in the physical space induce the following
transformation of the gap function matrix in the momentum space:Sam19-PRB

$g:\hat{\Delta}_{n}(\bm{k})\to\hat{\cal
D}_{n}(g)\hat{\Delta}_{n}(g^{-1}\bm{k})\hat{\cal D}^{\dagger}_{n}(g).$

In terms of the singlet and triplet components, this becomes

$g:\psi_{n}(\bm{k})\to\psi_{n}(g^{-1}\bm{k}),\ \bm{d}_{n}(\bm{k})\to{\cal
R}_{n}(g)\bm{d}_{n}(g^{-1}\bm{k}),$ (43)

where the $3\times 3$ orthogonal matrix $\hat{\cal R}_{n}$ is defined by Eq.
(13).

It follows from Eq. (43) that the singlet gap function always transforms under
the point group as a complex scalar, regardless of the band symmetry,
therefore the usual classification of the singlet superconducting statesSU-
review ; TheBook  is applicable. In contrast, the transformation properties of
the triplet gap function essentially depend on the band symmetry at the
$\Gamma$ point. In a pseudospin band, we have $\hat{\cal R}_{n}(g)=\hat{R}$
and $\bm{d}_{n}(\bm{k})$ transforms under $g$ into
$R\bm{d}_{n}(g^{-1}\bm{k})$, i.e., like a pseudovector. However, this is not
the case in general: if the $\Gamma$-point corep is such that $\hat{\cal
R}_{n}(g)\neq\hat{R}$, which happens in the exceptional bands in trigonal and
hexagonal crystals, see Eq. (17), then $\bm{d}_{n}(\bm{k})$ does not transform
under $g$ like a pseudovector, with profound consequences for its momentum
dependence and the nodal structure.Sam19-PRB

###  IV.1 Magnetic response

Writing the Zeeman Hamiltonian (38) in the form
$\hat{H}_{Z}=-\bm{H}\hat{\cal{\bm{M}}}$, where $\hat{\cal{\bm{M}}}$ is the
operator of the total spin magnetic moment of electrons, we define the
magnetization as $\bm{M}={\cal V}^{-1}\langle\hat{\cal{\bm{M}}}\rangle$ (the
angular brackets denote the thermodynamic average and ${\cal V}$ is the volume
of the system). Introducing the normal and anomalous Green’s functions, see
Appendix B, and taking the thermodynamic limit ${\cal V}\to\infty$, we obtain:

$\bm{M}=T\sum_{m}\int\frac{d^{3}\bm{k}}{(2\pi)^{3}}\sum_{nn^{\prime}}\,\mathrm{tr}\,[\hat{\bm{m}}_{nn^{\prime}}(\bm{k})\hat{G}_{n^{\prime}n}(\bm{k},\omega_{m})],$
(44)

where $\omega_{m}=(2m+1)\pi T$ is the fermionic Matsubara frequency and
“$\,\mathrm{tr}\,$” denotes the $2\times 2$ matrix trace with respect to the
conjugation indices. We would like to note that the non-pseudospin character
of the electron bands does not affect the way $\bm{M}$ transforms under the
symmetry operations. If the physical state of the system is acted upon by an
operation $g\in\mathbb{G}$, then $\langle\hat{\cal{\bm{M}}}\rangle$ is
transformed into $\langle g^{-1}\hat{\cal{\bm{M}}}g\rangle$. Using Eq. (II),
it is easy to show that $g^{-1}\hat{\cal M}_{i}g=\sum_{j}R_{ij}(g)\hat{\cal
M}_{j}$, where $R(g)$ is the rotational part of $g$, and therefore $\bm{M}\to
R(g)\bm{M}$, as expected.

In a weak field, expanding the Green’s function in Eq. (44) in powers of
$\bm{H}$, we find $\bm{M}=\bm{M}_{0}+\bm{M}_{1}+{\cal O}(H^{2})$. The first
term here is the spontaneous magnetization given by

$M_{0,i}=-T\sum_{m}\int\frac{d^{3}\bm{k}}{(2\pi)^{3}}\sum_{n}\bm{\mu}_{n,i}(\bm{k})\,\mathrm{tr}\,[\hat{\bm{\sigma}}\hat{G}_{n}(\bm{k},\omega_{m})],$
(45)

according to Eq. (10), whereas the second term describes the linear magnetic
response and can be written as $M_{1,i}=\sum_{j}\chi^{(spin)}_{ij}H_{j}$,
where

$\displaystyle\chi^{(spin)}_{ij}=-T\sum_{m}\int\frac{d^{3}\bm{k}}{(2\pi)^{3}}\sum_{nn^{\prime}}\,\mathrm{tr}\,\bigl{[}\hat{m}_{nn^{\prime},i}(\bm{k})\hat{G}_{n^{\prime}}(\bm{k},\omega_{m})\hat{m}_{n^{\prime}n,j}(\bm{k})\hat{G}_{n}(\bm{k},\omega_{m})$
$\displaystyle+\hat{m}_{nn^{\prime},i}(\bm{k})\hat{F}_{n^{\prime}}(\bm{k},\omega_{m})\hat{m}_{n^{\prime}n,j}(\bm{k})\hat{F}^{\dagger}_{n}(\bm{k},\omega_{m})\bigr{]}$
(46)

is the spin susceptibility tensor. The Green’s functions in Eqs. (45) and
(IV.1) correspond to zero field and are band-diagonal, see Appendix B. The
expressions (45) and (IV.1) are the multiband generalizations of the standard
formulas for the spontaneous magnetization and the spin susceptibility in
unconventional superconductors.

In this paper, we consider only unitary superconducting states, which satisfy
$\hat{\Delta}_{n}\hat{\Delta}_{n}^{\dagger}\propto\hat{\sigma}_{0}$. According
to Eqs. (41) and (42), any singlet state is unitary, whereas a triplet state
is unitary only if $\bm{d}_{n}\times\bm{d}_{n}^{*}=\bm{0}$. The Green’s
functions have the following form:

$\begin{array}[]{l}\hat{G}_{n}(\bm{k},\omega_{m})=-\dfrac{i\omega_{m}+\xi_{n}(\bm{k})}{\omega_{m}^{2}+E_{n}^{2}(\bm{k})}\hat{\sigma}_{0},\vskip
6.0pt plus 2.0pt minus 2.0pt\\\
\hat{F}_{n}(\bm{k},\omega_{m})=\dfrac{\hat{\Delta}_{n}(\bm{k})}{\omega_{m}^{2}+E_{n}^{2}(\bm{k})},\end{array}$
(47)

where

$E_{n}(\bm{k})=\sqrt{\xi_{n}^{2}(\bm{k})+|\Delta_{n}(\bm{k})|^{2}}$ (48)

is the energy of quasiparticle excitations. The energy gap is given by
$|\Delta_{n}|=|\psi_{n}|$ in the singlet case and $|\Delta_{n}|=|\bm{d}_{n}|$
in the unitary triplet case. The zeros of the energy (the “gap nodes”) in the
$n$th band correspond to the intersections of the Fermi surface, defined by
$\xi_{n}(\bm{k})=0$, with the manifold of zeros of $\psi_{n}(\bm{k})$ or
$\bm{d}_{n}(\bm{k})$.

We obtain from Eq. (45) that $\bm{M}_{0}=\bm{0}$, i.e., there is no
spontaneous spin magnetization in a unitary state, regardless of the symmetry
of the electron bands involved in the pairing. The spin susceptibility tensor
(IV.1) can be represented in the following form:

$\chi^{(spin)}_{ij}=\tilde{\chi}_{ij}+\sum_{n}\chi_{n,ij},$ (49)

where the first term contains the interband contributions (corresponding to
$n\neq n^{\prime}$) and the second term is the sum of the intraband
susceptibilities:

$\displaystyle\chi_{n,ij}=-T\sum_{m}\int\frac{d^{3}\bm{k}}{(2\pi)^{3}}\,\mathrm{tr}\,\bigl{[}(\bm{\mu}_{n,i}\hat{\bm{\sigma}})\hat{G}_{n}(\bm{\mu}_{n,j}\hat{\bm{\sigma}})\hat{G}_{n}$
$\displaystyle+(\bm{\mu}_{n,i}\hat{\bm{\sigma}})\hat{F}_{n}(\bm{\mu}_{n,j}\hat{\bm{\sigma}})\hat{F}^{\dagger}_{n}\bigr{]}.$
(50)

Below we show that the temperature dependence of the susceptibility in the
superconducting state is determined by the intraband terms, while the
interband contribution is essentially unaffected by the superconducting
transition.

Let us start with the normal (“$N$”) state. Using the fact that the intraband
Zeeman coupling $\hat{\mu}_{n}$ is real, we obtain from Eq. (IV.1) the
following expression:

$\displaystyle\chi^{N}_{n,ij}=2\int\frac{d^{3}\bm{k}}{(2\pi)^{3}}\left(-\frac{\partial
f}{\partial\xi_{n}}\right)\bigl{[}\hat{\mu}_{n}(\bm{k})\hat{\mu}_{n}^{\dagger}(\bm{k})\bigr{]}_{ij}$
$\displaystyle=2N_{F,n}\left\langle\bigl{[}\hat{\mu}_{n}(\bm{k})\hat{\mu}_{n}^{\dagger}(\bm{k})\bigr{]}_{ij}\right\rangle_{FS,n},$
(51)

where $f(\epsilon)=1/(e^{\epsilon/T}+1)$ is the Fermi function, $N_{F,n}$ is
the Fermi-level density of states in the $n$th band, and the angular brackets
denote the Fermi-surface average. Thus, the intraband susceptibility is
determined by the quasiparticles near the Fermi surface and is essentially
temperature independent. The “usual” expression is recovered if one neglects
the anisotropy of the Zeeman coupling: putting
$\mu_{n,i\nu}(\bm{k})=\mu_{B}\delta_{i\nu}$, we obtain:
$\chi^{N}_{n,ij}=\chi_{P}\delta_{ij}$, where $\chi_{P}=2N_{F,n}\mu_{B}^{2}$ is
the Pauli susceptibility.

Regarding the interband susceptibility, we focus on just one pair of bands
($n,n^{\prime}=1,2$) with $\xi_{2}(\bm{k})>\xi_{1}(\bm{k})$, so that the band
splitting is given by ${\cal E}(\bm{k})=\xi_{2}(\bm{k})-\xi_{1}(\bm{k})$.
Introducing the notations $A_{12,i}=-A_{21,i}=a_{i}$ and
$\bm{B}_{12,i}=\bm{B}_{21,i}=\bm{b}_{i}$, see Eqs. (7) and (8), the interband
components of the magnetic moment take the following form:

$\begin{array}[]{l}\hat{m}_{12,i}(\bm{k})=ia_{i}(\bm{k})+\bm{b}_{i}(\bm{k})\hat{\bm{\sigma}},\vskip
6.0pt plus 2.0pt minus 2.0pt\\\
\hat{m}_{21,i}(\bm{k})=-ia_{i}(\bm{k})+\bm{b}_{i}(\bm{k})\hat{\bm{\sigma}}.\end{array}$
(52)

Here $a_{i}$ and $\bm{b}_{i}$ are real functions of $\bm{k}$, whose parity is
determined by the relative parity of the bands. From Eq. (IV.1), we obtain:

$\tilde{\chi}^{N}_{ij}=4\int\frac{d^{3}\bm{k}}{(2\pi)^{3}}(a_{i}a_{j}+\bm{b}_{i}\bm{b}_{j})\frac{f(\xi_{1})-f(\xi_{2})}{\xi_{2}-\xi_{1}}.$
(53)

We see that the interband susceptibility is determined by all quasiparticles
in the momentum-space shell “sandwiched” between the two Fermi surfaces and
its temperature dependence is negligible at temperatures small compared to the
band splitting.

One can also expect that the interband susceptibility is not affected by the
superconducting transition. Indeed, the difference between its values in the
superconducting (“$S$”) and normal states is very small: even at $T=0$, it can
be estimated as
$\tilde{\chi}^{S}_{ij}-\tilde{\chi}^{N}_{ij}\propto(|\Delta|^{2}/{\cal
E}^{2})\ln(\epsilon_{c}/|\Delta|)$, see Appendix C. In weak-coupling
superconductors, the band splitting far exceeds the energy scales associated
with superconductivity: $|\Delta_{1,2}|\ll\epsilon_{c}\ll{\cal E}$, therefore
one can put $\tilde{\chi}^{S}_{ij}=\tilde{\chi}^{N}_{ij}$.

Summarizing our findings so far, the total spin susceptibility (49) can be
written in the form

$\chi^{(spin)}_{ij}(T)=\tilde{\chi}_{ij}^{N}+\sum_{n}\chi_{n,ij}(T),$ (54)

where the first term is the temperature-independent interband contribution
(53). The second term is given by Eq. (IV.1) and can be calculated in each
band individually.

###  IV.2 Intraband susceptibility

The normal and anomalous Green’s functions are given by Eq. (47). In a singlet
state, dropping the band index $n$ and summing over the Matsubara frequencies
in Eq. (IV.1), we obtain:

$\chi_{ij}(T)=2N_{F}\left\langle\bigl{[}\hat{\mu}(\bm{k})\hat{\mu}^{\dagger}(\bm{k})\bigr{]}_{ij}{\cal
Y}(\bm{k},T)\right\rangle_{FS},$ (55)

where

${\cal
Y}(\bm{k},T)=\frac{1}{2T}\int_{0}^{\infty}\frac{d\xi}{\cosh^{2}(\sqrt{\xi^{2}+|\Delta(\bm{k})|^{2}}/2T)}$
(56)

is the angle-resolved Yosida function (Ref. TheBook, ) and
$|\Delta(\bm{k})|=|\psi(\bm{k})|$ is the energy gap. Being a measure of the
number of the thermally excited Bogoliubov quasiparticles in the
superconducting state, the Yosida function vanishes at $T\to 0$. In the normal
state, we have ${\cal Y}=1$ and Eq. (IV.1) is recovered.

In the triplet case, the calculation is somewhat more involved. Introducing
the parallel and transverse projectors onto the triplet gap function
$\bm{d}(\bm{k})$,

$\Pi_{\parallel,\mu\nu}=\frac{\,\mathrm{Re}\,(d_{\mu}d_{\nu}^{*})}{|\bm{d}|^{2}},\quad\Pi_{\perp,\mu\nu}=\delta_{\mu\nu}-\frac{\,\mathrm{Re}\,(d_{\mu}d_{\nu}^{*})}{|\bm{d}|^{2}},$

we obtain:

$\displaystyle\chi_{ij}(T)=2N_{F}\left\langle\bigl{[}\hat{\mu}(\bm{k})\hat{\Pi}_{\parallel}(\bm{k})\hat{\mu}^{\dagger}(\bm{k})\bigr{]}_{ij}{\cal
Y}(\bm{k},T)\right\rangle_{FS}$
$\displaystyle+2N_{F}\left\langle\bigl{[}\hat{\mu}(\bm{k})\hat{\Pi}_{\perp}(\bm{k})\hat{\mu}^{\dagger}(\bm{k})\bigr{]}_{ij}\right\rangle_{FS},$
(57)

see Appendix D for details. Here ${\cal Y}$ is the Yosida function (56) with
$|\Delta(\bm{k})|=|\bm{d}(\bm{k})|$.

As explained in Sec. III, in most bands one can use the isotropic Zeeman
coupling of the form $\mu_{i\nu}(\bm{k})=\mu_{B}\delta_{i\nu}$, see Eq. (16)
with the Bohr magneton used for simplicity. This leads to the standard
expressions for the spin susceptibility. Namely, in the singlet case Eq. (55)
yields $\chi_{ij}(T)=\chi(T)\delta_{ij}$, with

$\displaystyle\chi(T)$ $\displaystyle=$
$\displaystyle\chi_{P}\left\langle{\cal Y}(\bm{k},T)\right\rangle_{FS}$ (58)
$\displaystyle=$
$\displaystyle\chi_{P}\frac{1}{2T}\int_{0}^{\infty}\frac{\nu(E)\,dE}{\cosh^{2}(E/2T)}.$

We introduced the dimensionless density of states (DoS) of the Bogoliubov
quasiparticles:

$\nu(E)=\left\langle\frac{E}{\sqrt{E^{2}-|\Delta(\bm{k})|^{2}}}\right\rangle_{FS},$
(59)

where the angular integration is performed over the directions satisfying
$|\Delta(\bm{k})|\leq E$. Thus, at low temperatures the susceptibility is
entirely determined by the thermally excited quasiparticles near the gap nodes
and has the following asymptotics:SU-review ; TheBook

$\chi(T)\propto\left\\{\begin{array}[]{ll}e^{-\Delta/T},&\mathrm{fully\
gapped},\\\ T^{2/k},&k^{\mathrm{th}}\ \mathrm{order\ point\ nodes},\\\
T,&\mathrm{line\ nodes}.\end{array}\right.$ (60)

Here $\Delta$ denotes the value of the gap in an isotropic superconducting
state or the minimum gap in a nodeless anisotropic state.

In the triplet case with the isotropic Zeeman coupling
$\mu_{i\nu}(\bm{k})=\mu_{B}\delta_{i\nu}$, we obtain from Eq. (IV.2):

$\displaystyle\chi_{ij}(T)=\chi_{P}\left\langle\Pi_{\parallel,ij}(\bm{k}){\cal
Y}(\bm{k},T)\right\rangle_{FS}$
$\displaystyle+\chi_{P}\left\langle\Pi_{\perp,ij}(\bm{k})\right\rangle_{FS},$
(61)

i.e., the magnetic response depends on the relative orientation of $\bm{d}$
and $\bm{H}$ (Ref. Leggett75, ). If $\bm{d}(\bm{k})\parallel\bm{H}$, then the
susceptibility is determined only by the excitations, see the first term in
Eq. (IV.2), and vanishes at $T\to 0$. If $\bm{d}(\bm{k})\perp\bm{H}$, then
both the Cooper pairs and the excitations contribute to the susceptibility,
see the second term in Eq. (IV.2), leading to a nonzero residual
susceptibility at $T=0$.

In the exceptional bands, the intraband Zeeman coupling has a nontrivial
matrix structure and a complicated momentum dependence, see Sec. III.
Therefore, Eqs. (58) and (IV.2) are no longer applicable. In particular, the
quasiparticle contribution to the susceptibility is determined not only by the
gap nodal structure through the quasiparticle DoS (59), but also by the zeros
of the Zeeman coupling. Namely, if the nodal excitations are weakly coupled
with the external magnetic field, then they do not contribute to the
magnetization, which leads to significant changes in the temperature
dependence of $\chi_{ij}$.

##  V Application to hexagonal superconductors

As an application of the general theory developed above, in this section we
calculate the spin susceptibility of a hexagonal superconductor described by
the point group $\mathbf{D}_{6h}=\mathbf{D}_{6}\times\\{E,I\\}$. This group
has six double-valued coreps of even or odd parity,Lax-book  corresponding to
the $\Gamma^{\pm}_{7}$, $\Gamma^{\pm}_{8}$, or $\Gamma^{\pm}_{9}$ bands, with
only the $\Gamma^{+}_{7}$ bands being pseudospin ones. The parity superscript
can be omitted, because neither the intraband pairing nor the intraband Zeeman
coupling are affected by the band parity. We assume that there is just one
isotropic band of $\Gamma_{7}$, $\Gamma_{8}$, or $\Gamma_{9}$ symmetry
crossing the Fermi level and participating in superconductivity.

According to the Landau theory of phase transitions, the superconducting gap
functions are classified according to single-valued irreducible
representations (irreps) $\gamma$ of the point group $\mathbb{G}$, called the
pairing channels. In a $d$-dimensional pairing channel, the singlet gap
function can be represented in the form
$\psi(\bm{k})=\sum_{a=1}^{d}\eta_{a}\varphi_{a}(\bm{k})$, where $a$ labels the
scalar basis functions $\varphi$ of an even irrep $\gamma$. Similarly, the
triplet gap function can be represented as
$\bm{d}(\bm{k})=\sum_{a=1}^{d}\eta_{a}\bm{\varphi}_{a}(\bm{k})$, where the
basis functions $\bm{\varphi}$ of an odd irrep $\gamma$ have different form in
the pseudospin and non-pseudospin bands, see Ref. Sam19-PRB, . In both singlet
and triplet cases, the components of the superconducting order parameter
$\eta_{1},..,\eta_{d}$ are found by minimizing the free energy of the
superconductor.TheBook ; SU-review  The point group $\mathbf{D}_{6}$ has four
one-dimensional (1D) irreps $A_{1}$, $A_{2}$, $B_{1}$, $B_{2}$, and two 2D
irreps $E_{1}$, $E_{2}$ (we use the “chemical” notationLL-3  for the pairing
channels, reserving the $\Gamma$ notationLax-book  for the double-valued
coreps describing the symmetry of the Bloch bands). Taking into account the
parity, there are six singlet and six triplet pairing channels.

###  V.1 Singlet pairing

Examples of the basis functions for the even irreps of $\mathbf{D}_{6h}$ are
shown in Table 2. The basis functions in all channels except $A_{1g}$ (which
corresponds to the trivial, or “$s$-wave”, pairing) have zeros imposed by
symmetry. The $A_{2g}$ gap has twelve vertical lines of nodes, whereas the
$B_{1g}$ and $B_{2g}$ gaps have six vertical lines of nodes as well as a
horizontal line of nodes at $k_{z}=0$. For the 2D pairing channels, the nodal
structure depends on the state, i.e., on the order parameter components
$(\eta_{1},\eta_{2})$. The stable states correspond to
$(\eta_{1},\eta_{2})\propto(1,0)$ or $(1,1)$ (Ref. SU-review, ).

Table 2: Examples of the basis functions $\varphi(\bm{k})$ in the even pairing channels for $\mathbb{G}=\mathbf{D}_{6h}$ ($a$ is a real constant and $k_{\pm}=k_{x}\pm ik_{y}$). $\gamma$ |  $\Gamma_{7}$, $\Gamma_{8}$, and $\Gamma_{9}$ bands
---|---
$A_{1g}$ | $k_{x}^{2}+k_{y}^{2}+ak_{z}^{2}$
$A_{2g}$ | $i(k_{+}^{6}-k_{-}^{6})$
$B_{1g}$ | $(k_{+}^{3}+k_{-}^{3})k_{z}$
$B_{2g}$ | $i(k_{+}^{3}-k_{-}^{3})k_{z}$
$E_{1g}$ |  $k_{+}k_{z}$, $k_{-}k_{z}$
$E_{2g}$ |  $k_{+}^{2}$, $k_{-}^{2}$

####  V.1.1 $\Gamma_{7}$ and $\Gamma_{8}$ bands

In the non-exceptional $\Gamma_{7}$ and $\Gamma_{8}$ bands, one can use the
usual expression for the Zeeman coupling,
$\mu_{i\nu}(\bm{k})=\mu_{B}\delta_{i\nu}$. The susceptibility is determined by
the quasiparticle DoS, see Eq. (58), with the standard low-temperature
asymptotics (60).

####  V.1.2 $\Gamma_{9}$ bands

Substituting the Zeeman coupling (29) into Eq. (55) and using the fact that
the Yosida function has the same symmetry as the energy gap, one can show that
the off-diagonal elements of the susceptibility tensor vanish after the Fermi-
surface averaging, whereas the diagonal elements have the following form:

$\displaystyle\chi_{xx}=2N_{F}\bigl{\langle}\bigl{[}\beta_{1}^{2}(k_{x}^{2}-k_{y}^{2})^{2}+4\beta_{2}^{2}k_{x}^{2}k_{y}^{2}\bigr{]}$
$\displaystyle\hskip 119.50148pt\times{\cal Y}(\bm{k},T)\bigr{\rangle}_{FS},$
$\displaystyle\chi_{yy}=2N_{F}\bigl{\langle}\bigl{[}4\beta_{1}^{2}k_{x}^{2}k_{y}^{2}+\beta_{2}^{2}(k_{x}^{2}-k_{y}^{2})^{2}\bigr{]}$
(62) $\displaystyle\hskip 119.50148pt\times{\cal
Y}(\bm{k},T)\bigr{\rangle}_{FS},$
$\displaystyle\chi_{zz}=2N_{F}\mu_{z}^{2}\left\langle{\cal
Y}(\bm{k},T)\right\rangle_{FS}.$

In the normal state, we have ${\cal Y}=1$ and the susceptibility components
are given by

$\begin{array}[]{l}\chi_{xx}^{N}=\chi_{yy}^{N}=\dfrac{8}{15}N_{F}k_{F}^{4}(\beta_{1}^{2}+\beta_{2}^{2})\equiv\chi_{\perp}^{N},\vskip
6.0pt plus 2.0pt minus 2.0pt\\\ \chi_{zz}^{N}=2N_{F}\mu_{z}^{2},\end{array}$
(63)

assuming a spherical Fermi surface of radius $k_{F}$. In all superconducting
states, the temperature dependence of $\chi_{zz}$ is entirely determined by
the quasiparticle DoS, see Eq. (58), and has the asymptotics (60).

In contrast, the basal-plane components $\chi_{xx}$ and $\chi_{yy}$ contain
the Yosida function multiplied by strongly anisotropic factors, which vanish
along the sixfold symmetry axis, therefore the contribution of the
quasiparticles near the main axis to the Fermi-surface averages is suppressed.
The effect of this suppression is most pronounced in the superconducting
states which have only isolated point nodes at $k_{x}=k_{y}=0$. Using Table 2,
it is easy to see that there is only one such singlet state, namely, the
“chiral” $d$-wave state of the form $\psi(\bm{k})\propto k_{+}^{2}$ (or
$k_{-}^{2}$), which corresponds to the $E_{2g}$ irrep and has isolated second-
order point nodes along the main axis. In this state, the expressions (V.1.2)
yield

$\chi_{xx}(T)=\chi_{yy}(T)\propto T^{3},\quad\chi_{zz}(T)\propto T,$ (64)

at $T\to 0$. In all other singlet states, $\chi_{xx}$ and $\chi_{yy}$ have the
asymptotics (60).

###  V.2 Triplet pairing

In the triplet case, the gap function $\bm{d}(\bm{k})$ and the basis functions
$\bm{\varphi}(\bm{k})$ can be represented as linear combinations of three
orthonormal pseudovectors $\bm{e}_{1}$, $\bm{e}_{2}$, and $\bm{e}_{3}$, whose
transformations induced by the point group operations follow from Eq. (43),
namely, $\bm{e}_{\mu}\to{\cal R}(g)\bm{e}_{\mu}=\sum_{\nu}\bm{e}_{\nu}{\cal
R}_{\nu\mu}(g)$. Note that $\bm{e}_{1}$, $\bm{e}_{2}$, $\bm{e}_{3}$ are not
the same as the Cartesian basis vectors $\hat{\bm{x}}$, $\hat{\bm{y}}$,
$\hat{\bm{z}}$ in the physical space, which is easily seen from their
different response to inversion. However, they are not entirely independent:
if the spin quantization axis is chosen along $\hat{\bm{z}}$, then all corep
matrices for the rotations about $\hat{\bm{z}}$, see Table 1, commute with
$\hat{\sigma}_{3}$ and we obtain from Eq. (13) that $\bm{e}_{3}$ is unchanged
by these rotations. It is in this sense that $\bm{e}_{3}$ is “parallel” to
$\hat{\bm{z}}$.

The triplet basis functions for the odd pairing channels are shown in Table 3.
In the $\Gamma_{7}$ and $\Gamma_{8}$ bands we have $\hat{\cal R}(g)=\hat{R}$
for all $g$, therefore $\bm{d}$ and $\bm{e}_{\mu}$ transform under rotations
like vectors and the standard expressions for the basis functionsSU-review ;
TheBook  are applicable. In contrast, in the $\Gamma_{9}$ bands we have
$\hat{\cal R}(g)\neq\hat{R}$ for some $g$, therefore $\bm{d}$ and
$\bm{e}_{\mu}$ do not transform under all rotations like vectors, which
changes the momentum dependence of the basis functions.Sam19-PRB  Note that
all three components of $\bm{d}$ never vanish simultaneously in a whole plane
in the momentum space, so that Blount’s theorem about the absence of line
nodes in a generic triplet stateBlount85  holds in the exceptional non-
pseudospin bands as well.

The precise form of the basis functions in a given material, in particular,
the relative values of the coefficients in Table 3, is dictated by the
microscopic details. We will consider two cases:
$\bm{d}\parallel\hat{\bm{z}}$, which corresponds to $a_{1}=a_{2}=0$ in all
basis functions, and $\bm{d}\perp\hat{\bm{z}}$. In both cases, the triplet
states have isolated point nodes and/or accidental lines of nodes.

Table 3: Examples of the basis functions $\bm{\varphi}(\bm{k})$ in the odd pairing channels for $\mathbb{G}=\mathbf{D}_{6h}$, where $a_{1,2,3}$ are real constants, $k_{\pm}=k_{x}\pm ik_{y}$, and $\bm{e}_{\pm}=(\bm{e}_{1}\pm i\bm{e}_{2})/\sqrt{2}$. Second column: the bands in which $\bm{d}$ transforms under rotations like a vector. Third column: the bands in which $\bm{d}$ does not transform like a vector. $\gamma$ |  $\Gamma_{7}$ and $\Gamma_{8}$ bands |  $\Gamma_{9}$ bands
---|---|---
$A_{1u}$ | $a_{1}(k_{x}\bm{e}_{1}+k_{y}\bm{e}_{2})+a_{3}k_{z}\bm{e}_{3}$ | $a_{1}(k_{+}^{3}+k_{-}^{3})\bm{e}_{1}+ia_{2}(k_{+}^{3}-k_{-}^{3})\bm{e}_{2}+a_{3}k_{z}\bm{e}_{3}$
$A_{2u}$ | $a_{1}(k_{y}\bm{e}_{1}-k_{x}\bm{e}_{2})+ia_{3}(k_{+}^{6}-k_{-}^{6})k_{z}\bm{e}_{3}$ | $ia_{1}(k_{+}^{3}-k_{-}^{3})\bm{e}_{1}+a_{2}(k_{+}^{3}+k_{-}^{3})\bm{e}_{2}+ia_{3}(k_{+}^{6}-k_{-}^{6})k_{z}\bm{e}_{3}$
$B_{1u}$ | $a_{1}(k_{+}^{2}k_{z}\bm{e}_{+}+k_{-}^{2}k_{z}\bm{e}_{-})+a_{3}(k_{+}^{3}+k_{-}^{3})\bm{e}_{3}$ | $a_{1}k_{z}\bm{e}_{1}+ia_{2}(k_{+}^{6}-k_{-}^{6})k_{z}\bm{e}_{2}+a_{3}(k_{+}^{3}+k_{-}^{3})\bm{e}_{3}$
$B_{2u}$ | $ia_{1}(k_{+}^{2}k_{z}\bm{e}_{+}-k_{-}^{2}k_{z}\bm{e}_{-})+ia_{3}(k_{+}^{3}-k_{-}^{3})\bm{e}_{3}$ | $ia_{1}(k_{+}^{6}-k_{-}^{6})k_{z}\bm{e}_{1}+a_{2}k_{z}\bm{e}_{2}+ia_{3}(k_{+}^{3}-k_{-}^{3})\bm{e}_{3}$
$E_{1u}$ |  $a_{1}k_{z}\bm{e}_{+}+a_{3}k_{+}\bm{e}_{3}$, |  $a_{1}k_{-}^{2}k_{z}\bm{e}_{1}+a_{2}k_{-}^{2}k_{z}\bm{e}_{2}+a_{3}k_{+}\bm{e}_{3}$,
| $a_{1}k_{z}\bm{e}_{-}+a_{3}k_{-}\bm{e}_{3}$ | $a_{1}k_{+}^{2}k_{z}\bm{e}_{1}-a_{2}k_{+}^{2}k_{z}\bm{e}_{2}+a_{3}k_{-}\bm{e}_{3}$
$E_{2u}$ |  $a_{1}k_{+}\bm{e}_{+}+a_{3}k_{+}^{2}k_{z}\bm{e}_{3}$, |  $a_{1}k_{+}\bm{e}_{1}+a_{2}k_{+}\bm{e}_{2}+a_{3}k_{-}^{2}k_{z}\bm{e}_{3}$,
| $a_{1}k_{-}\bm{e}_{-}+a_{3}k_{-}^{2}k_{z}\bm{e}_{3}$ | $a_{1}k_{-}\bm{e}_{1}-a_{2}k_{-}\bm{e}_{2}+a_{3}k_{+}^{2}k_{z}\bm{e}_{3}$

####  V.2.1 $\Gamma_{7}$ and $\Gamma_{8}$ bands

In the non-exceptional $\Gamma_{7}$ and $\Gamma_{8}$ bands, one can put
$\mu_{i\nu}(\bm{k})=\mu_{B}\delta_{i\nu}$. Therefore, the susceptibility is
given by Eq. (IV.2) and has the usual temperature dependence. For instance, if
$\bm{d}\parallel\hat{\bm{z}}$, then $\chi_{xx}=\chi_{yy}=\chi_{P}$ and
$\chi_{zz}=\chi_{P}\left\langle{\cal Y}(\bm{k},T)\right\rangle_{FS}$, i.e.,
the basal-plane susceptibility is the same as in the normal state, whereas the
$z$-axis susceptibility is determined by the nodal quasiparticles and has the
asymptotics (60) at $T\to 0$.

####  V.2.2 $\Gamma_{9}$ bands

The Zeeman coupling in the $\Gamma_{9}$ bands has the form (29). If
$\bm{d}\parallel\hat{\bm{z}}$, then we obtain from Eq. (IV.2) that the off-
diagonal components vanish after the Fermi-surface averaging, the basal-plane
susceptibility is not affected by the superconductivity, and the $z$-axis
susceptibility is determined by the nodal quasiparticles:

$\begin{array}[]{c}\chi_{xx}(T)=\chi_{yy}(T)=\chi_{\perp}^{N},\vskip 6.0pt
plus 2.0pt minus 2.0pt\\\ \chi_{zz}(T)=\chi_{zz}^{N}\left\langle{\cal
Y}(\bm{k},T)\right\rangle_{FS},\end{array}$ (65)

where $\chi_{\perp}^{N}$ and $\chi_{zz}^{N}$ are given by Eq. (63). The low-
temperature asymptotics of $\chi_{zz}$ are the same as in Eq. (60). Therefore,
the non-pseudospin character of the $\Gamma_{9}$ bands does not qualitatively
affect the spin response for $\bm{d}\parallel\hat{\bm{z}}$.

In the case $\bm{d}\perp\hat{\bm{z}}$, let us first put $a_{2}=a_{3}=0$ in all
basis functions in the third column of Table 3, which corresponds to
$\bm{d}\parallel\hat{\bm{x}}$. In all stable states, the off-diagonal elements
of the susceptibility tensor vanish and we obtain:

$\displaystyle\chi_{xx}=2N_{F}\bigl{[}\beta_{1}^{2}\left\langle(k_{x}^{2}-k_{y}^{2})^{2}{\cal
Y}(\bm{k},T)\right\rangle_{FS}$ $\displaystyle\hskip
96.73918pt+4\beta_{2}^{2}\left\langle
k_{x}^{2}k_{y}^{2}\right\rangle_{FS}\bigr{]},$
$\displaystyle\chi_{yy}=2N_{F}\bigl{[}4\beta_{1}^{2}\left\langle
k_{x}^{2}k_{y}^{2}{\cal Y}(\bm{k},T)\right\rangle_{FS}$ (66)
$\displaystyle\hskip
76.82234pt+\beta_{2}^{2}\left\langle(k_{x}^{2}-k_{y}^{2})^{2}\right\rangle_{FS}\bigr{]},$
$\displaystyle\chi_{zz}=2N_{F}\mu_{z}^{2}.$

We see that, while $\chi_{zz}$ is not affected by the superconductivity:
$\chi_{zz}(T)=\chi_{zz}^{N}$, the basal-plane components are suppressed in the
superconducting state, but not completely, as both $\chi_{xx}$ and $\chi_{yy}$
attain a nonzero value at $T=0$:

$\chi_{xx}(T=0)=\chi_{yy}(T=0)=\frac{\beta_{2}^{2}}{\beta_{1}^{2}+\beta_{2}^{2}}\chi_{\perp}^{N}.$
(67)

In contrast, in the case of the usual isotropic Zeeman coupling one has
$\chi_{xx}=\chi_{P}\left\langle{\cal Y}(\bm{k},T)\right\rangle_{FS}$ and
$\chi_{yy}=\chi_{zz}=\chi_{P}$.

The temperature dependence of $\chi_{xx}$ and $\chi_{yy}$ is determined by the
Fermi-surface averages in the first terms in Eq. (V.2.2). Since the Yosida
functions are multiplied by strongly anisotropic factors, which vanish in the
vertical planes $|k_{x}|=|k_{y}|$ (in $\chi_{xx}$) or $k_{x}k_{y}=0$ (in
$\chi_{yy}$), the quasiparticle contribution to the susceptibility is
suppressed if the superconducting gap has zeros only in these planes. Using
Table 3, it is easy to see that this happens in the non-chiral “nematic”
$p$-wave state $\bm{d}(\bm{k})\propto k_{x}\bm{e}_{1}$ (or $k_{y}\bm{e}_{1}$),
which corresponds to the $E_{2u}$ irrep and has vertical line nodes. In this
state, we obtain:

$\chi_{xx}(T)-\chi_{xx}(0)\propto T,\quad\chi_{yy}(T)-\chi_{yy}(0)\propto
T^{3}.$ (68)

In the chiral $p$-wave state $\bm{d}(\bm{k})\propto k_{+}\bm{e}_{1}$ (or
$k_{-}\bm{e}_{1}$), the gap has isolated first-order point nodes along the
main axis. The anisotropic Zeeman factors in the first terms in $\chi_{xx}$
and $\chi_{yy}$ vanish along this axis and we obtain:

$\chi_{xx}(T)-\chi_{xx}(0)=\chi_{yy}(T)-\chi_{yy}(0)\propto T^{6}.$ (69)

We see that the spin susceptibility in both chiral and nonchiral $E_{2u}$
states in the $\Gamma_{9}$ bands exhibits the temperature dependence which is
never found in the $\Gamma_{7}$ and $\Gamma_{8}$ bands.

In a more general state with $\bm{d}\perp\hat{\bm{z}}$, the direction of
$\bm{d}$ may change as a function of momentum and the expressions (V.2.2) are
no longer applicable. As an example, let us consider an $f$-wave $A_{2u}$
state with $a_{1}=a_{2}\neq 0$ and $a_{3}=0$, see Table 3. The gap function

$\bm{d}(\bm{k})\propto
i(k_{+}^{3}-k_{-}^{3})\bm{e}_{1}+(k_{+}^{3}+k_{-}^{3})\bm{e}_{2}$ (70)

has two isolated third-order point nodes at the main symmetry axis. From Eq.
(IV.2), we obtain the following nonzero components of the susceptibility
tensor:

$\displaystyle\chi_{xx}=N_{F}\bigl{\langle}\bigl{[}\beta_{1}^{2}(k_{x}^{2}-k_{y}^{2})^{2}+4\beta_{2}^{2}k_{x}^{2}k_{y}^{2}\bigr{]}$
$\displaystyle\hskip 88.2037pt\times[1+{\cal
Y}(\bm{k},T)]\bigr{\rangle}_{FS},$
$\displaystyle\chi_{yy}=N_{F}\bigl{\langle}\bigl{[}4\beta_{1}^{2}k_{x}^{2}k_{y}^{2}+\beta_{2}^{2}(k_{x}^{2}-k_{y}^{2})^{2}\bigr{]}$
$\displaystyle\hskip 88.2037pt\times[1+{\cal
Y}(\bm{k},T)]\bigr{\rangle}_{FS},$ $\displaystyle\chi_{zz}=2N_{F}\mu_{z}^{2}.$

Therefore, $\chi_{zz}$ is the same as in the normal state, whereas the basal-
plane components are suppressed in the superconducting state, but have a
residual value at $T=0$:

$\chi_{xx}(T=0)=\chi_{yy}(T=0)=\frac{1}{2}\chi_{\perp}^{N}.$ (72)

At nonzero temperatures, we calculate the Fermi-surface averages in Eq.
(V.2.2) containing the Yosida function and obtain:

$\chi_{xx}(T)-\chi_{xx}(0)=\chi_{yy}(T)-\chi_{yy}(0)\propto T^{2},$ (73)

which is very different from the $T^{2/3}$ behaviour expected for isolated
third-order nodes in the isotropic Zeeman coupling case, see Eq. (60).

###  V.3 Discussion

The point group $\mathbf{D}_{6h}$ describes, for instance, the symmetry of the
heavy-fermion material UPt3 (Ref. UPt3-review, ), in which a variety of
thermodynamic and transport measurements have revealed an unconventional
superconducting state. Although there is still no consensus on the pairing
symmetry in UPt3, the most promising candidate model is based on the 2D irrep
$E_{2u}$ of $\mathbf{D}_{6h}$. The corresponding order parameter is real in
the high-temperature $A$ phase and complex (TR symmetry breaking) in the low-
temperature $B$ phase.Strand09 ; Schemm14 ; Avers20  The electron bands in
UPt3 originate from the $j_{z}=\pm 1/2$, $\pm 3/2$, and $\pm 5/2$ doublets at
the $\Gamma$ point, which correspond, respectively, to the $\Gamma_{7}$,
$\Gamma_{9}$, and $\Gamma_{8}$ coreps of $\mathbf{D}_{6h}$.

Unfortunately, the NMR Knight shift in UPt3 changes very little in the
superconducting state for all field orientations and the experimental results
“are not conclusive in favor of any theory” (Ref. UPt3-review, ), see also
Ref. Gannon17, . One possible expalanation is that in this material there are
five bands crossing the Fermi level, which alone would lead to a significant
temperature-independent interband contribution, see Eq. (54), regardless of
the pairing symmetry. Additional contribution to the residual susceptibility
can come, e.g., from impurity scattering.AG62

##  VI Conclusions

We have derived the general symmetry-constrained expressions for the Zeeman
coupling in crystals with the SO coupling. We showed that in some Bloch bands
in trigonal and hexagonal crystals the Zeeman interaction of the band
electrons with a magnetic field directed in the basal plane necessarily
vanishes along the main symmetry axis, either linearly or quadratically in
$k_{x}$ and $k_{y}$. This significantly changes the low-temperature behaviour
of the spin susceptibility in the superconducting states with the gap nodes.
For example, in a hexagonal crystal, the basal-plane intraband susceptibility
in the $\Gamma_{9}$ bands has the form $\chi_{xx}(T)=\chi_{yy}(T)\propto
T^{3}$ in the singlet chiral $d$-wave state (instead of the “usual” linear in
$T$ dependence, for the isotropic Zeeman coupling). In the triplet states, one
has $\chi_{xx}(T)-\chi_{xx}(0)\propto T$, $\chi_{yy}(T)-\chi_{yy}(0)\propto
T^{3}$ in the nematic $p$-wave state (instead of the linear in $T$ dependence
for $\chi_{xx}$ and a temperature-independent $\chi_{yy}$), and
$\chi_{xx}(T)-\chi_{xx}(0)=\chi_{yy}(T)-\chi_{yy}(0)\propto T^{6}$ in the
chiral $p$-wave state (instead of the $T^{2}$ dependence for $\chi_{xx}$ and a
temperature-independent $\chi_{yy}$).

###### Acknowledgements.

The author is grateful to V. P. Mineev for a stimulating correspondence. This
work was supported by a Discovery Grant 2015-06656 from the Natural Sciences
and Engineering Research Council of Canada.

##  Appendix A Symmetry of the Bloch states

The Bloch states $|\bm{k},n,1\rangle$ and $|\bm{k},n,2\rangle$ in a crystal
with the SO coupling form the basis of an irreducible double-valued corep of
the magnetic point group of the wave vector $\bm{k}$. The full symmetry group
of $\bm{k}$ is “magnetic”, because it contains the antiunitary conjugation
operation ${\cal C}$. A detailed review of magnetic groups and their coreps
can be found, for instance, in Refs. BD68,  and BC-book, . If the crystal
point group is $\mathbb{G}$, then the magnetic group at the $\Gamma$ point is
${\cal G}=\mathbb{G}+{\cal C}\mathbb{G}$. Since we consider only crystals with
a center of inversion, the coreps of ${\cal G}$ have a definite parity, i.e.,
are either inversion-even ($\Gamma^{+}$) or inversion-odd ($\Gamma^{-}$).

The double-valued coreps of ${\cal G}$ can be obtained from the double-valued
irreducible representations (irreps) of $\mathbb{G}$, with the results listed
in Table 4. All the coreps are 2D, except
$(\Gamma_{6}^{\pm},\Gamma_{7}^{\pm})$ for $\mathbb{G}=\mathbf{T}_{h}$ and
$\Gamma_{8}^{\pm}$ for $\mathbb{G}=\mathbf{O}_{h}$, which are four-
dimensional. These latter coreps correspond to the bands which are fourfold
degenerate at the $\Gamma$ point, e.g., the $\Gamma_{8}^{\pm}$ (“$j=3/2$”)
bands in cubic crystals (Ref. Lutt56, ), and are not considered here. Pairs of
complex conjugate irreps $(\Gamma,\Gamma^{*})$ produce coreps of the “pairing”
type of twice the dimension, while the 1D irreps $\Gamma_{2}^{\pm}$ of
$\mathbf{C}_{i}$ and $\Gamma_{6}^{\pm}$ of $\mathbf{C}_{3i}$ produce 2D coreps
of the “doubling” type.Lax-book

The fourth and fifth columns in Table 4 show the spinor basis functions
$(\phi_{1},\phi_{2})$ which reproduce the corep matrices in Table 1. For the
pseudospin coreps, the basis functions can be chosen to be the pure spin
states $\lvert\uparrow\rangle$ and $\lvert\downarrow\rangle$ (the eigenstates
of the spin operator $\hat{s}_{z}$), whereas the basis functions for the non-
pseudospin coreps contain the pure spin states multiplied by additional non-
removable coordinate-dependent factors. The phases are chosen to satisfy the
conjugation condition $\phi_{2}={\cal C}\phi_{1}$.

Transformation of the Bloch states at the $\Gamma$ point in a twofold
degenerate band corresponding to a 2D corep ${\cal D}_{n}$ is given by

$g|\bm{0},n,s\rangle=\sum_{s^{\prime}}|\bm{0},n,s^{\prime}\rangle{\cal
D}_{n,s^{\prime}s}(g).$ (74)

The explicit form of the corep matrices depends on the choice of the basis: an
arbitrary unitary rotation of the basis,
$|\bm{0},n,s\rangle\to|\bm{0},n,s\rangle^{\prime}=\sum_{s_{1}}|\bm{0},n,s_{1}\rangle
U_{n,s_{1}s}$, produces an equivalent corep with $\hat{\cal
D}^{\prime}_{n}(g)=\hat{U}_{n}^{-1}\hat{\cal D}_{n}(g)\hat{U}_{n}$ and
$\hat{\cal D}^{\prime}_{n}({\cal C})=\hat{U}_{n}^{-1}\hat{\cal D}_{n}({\cal
C})\hat{U}_{n}^{*}$. The “orientation” and the phases of the Bloch basis at
the $\Gamma$-point can always be chosen to reproduce both the matrices in
Table 1 and the matrix representation of the conjugation, $\hat{\cal
D}_{n}({\cal C})=-i\hat{\sigma}_{2}$. The prescription (4) is obtained from
Eq. (74) by using the fact that the Bloch states in the bands that are only
twofold degenerate due to the conjugation symmetry are analytic functions of
$\bm{k}$, at least in the vicinity of the $\Gamma$ point, which can be shown,
e.g., using the $\bm{k}\cdot\bm{p}$ perturbation theory.

Table 4: The double-valued coreps of the centrosymmetric magnetic point groups at the $\Gamma$ point, with examples of even and odd spinor basis functions for the 2D coreps; $f_{\Gamma}(\bm{r})$ is a real basis function of a single-valued irrep $\Gamma$ of $\mathbb{G}$, and $\rho_{\pm}=x\pm iy$. $\mathbb{G}$ | corep | $\dim$ | even basis | odd basis
---|---|---|---|---
$\mathbf{C}_{i}$ | $\Gamma_{2}$ | 2 | $(\lvert\uparrow\rangle,\lvert\downarrow\rangle)$ | $if_{\Gamma_{1}^{-}}(\bm{r})(\lvert\uparrow\rangle,\lvert\downarrow\rangle)$
$\mathbf{C}_{2h}$ | $(\Gamma_{3},\Gamma_{4})$ | 2 | $(\lvert\uparrow\rangle,\lvert\downarrow\rangle)$ | $if_{\Gamma_{1}^{-}}(\bm{r})(\lvert\uparrow\rangle,\lvert\downarrow\rangle)$
$\mathbf{D}_{2h}$ | $\Gamma_{5}$ | 2 | $(\lvert\uparrow\rangle,\lvert\downarrow\rangle)$ | $if_{\Gamma_{1}^{-}}(\bm{r})(\lvert\uparrow\rangle,\lvert\downarrow\rangle)$
$\mathbf{C}_{4h}$ | $(\Gamma_{5},\Gamma_{6})$ | 2 | $(\lvert\uparrow\rangle,\lvert\downarrow\rangle)$ | $if_{\Gamma_{1}^{-}}(\bm{r})(\lvert\uparrow\rangle,\lvert\downarrow\rangle)$
| $(\Gamma_{7},\Gamma_{8})$ | 2 | $(\rho_{-}^{2}\lvert\uparrow\rangle,\rho_{+}^{2}\lvert\downarrow\rangle)$ | $(\rho_{+}\lvert\uparrow\rangle,-\rho_{-}\lvert\downarrow\rangle)$
$\mathbf{D}_{4h}$ | $\Gamma_{6}$ | 2 | $(\lvert\uparrow\rangle,\lvert\downarrow\rangle)$ | $if_{\Gamma_{1}^{-}}(\bm{r})(\lvert\uparrow\rangle,\lvert\downarrow\rangle)$
| $\Gamma_{7}$ | 2 | $(\rho_{-}^{2}\lvert\uparrow\rangle,\rho_{+}^{2}\lvert\downarrow\rangle)$ | $(\rho_{+}\lvert\uparrow\rangle,-\rho_{-}\lvert\downarrow\rangle)$
$\mathbf{C}_{3i}$ | $(\Gamma_{4},\Gamma_{5})$ | 2 | $(\lvert\uparrow\rangle,\lvert\downarrow\rangle)$ | $if_{\Gamma_{1}^{-}}(\bm{r})(\lvert\uparrow\rangle,\lvert\downarrow\rangle)$
| $\Gamma_{6}$ | 2 | $(\rho_{-}^{2}\lvert\uparrow\rangle,\rho_{+}^{2}\lvert\downarrow\rangle)$ | $(\rho_{+}\lvert\uparrow\rangle,-\rho_{-}\lvert\downarrow\rangle)$
$\mathbf{D}_{3d}$ | $\Gamma_{4}$ | 2 | $(\lvert\uparrow\rangle,\lvert\downarrow\rangle)$ | $if_{\Gamma_{1}^{-}}(\bm{r})(\lvert\uparrow\rangle,\lvert\downarrow\rangle)$
| $(\Gamma_{5},\Gamma_{6})$ | 2 | $(\rho_{-}^{2}\lvert\uparrow\rangle,\rho_{+}^{2}\lvert\downarrow\rangle)$ | $(\rho_{+}\lvert\uparrow\rangle,-\rho_{-}\lvert\downarrow\rangle)$
$\mathbf{C}_{6h}$ | $(\Gamma_{7},\Gamma_{8})$ | 2 | $(\lvert\uparrow\rangle,\lvert\downarrow\rangle)$ | $if_{\Gamma_{1}^{-}}(\bm{r})(\lvert\uparrow\rangle,\lvert\downarrow\rangle)$
| $(\Gamma_{9},\Gamma_{10})$ | 2 | $(\rho_{+}^{2}\lvert\uparrow\rangle,\rho_{-}^{2}\lvert\downarrow\rangle)$ | $(\rho_{-}^{3}\lvert\uparrow\rangle,-\rho_{+}^{3}\lvert\downarrow\rangle)$
| $(\Gamma_{11},\Gamma_{12})$ | 2 | $(\rho_{-}^{2}\lvert\uparrow\rangle,\rho_{+}^{2}\lvert\downarrow\rangle)$ | $(\rho_{+}\lvert\uparrow\rangle,-\rho_{-}\lvert\downarrow\rangle)$
$\mathbf{D}_{6h}$ | $\Gamma_{7}$ | 2 | $(\lvert\uparrow\rangle,\lvert\downarrow\rangle)$ | $if_{\Gamma_{1}^{-}}(\bm{r})(\lvert\uparrow\rangle,\lvert\downarrow\rangle)$
| $\Gamma_{8}$ | 2 | $(\rho_{+}^{2}\lvert\uparrow\rangle,\rho_{-}^{2}\lvert\downarrow\rangle)$ | $(\rho_{-}^{3}\lvert\uparrow\rangle,-\rho_{+}^{3}\lvert\downarrow\rangle)$
| $\Gamma_{9}$ | 2 | $(\rho_{-}^{2}\lvert\uparrow\rangle,\rho_{+}^{2}\lvert\downarrow\rangle)$ | $(\rho_{+}\lvert\uparrow\rangle,-\rho_{-}\lvert\downarrow\rangle)$
$\mathbf{T}_{h}$ | $\Gamma_{5}$ | 2 | $(\lvert\uparrow\rangle,\lvert\downarrow\rangle)$ | $if_{\Gamma_{1}^{-}}(\bm{r})(\lvert\uparrow\rangle,\lvert\downarrow\rangle)$
| $(\Gamma_{6},\Gamma_{7})$ | 4 | - | -
$\mathbf{O}_{h}$ | $\Gamma_{6}$ | 2 | $(\lvert\uparrow\rangle,\lvert\downarrow\rangle)$ | $if_{\Gamma_{1}^{-}}(\bm{r})(\lvert\uparrow\rangle,\lvert\downarrow\rangle)$
| $\Gamma_{7}$ | 2 | $f_{\Gamma_{2}^{+}}(\bm{r})(\lvert\uparrow\rangle,\lvert\downarrow\rangle)$ | $if_{\Gamma_{2}^{-}}(\bm{r})(\lvert\uparrow\rangle,\lvert\downarrow\rangle)$
| $\Gamma_{8}$ | 4 | - | -

##  Appendix B Green’s functions

We introduce the normal and anomalous Green’s functions in the Matsubara
representation:AGD

$\displaystyle G_{nn^{\prime},ss^{\prime}}(\bm{k},\tau)=-\langle
T_{\tau}c_{\bm{k},n,s}(\tau)c^{\dagger}_{\bm{k},n^{\prime},s^{\prime}}(0)\rangle,$
$\displaystyle F_{nn^{\prime},ss^{\prime}}(\bm{k},\tau)=\langle
T_{\tau}c_{\bm{k},n,s}(\tau)\tilde{c}_{\bm{k},n^{\prime},s^{\prime}}(0)\rangle,$
(75) $\displaystyle\bar{F}_{nn^{\prime},ss^{\prime}}(\bm{k},\tau)=\langle
T_{\tau}\tilde{c}^{\dagger}_{\bm{k},n,s}(\tau)c^{\dagger}_{\bm{k},n^{\prime},s^{\prime}}(0)\rangle.$

These Green’s functions form $2N\times 2N$ matrices, denoted below by
$\check{G}$, etc, in the direct product of the band and Kramers spaces. We
reserve the notation $\hat{G}$ for $2\times 2$ matrices in the Kramers space.
Introducing the four-component Nambu operators in the $n$th band,

$C_{\bm{k},n}=\left(\begin{array}[]{c}c_{\bm{k},n,1}\\\ c_{\bm{k},n,2}\\\
\tilde{c}^{\dagger}_{\bm{k},n,1}\\\ \tilde{c}^{\dagger}_{\bm{k},n,2}\\\
\end{array}\right)=\left(\begin{array}[]{c}c_{\bm{k},n,1}\\\ c_{\bm{k},n,2}\\\
p_{n}c^{\dagger}_{-\bm{k},n,2}\\\ -p_{n}c^{\dagger}_{-\bm{k},n,1}\\\
\end{array}\right),$

one can combine the expressions (B) into

$\displaystyle{\cal G}_{nn^{\prime}}(\bm{k},\tau)$ $\displaystyle=$
$\displaystyle-\langle
T_{\tau}C_{\bm{k},n}(\tau)C^{\dagger}_{\bm{k},n^{\prime}}(0)\rangle$ (78)
$\displaystyle=$
$\displaystyle\left(\begin{array}[]{cc}\hat{G}_{nn^{\prime}}(\bm{k},\tau)&-\hat{F}_{nn^{\prime}}(\bm{k},\tau)\\\
-\hat{\bar{F}}_{nn^{\prime}}(\bm{k},\tau)&\hat{\bar{G}}_{nn^{\prime}}(\bm{k},\tau)\\\
\end{array}\right),$

where the auxiliary normal Green’s function is given by

$\displaystyle\bar{G}_{nn^{\prime},ss^{\prime}}(\bm{k},\tau)$ $\displaystyle=$
$\displaystyle-\langle
T_{\tau}\tilde{c}^{\dagger}_{\bm{k},n,s}(\tau)\tilde{c}_{\bm{k},n^{\prime},s^{\prime}}(0)\rangle$
$\displaystyle=$ $\displaystyle-
p_{n}p_{n^{\prime}}[\hat{\sigma}_{2}\hat{G}_{n^{\prime}n}(-\bm{k},-\tau)\hat{\sigma}_{2}]_{s^{\prime}s},$

using Eq. (40).

The Green’s functions satisfy the Gor’kov equations obtained from the
Hamiltonian (36) using the standard procedure:AGD

$\displaystyle-\frac{\partial{\cal G}(\bm{k},\tau)}{\partial\tau}$
$\displaystyle=$ $\displaystyle\delta(\tau)\mathbb{1}_{4N}$ (81)
$\displaystyle+\left(\begin{array}[]{cc}\check{\epsilon}(\bm{k})&\check{\Delta}(\bm{k})\\\
\check{\Delta}^{\dagger}(\bm{k})&-\check{\tilde{\epsilon}}(\bm{k})\end{array}\right){\cal
G}(\bm{k},\tau),$

where $\mathbb{1}_{4N}$ is the $4N\times 4N$ unit matrix,

$\displaystyle\epsilon_{nn^{\prime},ss^{\prime}}(\bm{k})=\xi_{n}(\bm{k})\delta_{nn^{\prime}}\delta_{ss^{\prime}}-\bm{H}\bm{m}_{nn^{\prime},ss^{\prime}}(\bm{k}),$
$\displaystyle\tilde{\epsilon}_{nn^{\prime},ss^{\prime}}(\bm{k})=\xi_{n}(\bm{k})\delta_{nn^{\prime}}\delta_{ss^{\prime}}+\bm{H}\bm{m}_{nn^{\prime},ss^{\prime}}(\bm{k}),$

and
$\Delta_{nn^{\prime},ss^{\prime}}(\bm{k})=\Delta_{n,ss^{\prime}}(\bm{k})\delta_{nn^{\prime}}$.
Note that the gap functions are band-diagonal (since we neglected the
interband pairing), whereas the single-particle energy $\check{\epsilon}$ and
its time-reversed counterpart $\check{\tilde{\epsilon}}$ are not. Using the
Matsubara frequency representation, ${\cal G}(\bm{k},\tau)=T\sum_{m}{\cal
G}(\bm{k},\omega_{m})e^{-i\omega_{m}\tau}$, where $\omega_{m}=(2m+1)\pi T$,
one can solve Eq. (81), with the following result:

$\displaystyle\check{G}=[i\omega_{m}-\check{\epsilon}-\check{\Delta}(i\omega_{m}+\check{\tilde{\epsilon}})^{-1}\check{\Delta}^{\dagger}]^{-1},$
(82)
$\displaystyle\check{F}=-(i\omega_{m}-\check{\epsilon})^{-1}\check{\Delta}\check{\bar{G}},$
$\displaystyle\check{\bar{F}}=-(i\omega_{m}+\check{\tilde{\epsilon}})^{-1}\check{\Delta}^{\dagger}\check{G},$
$\displaystyle\check{\bar{G}}=[i\omega_{m}+\check{\tilde{\epsilon}}-\check{\Delta}^{\dagger}(i\omega_{m}-\check{\epsilon})^{-1}\check{\Delta}]^{-1}.$

These expressions are valid at arbitrary magnetic field.

At zero field, $\check{\epsilon}$ and $\check{\tilde{\epsilon}}$ become band-
diagonal. Therefore, the Green’s functions are also band-diagonal:
$G^{(0)}_{nn^{\prime},ss^{\prime}}=\delta_{nn^{\prime}}G^{(0)}_{n,ss^{\prime}}$,
etc, where

$\displaystyle\hat{G}^{(0)}_{n}$ $\displaystyle=$
$\displaystyle-(i\omega_{m}+\xi_{n})(\omega_{m}^{2}+\xi_{n}^{2}+\hat{\Delta}_{n}\hat{\Delta}_{n}^{\dagger})^{-1},$
$\displaystyle\hat{F}^{(0)}_{n}$ $\displaystyle=$
$\displaystyle\hat{\Delta}_{n}(\omega_{m}^{2}+\xi_{n}^{2}+\hat{\Delta}_{n}^{\dagger}\hat{\Delta}_{n})^{-1}$
$\displaystyle=$
$\displaystyle(\omega_{m}^{2}+\xi_{n}^{2}+\hat{\Delta}_{n}\hat{\Delta}_{n}^{\dagger})^{-1}\hat{\Delta}_{n},$
$\displaystyle\hat{\bar{F}}^{(0)}_{n}$ $\displaystyle=$
$\displaystyle\hat{\Delta}_{n}^{\dagger}(\omega_{m}^{2}+\xi_{n}^{2}+\hat{\Delta}_{n}\hat{\Delta}_{n}^{\dagger})^{-1}=[\hat{F}^{(0)}_{n}]^{\dagger}.$

In a weak field, one can expand Eq. (82), with the following result:

$\check{G}=\check{G}^{(0)}-\check{G}^{(0)}(\bm{H}\check{\bm{m}})\check{G}^{(0)}-\check{F}^{(0)}(\bm{H}\check{\bm{m}})\check{\bar{F}}^{(0)}+{\cal
O}(H^{2}).$

Substituting this into Eq. (44) and dropping the superscript “$(0)$”, we
arrive at the expressions (45) and (IV.1) for the spontaneous magnetization
and the spin susceptibility, respectively.

##  Appendix C Interband susceptibility

We consider just one pair of bands ($n,n^{\prime}=1,2$) with
$\xi_{2}(\bm{k})>\xi_{1}(\bm{k})$ and focus on the terms with $n\neq
n^{\prime}$ in Eq. (IV.1). Using Eq. (52), we obtain the following expression
for the interband contribution to the spin susceptibility:

$\displaystyle\tilde{\chi}_{ij}=4T\sum_{m}\int\frac{d^{3}\bm{k}}{(2\pi)^{3}}\bigl{[}P_{ij}(\omega_{m}^{2}-\xi_{1}\xi_{2})-Q_{ij}\bigr{]}$
$\displaystyle\times\frac{1}{(\omega_{m}^{2}+\xi_{1}^{2}+|\Delta_{1}|^{2})(\omega_{m}^{2}+\xi_{2}^{2}+|\Delta_{2}|^{2})},$
(83)

where $P_{ij}=a_{i}a_{j}+\bm{b}_{i}\bm{b}_{j}$,

$Q_{ij}=(a_{i}a_{j}+\bm{b}_{i}\bm{b}_{j})\,\mathrm{Re}\,(\psi_{1}^{*}\psi_{2})$

in the singlet case, and

$\displaystyle Q_{ij}$ $\displaystyle=$
$\displaystyle(a_{i}a_{j}-\bm{b}_{i}\bm{b}_{j})\,\mathrm{Re}\,(\bm{d}_{1}^{*}\bm{d}_{2})$
$\displaystyle-(a_{i}\bm{b}_{j}+a_{j}\bm{b}_{i})\,\mathrm{Re}\,(\bm{d}_{1}^{*}\times\bm{d}_{2})$
$\displaystyle+\,\mathrm{Re}\,(\bm{b}_{i}\bm{d}_{1}^{*})(\bm{b}_{j}\bm{d}_{2})+\,\mathrm{Re}\,(\bm{b}_{i}\bm{d}_{2}^{*})(\bm{b}_{j}\bm{d}_{1})$

in the unitary triplet case. In the normal state, the interband susceptibility
(C) takes the form (53). To calculate the interband susceptibility in the
superconducting state, we assume the weak-coupling picture, so that the gap
functions are nonzero only inside the BCS shells near the Fermi surfaces:
$\psi_{n}(\bm{k}),\bm{d}_{n}(\bm{k})\propto\theta(\epsilon_{c}-|\xi_{n}(\bm{k})|)$.
For a large band splitting, ${\cal E}(\bm{k})\gg\epsilon_{c}$, the BCS shells
in the bands 1 and 2 do not overlap and therefore $Q_{ij}=0$ in both singlet
and triplet cases.

One can expect that the interband susceptibility is almost unchanged when the
system undergoes a superconducting phase transition, in which only the
electrons near the Fermi surfaces are affected. To find the upper bound on the
effect of superconductivity on $\tilde{\chi}_{ij}$, it is sufficient to
calculate the latter at $T=0$, where the Matsubara sum in Eq. (C) can be
replaced by a frequency integral. In this way, we obtain:

$\displaystyle\tilde{\chi}_{ij}^{S}(T=0)-\tilde{\chi}_{ij}^{N}=2\int\frac{d^{3}\bm{k}}{(2\pi)^{3}}P_{ij}$
$\displaystyle\quad\times\biggl{[}\left(1-\frac{\xi_{1}\xi_{2}}{E_{1}E_{2}}\right)\frac{1}{E_{1}+E_{2}}-\frac{1-\mathrm{sign}(\xi_{1}\xi_{2})}{|\xi_{1}|+|\xi_{2}|}\biggr{]},$

where the excitation energies $E_{1,2}$ are given by Eq. (48). The expression
in the square brackets is nonzero only near the two Fermi surfaces, therefore
the susceptibility deviation can be written as a sum of the independent
contributions from the two BCS shells:
$\tilde{\chi}_{ij}^{S}-\tilde{\chi}_{ij}^{N}={\cal I}_{ij}^{(1)}+{\cal
I}_{ij}^{(2)}$, where

$\displaystyle{\cal
I}_{ij}^{(1)}=2\int\frac{d^{3}\bm{k}}{(2\pi)^{3}}\theta(\epsilon_{c}-|\xi_{1}|)P_{ij}$
$\displaystyle\quad\times\left[\left(1-\frac{\xi_{1}}{E_{1}}\mathrm{sign}\,\xi_{2}\right)\frac{1}{E_{1}+|\xi_{2}|}-\frac{1-\mathrm{sign}(\xi_{1}\xi_{2})}{|\xi_{1}|+|\xi_{2}|}\right]$

and ${\cal I}_{ij}^{(2)}$ is obtained from ${\cal I}_{ij}^{(1)}$ by replacing
$\xi_{1}\leftrightarrow\xi_{2}$, $E_{1}\leftrightarrow E_{2}$. Note that
$\xi_{2}>0$ in the BCS shell in band $1$, while $\xi_{1}<0$ in the BCS shell
in band $2$. Next, we neglect the energy dependence of the single-particle
DoS, of $P_{ij}$, and of ${\cal E}$ near the Fermi surface and obtain:

${\cal I}_{ij}^{(1)}=2N_{F,1}\langle P_{ij}I_{1}\rangle_{FS,1},$ (84)

where

$I_{1}=\int_{-x_{m}}^{x_{m}}dx\left(\frac{1-x/\sqrt{x^{2}+\delta^{2}}}{1+x+\sqrt{x^{2}+\delta^{2}}}-\frac{1-\mathrm{sign}\,x}{1+x+|x|}\right),$

$x_{m}=\epsilon_{c}/{\cal E}$, and $\delta=|\Delta_{1}|/{\cal E}$. Since
$\delta\ll x_{m}\ll 1$, the last integral can be easily calculated with the
logarithmic accuracy: $I_{1}\simeq-2\delta^{2}\ln(x_{m}/\delta)$. Substituting
this into Eq. (84) and repeating the calculation for the second band, we
finally obtain:
$\tilde{\chi}_{ij}^{S}(T=0)-\tilde{\chi}_{ij}^{N}\propto(|\Delta|^{2}/{\cal
E}^{2})\ln(\epsilon_{c}/|\Delta|)$.

##  Appendix D Derivation of Eq. (IV.2)

Substituting the Green’s functions (47) in Eq. (IV.1), the intraband
susceptibility in the triplet state can be represented in the form
$\chi_{ij}=\chi^{(1)}_{ij}+\chi^{(2)}_{ij}$, where

$\chi^{(1)}_{ij}=2T\sum_{m}\int\frac{d^{3}\bm{k}}{(2\pi)^{3}}\frac{\omega_{m}^{2}-\xi^{2}-|\bm{d}|^{2}}{(\omega_{m}^{2}+\xi^{2}+|\bm{d}|^{2})^{2}}\bigl{(}\hat{\mu}\hat{\mu}^{\dagger}\bigr{)}_{ij}$

and

$\chi^{(2)}_{ij}=4T\sum_{m}\int\frac{d^{3}\bm{k}}{(2\pi)^{3}}\frac{|\bm{d}|^{2}}{(\omega_{m}^{2}+\xi^{2}+|\bm{d}|^{2})^{2}}\bigl{(}\hat{\mu}\hat{\Pi}_{\perp}\hat{\mu}^{\dagger}\bigr{)}_{ij},$

where
$\Pi_{\perp,\mu\nu}=\delta_{\mu\nu}-\,\mathrm{Re}\,(d_{\mu}d_{\nu}^{*})/|\bm{d}|^{2}$.
Summing over the Matsubara frequencies in $\chi^{(1)}$, we obtain:

$\chi^{(1)}_{ij}=2N_{F}\bigl{\langle}\bigl{(}\hat{\mu}\hat{\mu}^{\dagger}\bigr{)}_{ij}{\cal
Y}\bigr{\rangle}_{FS},$ (85)

where ${\cal Y}$ is the Yosida function (56) with
$|\Delta(\bm{k})|=|\bm{d}(\bm{k})|$.

To express $\chi^{(2)}$ in terms of the Yosida function, we use the fact that
$\bm{d}(\bm{k})$ is nonzero only inside the BCS shell near the Fermi surface,
which allows one to write

$\chi^{(2)}_{ij}=2N_{F}\langle(\hat{\mu}\hat{\Pi}_{\perp}\hat{\mu}^{\dagger})_{ij}(J_{+}-J_{-})\rangle_{FS},$

where

$J_{\pm}=T\sum_{m}\int_{-\infty}^{\infty}d\xi\frac{\omega_{m}^{2}-\xi^{2}\pm|\bm{d}|^{2}}{(\omega_{m}^{2}+\xi^{2}+|\bm{d}|^{2})^{2}}.$

It is easy to see that $J_{+}=1$ (this can be shown by adding and subtracting
the value of $J_{+}$ in the normal state and performing the energy integration
before the Matsubara summation)) and $J_{-}={\cal Y}$, therefore

$\chi^{(2)}_{ij}=2N_{F}\bigl{\langle}\bigl{(}\hat{\mu}\hat{\Pi}_{\perp}\hat{\mu}^{\dagger}\bigr{)}_{ij}(1-{\cal
Y})\bigr{\rangle}_{FS}.$ (86)

Combining this with Eq. (85), we obtain:

$\displaystyle\chi_{ij}(T)=2N_{F}\left\langle\bigl{[}\hat{\mu}(\bm{k})\hat{\mu}^{\dagger}(\bm{k})\bigr{]}_{ij}{\cal
Y}(\bm{k},T)\right\rangle_{FS}$
$\displaystyle+2N_{F}\left\langle\bigl{[}\hat{\mu}(\bm{k})\hat{\Pi}_{\perp}(\bm{k})\hat{\mu}^{\dagger}(\bm{k})\bigr{]}_{ij}\bigl{[}1-{\cal
Y}(\bm{k},T)\bigr{]}\right\rangle_{FS},$ (87)

which becomes Eq. (IV.2) after rearranging the terms.

## References

  * (1) M. Sigrist and K. Ueda, Rev. Mod. Phys. 63, 239 (1991).
  * (2) V. P. Mineev and K. V. Samokhin, Introduction to Unconventional Superconductivity (Gordon and Breach, London, 1999).
  * (3) C. Kittel, Quantum Theory of Solids (Wiley, 1987).
  * (4) M. H. Cohen and E. I. Blount, Phil. Mag. 5, 115 (1960).
  * (5) L. M. Roth, Phys. Rev. 145, 434 (1966).
  * (6) E. M. Lifshitz and L. P. Pitaevskii, Statistical Physics, Part 2, Sec. 56 (Butterworth-Heinemann, Oxford, 2002).
  * (7) R. Winkler, Spin-orbit Coupling Effects in Two-Dimensional Electron and Hole Systems (Springer, Berlin, 2003).
  * (8) Y. Wan and Q.-H. Wang, Europhys. Lett. 85, 57007 (2009).
  * (9) M. H. Fischer, New J. Phys. 15, 073006 (2013).
  * (10) A. Ramires and M. Sigrist, Phys. Rev. B 94, 104501 (2016).
  * (11) T. Nomoto, K. Hattori, and H. Ikeda, Phys. Rev. B 94, 174513 (2016).
  * (12) H. Kim, K. Wang, Y. Nakajima, R. Hu, S. Ziemak, P. Syers, L. Wang, H. Hodovanets, J. D. Denlinger, P. M. R. Brydon, D. F. Agterberg, M. A. Tanatar, R. Prozorov, and J. Paglione, Sci. Adv. 4, eaao4513 (2018).
  * (13) K. V. Samokhin, Ann. Phys. (N. Y.) 407, 179 (2019).
  * (14) K. V. Samokhin, Phys. Rev. B 100, 054501 (2019).
  * (15) P. W. Anderson, Phys. Rev. B 30, 4000 (1984).
  * (16) K. Ueda and T. M. Rice, Phys. Rev. B 31, 7114 (1985).
  * (17) Non-centrosymmetric Superconductors: Introduction and Overview, ed. by E. Bauer and M. Sigrist, Lecture Notes in Physics 847 (Springer, Heidelberg, 2012).
  * (18) L. D. Landau and E. M. Lifshitz, Quantum Mechanics (Butterworth-Heinemann, Oxford, 2002).
  * (19) G. L. Bir and G. E. Pikus, Symmetry and Strain-induced Effects in Semiconductors, Sec. 25 (Wiley, New York, 1974).
  * (20) P. W. Anderson, J. Phys. Chem. Solids 11, 26 (1959).
  * (21) E. I. Blount, Phys. Rev. B 32, 2935 (1985).
  * (22) A. J. Leggett, Rev. Mod. Phys. 47, 331 (1975).
  * (23) M. Lax, Symmetry Principles in Solid State and Molecular Physics (Dover Publications, New York, 2001).
  * (24) R. Joynt and L. Taillefer, Rev. Mod. Phys. 74, 235 (2002).
  * (25) J. D. Strand, D. J. Van Harlingen, J. B. Kycia, and W. P. Halperin, Phys. Rev. Lett. 103, 197002 (2009).
  * (26) E. R. Schemm, W. J. Gannon, C. M. Wishne, W. P. Halperin, and A. Kapitulnik, Science 345, 190 (2014).
  * (27) K. E. Avers, W. J. Gannon, S. J. Kuhn, W. P. Halperin, J. A. Sauls, L. DeBeer-Schmitt, C. D. Dewhurst, J. Gavilano, G. Nagy, U. Gasser, and M. R. Eskildsen, Nat. Phys. 16, 531 (2020).
  * (28) W. J. Gannon, W. P. Halperin, M. R. Eskildsen, P. Dai, U. B. Hansen, K. Lefmann, and A. Stunault, Phys. Rev. B 96, 041111(R) (2017).
  * (29) A. A. Abrikosov and L. P. Gor’kov, Zh. Eksp. Teor. Fiz. 42, 1088 (1962) [Sov. Phys. JETP 15, 752 (1962)].
  * (30) C. J. Bradley and B. L. Davies, Rev. Mod. Phys. 40, 359 (1968).
  * (31) C. J. Bradley and A. P. Cracknell, The Mathematical Theory of Symmetry in Solids (Oxford University Press, Oxford, 2010).
  * (32) J. M. Luttinger, Phys. Rev. 102, 1030 (1956).
  * (33) A. A. Abrikosov, L. P. Gorkov, and I. E. Dzyaloshinski, Methods of Quantum Field Theory in Statistical Physics (Dover Publications, New York, 1975).