# Pair-breaking in superconductors with strong spin-orbit coupling

D. C. Cavanagh  Department of Physics, University of Otago, P.O. Box 56,
Dunedin 9054, New Zealand  Daniel F. Agterberg<EMAIL_ADDRESS>Department of
Physics, University of Wisconsin, Milwaukee, Wisconsin 53201, USA  P. M. R.
Brydon<EMAIL_ADDRESS>Department of Physics and MacDiarmid
Institute for Advanced Materials and Nanotechnology, University of Otago, P.O.
Box 56, Dunedin 9054, New Zealand

(July 4, 2022)

###### Abstract

We study the influence of symmetry-breaking perturbations on superconductivity
in multiorbital materials, with a particular focus on an external magnetic
field. We introduce the field-fitness function which characterizes the pair-
breaking effects of the perturbation on a given superconducting state. For
even parity superconductors we find that this field-fitness function for an
external magnetic field is one, implying that the paramagnetic response is
controlled only by a generalized effective $g$-factor. For odd parity
superconductors, the interplay of the effective $g$-factor and the field-
fitness function can lead to counter-intuitive results. We demonstrate this
for $p$-wave pairing in the effective $j=3/2$ electronic states of the
Luttinger-Kohn model.

Introduction. A diverse variety of superconductors have recently been found to
exhibit critical fields far exceeding the Pauli limiting field, e.g. UTe2 [1],
CeRh2As2 [2], UCoGe[3], URhGe[3], and YbRh2Si2 [4]. In some of these materials
superconductivity even appears as a re-entrant phase [3, 1], again far above
the Pauli field. The origin of this anomalous high-field behaviour has been
attributed to spin-triplet superconductivity. However, the fermiology of these
materials is complicated, with multiple bands crossing the Fermi surface.
Moreover, spin-orbit coupling is expected to be large due to the presence of
heavy elements. Due to the interplay of the normal-state band structure, the
structure of the odd-parity pairing potential, and the applied magnetic field,
it is not clear that the established theory for the magnetic response of a
triplet superconductor is applicable [5, 6].

Concurrent to these experimental developments, it has been realized that
internal degrees of freedom of the band electrons, e.g. sublattice or orbital,
can profoundly impact the magnetic response of even-parity superconductors. In
particular, the large critical fields observed in artificial Rashba
heterostructures [7, 8], CeRh2As2 [2, 9], and WTe2 [10] are believed to arise
from a “hidden” antisymmetric spin-orbit coupling (ASOC) [11]. This ASOC is
odd in momentum and has opposite signs for internal degrees of freedom that
are related by inversion symmetry (IS). This preserves IS and the two-fold
degeneracy of the band electron states. Similar to noncentrosymmetric
materials, [12] however, the ASOC reduces the Zeeman splitting of the band
states and so enhances the Pauli limit [13]. For odd-parity superconductivity,
however, the effect of the ASOC on the critical fields is not as well
explored. However, in one remarkable example, CeRh2As2, spin-singlet pairing
interactions give way to extremely high critical fields due to the formation
of an odd-parity superconducting state that is stabilized by the ASOC [14, 2].

There also exist materials where a symmetric spin-orbit coupling (SSOC), i.e.
the spin-orbit coupling is even in momentum, is important. This influence of
this SSOC on superconductivity has been studied in the context of the iron
pnictides [15], Sr2RuO4 [16, 17], half-Heusler materials, [18, 19] and the
pyrochlore lattice [20]. The latter two cases support effective $j=3/2$
electronic states which exhibit properties that are quite different from the
more usual $j=1/2$ states [21]. The magnetic response of even-parity
superconducting states in such materials shows similar features to the case of
the ASOC [22]; the response of odd-parity superconductivity also remains
poorly understood. Indeed, one of the main results of this work is to reveal
the counter-intuitive response of odd-parity states in $j=3/2$ materials to
applied fields.

In this article we examine the influence of the spin-orbit coupling on the
response of a superconducting state to a perturbation which breaks either IS
or time-reversal symmetry (TRS), with a focus on the familiar example of an
applied magnetic field. Within a general minimal model for systems with both
antisymmetric and symmetric spin-orbit coupling, we show that the response of
the superconducting state to the perturbation is fully determined by two
basis-independent quantities: a generalized effective $g$-factor, and a
parameter that quantifies the pair-breaking due to the field, which we term
the field-fitness function $F_{h}$ by analogy with the superconducting fitness
[16, 23]. In the case of an external magnetic field, these quantities also
control the spin susceptibility in the superconducting state. For even-parity
superconductors, $F_{h}=1$, so that the response is given solely by the
effective $g$-factor; in contrast, odd-parity superconducting states display a
complicated interplay of the field-fitness and the effective $g$-factor. We
apply our general theory to odd-parity superconductivity in $j=3/2$ materials,
where we find that the SSOC leads to a magnetic response very unlike that of
$j=1/2$ superconductors.

General theory. We consider a system described by the Bogoliubov-de Gennes
Hamiltonian

$H=\sum_{\bm{k}}\vec{\Psi}_{\bm{k}}^{\dagger}\begin{pmatrix}\mathcal{H}_{0,\bm{k}}&\Delta_{\bm{k}}\\\
\Delta_{\bm{k}}^{\dagger}&-\mathcal{H}_{0,-\bm{k}}^{T}\end{pmatrix}\vec{\Psi}_{\bm{k}},$
(1)

where $\vec{\Psi}_{\bm{k}}=(\vec{c}_{\bm{k}},\vec{c}^{\dagger}_{-\bm{k}})$,
with $\vec{c}_{\bm{k}}$ representing a spinor of annihilation operators for
fermions with four internal degrees of freedom. In the presence of both TRS
and IS, the most general form of the Hamiltonian matrix ${\cal H}_{0,\bm{k}}$
is [24, 25]

$\mathcal{H}_{0,\bm{k}}=\varepsilon_{0,\bm{k}}\mathbb{1}+\vec{\varepsilon}_{\bm{k}}\cdot\vec{\gamma}=\varepsilon_{0,\bm{k}}\mathbb{1}+\bar{\mathcal{H}}_{0,\bm{k}},$
(2)

where $\vec{\gamma}=(\gamma^{1},\gamma^{2},\gamma^{3},\gamma^{4},\gamma^{5})$
are the mutually anticommuting Euclidean Dirac matrices with coefficients
$\vec{\varepsilon}_{\bm{k}}=(\varepsilon_{1,\bm{k}},\varepsilon_{2,\bm{k}},\varepsilon_{3,\bm{k}},\varepsilon_{4,\bm{k}},\varepsilon_{5,\bm{k}})$.
The Hamiltonian has doubly-degenerate eigenenergies
$\xi_{\pm,{\bm{k}}}=\varepsilon_{0,\bm{k}}\pm|\vec{\varepsilon}_{\bm{k}}|$,
where
$|\vec{\varepsilon}_{\bm{k}}|=\sqrt{\sum_{l=1}^{5}\varepsilon_{l,\bm{k}}^{2}}$.
In Eq. 2 we have introduced
$\bar{\mathcal{H}}_{0,\bm{k}}=\vec{\epsilon}_{\bm{k}}\cdot\vec{\gamma}$ to
denote the part of the Hamiltonian which depends nontrivially on the internal
degrees of freedom.

Hamiltonians with the form of Eq. 2 describe a diverse range of two-band
systems [26, 27, 15, 28, 14, 29, 18, 24, 10], and so the exact form of the
$\gamma$ matrices depends on the system under consideration [30]. For
spin-$\frac{1}{2}$ systems, we generally construct the $\gamma$ matrices as
Kronecker products of Pauli matrices acting in the orbital and spin spaces. In
this case we choose $\gamma^{1}$ and $\gamma^{2}$ to be trivial in the spin
space, with these contributions to the Hamiltonian describing purely orbital
effects, while the remaining matrices couple the spin and orbital degrees of
freedom, accounting for SOC in the system. As we shall see, this
parameterization is also valid for spin-$\frac{3}{2}$ systems, even though the
underlying orbital and spin degrees of freedom cannot be factorized.

The pairing potential appearing in Eq. 1 is written
$\Delta_{\bm{k}}=\Delta_{0}\tilde{\Delta}_{{\bm{k}}}U_{T}$, where $\Delta_{0}$
is the magnitude, $\tilde{\Delta}_{\bm{k}}$ encodes the dependence on the
momentum and the internal degrees of freedom, and $U_{T}$ is the unitary part
of the time-reversal operator. The general form for even- ($e$) and odd-parity
($o$) states is

$\displaystyle\tilde{\Delta}^{(e)}_{\bm{k}}$
$\displaystyle=\sum_{a=0}^{5}e^{a}_{\bm{k}}\gamma^{a}$ (3)
$\displaystyle\tilde{\Delta}^{(o)}_{\bm{k}}$
$\displaystyle=\sum_{a=1}^{4}\sum_{b>a}o^{ab}_{\bm{k}}i\gamma^{a}\gamma^{b}$
(4)

where $e^{a}_{\bm{k}}$ and $o^{ab}_{\bm{k}}$ are normalized form factors. Note
that only when the internal degrees of freedom transform trivially under
inversion are the functions $e^{a}_{\bm{k}}$ and $o^{ab}_{\bm{k}}$ necessarily
even and odd in momentum, respectively.

Due to the mixing of orbital and spin, the pairing potential in the band basis
typically has both intraband and interband matrix elements. The intraband gap
is particularly important as it is responsible for the Cooper instability.
Assuming that $\Delta_{0}$ is small compared to the band separation, the gap
in band $a$ is given by

$\left|\Delta_{a,\bm{k}}\right|^{2}=\Delta_{0}^{2}\frac{\text{Tr}\\{|\\{\bar{\mathcal{H}}_{0,\bm{k}},\tilde{\Delta}_{\bm{k}}\\}|^{2}\mathcal{P}_{a,{\bm{k}}}\\}}{8\left|\vec{\varepsilon}_{\bm{k}}\right|^{2}},$
(5)

where $\\{\bar{\mathcal{H}}_{0,\bm{k}},\tilde{\Delta}_{\bm{k}}\\}U_{T}=F_{A}$
is the superconducting fitness as defined in Refs. [16, 23], and
$\mathcal{P}_{a,{\bm{k}}}=\frac{1}{2}(\mathbb{1}+a\bar{\mathcal{H}}_{0,\bm{k}}/|\vec{\epsilon}|)$
projects into the $a=\pm$ band. The projection operator is necessary to
account for band-dependence of the intraband pairing, which can arise when
$\tilde{\Delta}_{\bm{k}}\tilde{\Delta}_{\bm{k}}^{\dagger}\not\propto\mathbb{1}$.

To investigate the effect of symmetry-breaking, we introduce the perturbation
Hamiltonian

$\delta
H=\sum_{\bm{k}}\vec{\Psi}_{\bm{k}}^{T}\begin{pmatrix}\mathcal{H}_{h,\bm{k}}&0\\\
0&-\mathcal{H}_{h,-\bm{k}}^{T}\end{pmatrix}\vec{\Psi}_{\bm{k}}.$ (6)

We adopt a general form of $\mathcal{H}_{h,\bm{k}}$

$\mathcal{H}_{h,\bm{k}}=\sum_{\alpha=0}^{5}h_{0\alpha,\bm{k}}\gamma^{\alpha}+\sum_{\alpha=1}^{4}\sum_{\beta>\alpha}h_{\alpha\beta,\bm{k}}i\gamma^{\alpha}\gamma^{\beta},$
(7)

The perturbation lifts the twofold degeneracy of the normal state spectrum:
For sufficiently well-separated bands, the perturbed energies of band $a$ are
$\xi_{a,\pm,{\bm{k}}}\approx\xi_{a,{\bm{k}}}\pm\tilde{g}_{a,{\bm{k}}}h_{\bm{k}}$,
where
$h^{2}_{\bm{k}}=\sum_{\alpha,\beta}h_{\alpha\beta,\bm{k}}^{2}=\text{Tr}\\{|\mathcal{H}_{h,\bm{k}}|^{2}\\}/4$,
and the effective $g$-factor in band $a$ is

$\tilde{g}_{a,{\bm{k}}}^{2}=\frac{\text{Tr}\\{|\\{\bar{\mathcal{H}}_{0,\bm{k}},\mathcal{H}_{h,\bm{k}}\\}|^{2}\mathcal{P}_{a,{\bm{k}}}\\}}{8\left|\vec{\varepsilon}_{\bm{k}}\right|^{2}h^{2}_{\bm{k}}}.$
(8)

Equation 8 resembles the expression for the intraband superconducting gap Eq.
5, and can be similarly interpreted as giving the splitting of the bands due
to the projection of the perturbation onto band $a$. When
$|\mathcal{H}_{h,\bm{k}}|^{2}\propto\mathbb{1}$ the splitting of the energy
spectrum is independent of the band index, i.e. $g_{a,\bm{k}}=g_{\bm{k}}$;
more generally, when $|\mathcal{H}_{h,\bm{k}}|^{2}\not\propto\mathbb{1}$ the
effective $g$-factors are different in each band, which is accounted for by
the projection operator in Eq. 8. For definiteness, in the following we
consider only perturbations which break TRS but preserve IS; the results for
perturbations which break IS but preserve TRS are similar and provided in
appendix A.4.

Pair-breaking. The lifting of the two-fold degeneracy of the band states by
the perturbation generally suppresses the superconductivity. The central
result of our work is that the pair-breaking effects of the perturbation in
band $a$ can be quantified by the field-fitness function

$\tilde{F}_{h,{\bm{k}}}^{(a)}=\frac{\text{Tr}\\{|\\{\\{\bar{\mathcal{H}}_{0,\bm{k}},\tilde{\Delta}_{\bm{k}}\\},\\{\bar{\mathcal{H}}_{0,\bm{k}},\mathcal{H}_{h,\bm{k}}\\}\\}|^{2}\mathcal{P}_{a,{\bm{k}}}\\}}{2\text{Tr}\\{|\\{\bar{\mathcal{H}}_{0,\bm{k}},\mathcal{H}_{h,\bm{k}}\\}|^{2}\mathcal{P}_{a,{\bm{k}}}\\}\text{Tr}\\{|\\{\bar{\mathcal{H}}_{0,\bm{k}},\tilde{\Delta}_{\bm{k}}\\}|^{2}\mathcal{P}_{a,{\bm{k}}}\\}}\,.$
(9)

The field-fitness function ranges in value from zero to one. For
$\tilde{F}^{(a)}_{h,{\bm{k}}}=0$, the states at ${\bm{k}}$ and $-{\bm{k}}$
involved in the intraband pairing remain degenerate, and so there is no pair-
breaking effect. On the other hand, the perturbation is maximally pair-
breaking for $\tilde{F}^{(a)}_{h,{\bm{k}}}=1$, i.e. the states paired by the
intraband pairing potential are split by the perturbation. An intermediate
value $0<\tilde{F}^{(a)}_{h,{\bm{k}}}<1$ indicates that the intraband pairing
potential pairs electrons in a superposition of the perturbed states, and
there will be some pair-breaking effect. Inserting Eq. 3 into Eq. 9, we find
that for any even-parity state the field fitness
$\tilde{F}_{h,{\bm{k}}}^{(a)}=1$, as the numerator can be factored to give the
denominator. Since even-parity superconductors always pair time-reversed
partners within the same band, any TRS-breaking perturbation is maximally
pair-breaking. On the other hand, odd-parity superconducting states do not
necessarily pair time-reversed states in the same band, and so they may
experience less or no pair-breaking due to broken TRS, i.e.
$0\leq\tilde{F}_{h}^{(a)}\leq 1$.

Solving the linearized gap equation in the presence of the TRS-breaking
perturbation gives the critical temperature $T_{c}$ in terms of the
unperturbed value $T_{c,0}$,

$\log\left(\frac{T_{c}}{T_{c,0}}\right)=\sum_{a=\pm}\left\langle\left[\frac{\mathcal{D}_{a,\bm{k}}\left|\Delta_{a,\bm{k}}\right|^{2}}{\sum_{a^{\prime}=\pm}\langle\mathcal{D}_{a^{\prime},\bm{k}^{\prime}}\left|\Delta_{a^{\prime},\bm{k}^{\prime}}\right|^{2}\rangle_{a^{\prime}}}\right]\tilde{F}_{h,{\bm{k}}}^{(a)}\text{Re}\left\\{\psi\left(\frac{1}{2}\right)-\psi\left(\frac{1}{2}+i\frac{\tilde{g}_{a,{\bm{k}}}h_{\bm{k}}}{2\pi
k_{B}T_{c}}\right)\right\\}\right\rangle_{a},$ (10)

where $\psi(x)$ is the digamma function, $\langle\ldots\rangle_{a}$ indicates
the average over the Fermi surface of band $a$,
$\mathcal{D}_{a,\bm{k}}=|\vec{\nabla}_{\bm{k}}\xi_{a,\bm{k}}|^{-1}$, and the
factor in the square brackets defines the fraction of the total condensation
energy due to the gap on each band. The suppression of the critical
temperature by a TRS-breaking perturbation is controlled by both the field
fitness function and the effective $g$-factor, which tune the degree of pair-
breaking and the magnitude of the band splitting, respectively. A brief
derivation of Eq. 10 is presented in appendix A.

Magnetic susceptibility. We now turn to the important case where the
perturbation is an applied magnetic field, which couples to the electron
states via the Zeeman effect. A key experimental quantity is the magnetic
susceptibility, which in a multiband system can be divided into components due
to intraband (‘Pauli’) and interband (‘van Vleck’) transitions. The latter is
negligibly affected by superconductivity, as the pairing potential is
typically much smaller than the band separation. On the other hand, the Pauli
contribution carries clear signatures of the pair-breaking effect. For a field
applied along the $i$-axis, the Pauli susceptibility below the critical
temperature is given by

$\chi_{ii}=\sum_{a=\pm}\left\langle
2\mu_{B}^{2}\mathcal{D}_{a,\bm{k}}\tilde{g}^{(i)\,2}_{a,\bm{k}}\left\\{1+\tilde{F}_{i,{\bm{k}}}^{(a)}\left[Y_{a}(\hat{k},T)-1\right]\right\\}\right\rangle_{a},$
(11)

where $\tilde{g}^{(i)}_{a,\bm{k}}$ is the effective $g$-factor for an $i$-axis
field, and $Y_{a}(\hat{k},T)$ is the angle-dependent Yosida function for the
intraband gap Eq. 5. Explicit expressions for the susceptibility in both the
normal and superconducting states are provided in the appendix A. For even-
parity gaps $\tilde{F}_{h,{\bm{k}}}^{(a)}=1$, and the pairing suppresses the
Pauli susceptibility, with it vanishing at zero temperature. On the other
hand, an odd-parity state typically gives only a partial suppression of the
Pauli susceptibility; In the extreme case $\tilde{F}_{h,{\bm{k}}}^{(a)}=0$ the
susceptibility is unaffected by the superconductivity.

$j=\frac{3}{2}$ superconductors. As a concrete application of our approach we
consider a system of electrons with an effective $j=\frac{3}{2}$. A minimal
model is given by the Luttinger-Kohn Hamiltonian

$\mathcal{H}_{0,\bm{k}}=\left(\alpha\left|\bm{k}\right|^{2}-\mu\right)\mathbb{1}+\beta_{1}\sum_{i}k_{i}^{2}J_{i}+\beta_{2}\sum_{i^{\prime}\neq
i}k_{i}k_{i^{\prime}}J_{i}J_{i^{\prime}},$ (12)

where $J_{i}$ are the spin-$\frac{3}{2}$ matrices and the indices $i$,
$i^{\prime}$ run over Cartesian coordinates. The Hamiltonian can be cast into
the form of Eq. 2 by defining the $\gamma$ matrices
$\vec{\gamma}=(\frac{1}{\sqrt{3}}[J_{x}^{2}-J_{y}^{2}],\frac{1}{3}[2J_{z}^{2}-J_{x}^{2}-J_{y}^{2}],\frac{1}{\sqrt{3}}\\{J_{y},J_{z}\\},\frac{1}{\sqrt{3}}\\{J_{x},J_{y}\\},\frac{1}{\sqrt{3}}\\{J_{x},J_{z}\\})$,
with corresponding coefficients
$\vec{\varepsilon}=(\sqrt{3}\beta_{1}(k_{x}^{2}-k_{y}^{2})/2,\sqrt{3}\beta_{1}(3k_{z}^{2}-|\bm{k}|^{2})/2,\sqrt{3}\beta_{2}k_{y}k_{z},\sqrt{3}\beta_{2}k_{x}k_{y},\sqrt{3}\beta_{2}k_{x}k_{z})$.
The strong SOC in the Luttinger-Kohn model leads to a pronounced spin-momentum
locking. This is most clearly observed in the limit $\beta_{1}=\beta_{2}$,
where the model has full spherical symmetry and the quantity
${\bm{k}}\cdot{\bf J}$ commutes with the Hamiltonian, i.e. the projection of
the spin along ${\bm{k}}$ is a good quantum number. The eigenstates are
classified by the “helical” index $\sigma$ such that $\hat{\bm{k}}\cdot{\bf
J}|\sigma\rangle_{\bm{k}}=\sigma|\sigma\rangle_{\bm{k}}$, defining a
spin-$\frac{3}{2}$ band ($\sigma=\pm\frac{3}{2}$) and spin-$\frac{1}{2}$ band
($\sigma=\pm\frac{1}{2}$). This argument remains valid along high-symmetry
directions in the presence of cubic anisotropy, allowing us to identify
spin-$\frac{1}{2}$ and -$\frac{3}{2}$ bands when $\beta_{1}\neq\beta_{2}$. As
we will show below, compared to the spin-$\frac{1}{2}$-band, the
spin-$\frac{3}{2}$ band exhibits counter-intuitive properties.

The spin-momentum locking produces a highly anisotropic Zeeman splitting of
the bands by an applied magnetic field $H_{h}=g\mu_{B}{\bm{h}}\cdot{\bf J}$.
Although the general result using Eq. 8 is complicated, in the rotationally-
symmetric limit the effective $g$-factors take the compact forms
$\tilde{g}_{1/2}=\sqrt{1-3|\hat{\bm{h}}\cdot\hat{\bm{k}}|^{2}/4}$ and
$\tilde{g}_{3/2}=3|\hat{\bm{h}}\cdot\hat{\bm{k}}|/2$ in the spin-$\frac{1}{2}$
and -$\frac{3}{2}$ bands, respectively; these effective $g$-factors remain
approximately valid in the general case. Note that there is no splitting of
the spin-$\frac{3}{2}$ band in the direction perpendicular to the field,
whereas the splitting of the spin-$\frac{1}{2}$ states is maximal in this
plane. To understand this, we introduce the momentum-dependent angular
momentum operators $J_{z}^{\bm{k}}=\hat{\bm{k}}\cdot{\bf J}$, and raising and
lowering operators $J_{\pm}^{\bm{k}}$ which satisfy
$[J_{z}^{\bm{k}},{J}_{\pm}^{\bm{k}}]=\pm{J}_{\pm}^{\bm{k}}$. The band states
are eigenstates of $J_{z}^{\bm{k}}$, i.e.
$J_{z}^{\bm{k}}|\sigma\rangle_{\bm{k}}=\sigma|\sigma\rangle_{\bm{k}}$, and the
action of the raising and lowering operators is
$J_{\pm}^{\bm{k}}|\sigma\rangle_{\bm{k}}=\sqrt{(\frac{3}{2}\mp\sigma)(\frac{3}{2}\pm\sigma+1)}|\sigma\pm
1\rangle_{\bm{k}}$. Expressing the Zeeman Hamiltonian in terms of these
operators we have $H_{h}=g\mu_{B}\left({\bf
h}\cdot\hat{\bm{k}}J_{z}^{\bm{k}}+h_{+}^{\bm{k}}J_{+}^{\bm{k}}+h_{-}^{\bm{k}}J_{-}^{\bm{k}}\right)$;
the functions $h_{\pm}^{\bm{k}}=-h_{\pm}^{-{\bm{k}}}$ are vanishing for
$\hat{\bm{k}}\parallel{\bf h}$ and take maximal magnitude in the plane
perpendicular to ${\bm{k}}$; explicit expressions are given in the appendix B.
Since the raising and lowering operators do not couple the
$\sigma=\pm\frac{3}{2}$ states, this immediately explains why there is no
splitting for ${\bm{k}}\perp{\bf h}$ in the spin-$\frac{3}{2}$ band;
conversely, the larger matrix elements for the raising and lowering operators
compared to the $z$-spin operator gives a maximum of the Zeeman splitting for
the spin-$\frac{1}{2}$ band in this plane.

The spin-momentum locking also profoundly influences the superconducting gap
structure, which we illustrate by two odd-parity states: the fully-gapped
$p$-wave triplet $A_{1u}$ state with
$\tilde{\Delta}_{1u}=\frac{2}{\sqrt{5}}({\bm{k}}\cdot{\bf J})U_{T}$, and a
nodal $p$-wave septet $A_{2u}$ state
$\tilde{\Delta}_{2u}=\frac{1}{\sqrt{3}}({\bm{k}}\cdot{\bf T})U_{T}$, where
$T_{i}=\\{J_{i},J_{i+1}^{2}-J_{i+2}^{2}\\}$ for $i$ defined cyclically (e.g.
$J_{y+2}=J_{x}$). These two states are the leading $p$-wave instabilities in
the spin-$\frac{3}{2}$ band in the rotationally-symmetric limit [31]. Using
Eq. 5 we find that the gaps opened by these pairing states have nontrivial
band dependence. This is most apparent in the $A_{2u}$ case, which in the
spin-$\frac{1}{2}$ band has fourteen nodes along the $(100)$ and $(111)$ and
equivalent directions, whereas in the spin-$\frac{3}{2}$ band it only has six
nodes along the the $(100)$ directions [32]. These nodes result from the spin-
momentum locking: along the $(100)$ directions the $A_{2u}$ state pairs
electrons with helicity differing by $\pm 2$ which cannot be satisfied in
either band, and hence implies purely interband pairing and a node in the
intraband gap. On the other hand, along the $(111)$ direction the $A_{2u}$
state pairs electrons with helicity differing by $\pm 3$, and so it only opens
a gap in the spin-$\frac{3}{2}$ band. In contrast, the $A_{1u}$ state pairs
electrons with the same helicity, and hence it opens a full gap in both bands.
Further details on these pairing states are provided in the appendix B.

Figure 1: Comparison of $A_{1u}$ and $A_{2u}$ states. (a), (b): Normalized
field-fitness $\tilde{F}_{h}^{(a)}$ on the spin-$3/2$ (black solid line,
$\mu=-20$meV) and spin-$1/2$ (blue dashed line, $\mu=20$meV) band Fermi
surface for $A_{1u}$ (a) and $A_{2u}$ (b) pairing. (c), (d): Pauli
susceptibility as a function of temperature at $\beta_{2}=\beta_{1}$ for
$A_{1u}$ (c) and $A_{2u}$ (d) pairing. (e), (f) Upper critical field excluding
orbital effects for each band at $\beta_{2}=\beta_{1}$ for $A_{1u}$ (e) and
$A_{2u}$ (f) pairing. In all plots we take $\alpha=20(a/\pi)^{2}$eV,
$\beta_{1}=-15(a/\pi)^{2}$eV.

The response of the superconductivity to an applied Zeeman field displays
major differences between the two bands. For a field applied along the
$z$-axis, the averaged field-fitness function
$\tilde{F}_{h}^{(a)}=\langle\tilde{F}_{h,{\bm{k}}}^{(a)}\rangle_{a}$ is
plotted in Fig. 1(a) and (b) as a function of $\beta_{2}/\beta_{1}$ for the
$A_{1u}$ and $A_{2u}$ states, respectively. In the spin-$\frac{1}{2}$ band,
the nonzero value of $\tilde{F}_{h}^{(a)}$ indicates that a magnetic field is
partially pair-breaking for both states; in the spin-$\frac{3}{2}$ band,
however, the field is almost completely pair-breaking for the $A_{1u}$ state,
whereas the $A_{2u}$ state is almost insensitive to its effects. These results
are only weakly dependent upon the ratio $\beta_{2}/\beta_{1}$, and become
exact in the rotationally-symmetric limit. This behaviour is reflected in the
susceptibility, as shown in Fig. 1(c) and (d): in the spin-$\frac{3}{2}$ band
the Pauli contribution to the susceptibility is completely suppressed for the
$A_{1u}$ state, whereas we see no change in the susceptibility for the
$A_{2u}$ state. In contrast, in both cases there is moderate suppression of
the susceptibility in the spin-$\frac{1}{2}$ band. Likewise, as shown in Fig.
1(e) and (f), in the spin-$\frac{3}{2}$ band the $A_{1u}$ state is Pauli-
limited whereas the $A_{2u}$ state is only limited by orbital effects.

The remarkable response to a magnetic field can be understood in the
rotationally-symmetric limit in terms of the helicity quantum number.
Projected onto the spin-$\frac{3}{2}$ band, the $A_{1u}$ state pairs electrons
with the same helicity whereas the $A_{2u}$ state pairs electrons with
opposite helicity. In the projection of the Zeeman Hamiltonian only the
$J^{z}_{\bm{k}}$ operator has nonzero matrix elements, and so the energy shift
of the state $|\sigma\rangle_{\bm{k}}$ is $\sigma{\bf h}\cdot\hat{\bm{k}}$.
Since the $A_{2u}$ state pairs electrons which have the same Zeeman shift, it
does not experience pair-breaking; in contrast, the $A_{1u}$ state pairs
electrons with opposite Zeeman shift, and the pair-breaking is maximal. The
situation in the spin-$\frac{1}{2}$ band is more complicated, since helicity
is not a good quantum number for the projected Zeeman Hamiltonian, and the
$A_{2u}$ state has both equal- and opposite-helicity pairing matrix elements.
This gives the intermediate values of $\tilde{F}_{h}^{(a)}$ for the
spin-$\frac{1}{2}$ band. Although this argument is only rigorously valid in
the rotationally-symmetric limit, the plot of $\tilde{F}_{h}^{(a)}$ in Fig.
1(a) and (b) indicates that it remains valid more generally.

Conclusions. In this paper we have developed a basis-independent framework to
understand the interplay of superconductivity and symmetry-breaking
perturbations in a multiband system. Using a generic minimal model, we have
shown that that coupling of orbital and spin degrees of freedom typically
reduces the effect of the perturbation via the appearance of an effective
$g$-factor Eq. 8. Moreover, the pair-breaking effect of this perturbation can
be formulated in terms of the field-fitness function Eq. 9. Together, these
quantities control the suppression of the critical temperature by the
perturbation, and in the case of a Zeeman field determine the magnetic
susceptibility below the critical temperature. To illustrate these effects, we
have examined $p$-wave pairing of effective spin-$\frac{3}{2}$ electrons with
normal state described by the Luttinger-Kohn model. The characteristic spin-
momentum locking in this model leads to the remarkable result that in the
spin-$\frac{3}{2}$ band the triplet $A_{1u}$ state is strongly suppressed by a
Zeeman field, whereas the septet $A_{2u}$ state is largely immune to it,
giving dramatically different behaviour in the upper critical field and the
spin susceptibility.

Acknowledgements D.C.C. and P.M.R.B were supported by the Marsden Fund Council
from Government funding, managed by Royal Society Te Apārangi. D.F.A. was
supported by the U.S. Department of Energy, Office of Basic Energy Sciences,
Division of Materials Sciences and Engineering, under Award No. DE-SC0021971.

## References

  * Ran  _et al._ [2019] S. Ran, C. Eckberg, Q.-P. Ding, Y. Furukawa, T. Metz, S. R. Saha, I.-L. Liu, M. Zic, H. Kim, J. Paglione, and N. P. Butch, Nearly ferromagnetic spin-triplet superconductivity, Science 365, 684 (2019).
  * Khim  _et al._ [2021] S. Khim, J. F. Landaeta, J. Banda, N. Bannor, M. Brando, P. M. R. Brydon, D. Hafner, R. Küchler, R. Cardoso-Gil, U. Stockert, A. P. Mackenzie, D. F. Agterberg, C. Geibel, and E. Hassinger, Field-induced transition within the superconducting state of CeRh2As2, Science 373, 1012 (2021).
  * Aoki  _et al._ [2019] D. Aoki, A. Nakamura, F. Honda, D. Li, Y. Homma, Y. Shimizu, Y. J. Sato, G. Knebel, J.-P. Brison, A. Pourret, D. Braithwaite, G. Lapertot, Q. Niu, M. Vališka, H. Harima, and J. Flouquet, Unconventional superconductivity in heavy fermion UTe2, Journal of the Physical Society of Japan 88, 043702 (2019).
  * Nguyen  _et al._ [2021] D. H. Nguyen, A. Sidorenko, M. Taupin, G. Knebel, G. Lapertot, E. Schuberth, and S. Paschen, Superconductivity in an extreme strange metal, Nature Communications 12, 4341 (2021).
  * Mineev and Samokhin [1999] V. P. Mineev and K. V. Samokhin,  _Introduction to Unconventional Superconductivity_ (Gordon and Breach Science Publishers, 1999).
  * Sigrist [2005] M. Sigrist, Introduction to unconventional superconductivity, AIP Conference Proceedings 789, 165 (2005).
  * Shimozawa  _et al._ [2014] M. Shimozawa, S. K. Goh, R. Endo, R. Kobayashi, T. Watashige, Y. Mizukami, H. Ikeda, H. Shishido, Y. Yanase, T. Terashima, T. Shibauchi, and Y. Matsuda, Controllable rashba spin-orbit interaction in artificially engineered superlattices involving the heavy-fermion superconductor ${\mathrm{cecoin}}_{5}$, Phys. Rev. Lett. 112, 156404 (2014).
  * Watanabe  _et al._ [2015] T. Watanabe, T. Yoshida, and Y. Yanase, Odd-parity superconductivity by competing spin-orbit coupling and orbital effect in artificial heterostructures, Phys. Rev. B 92, 174502 (2015).
  * Cavanagh  _et al._ [2022] D. C. Cavanagh, T. Shishidou, M. Weinert, P. M. R. Brydon, and D. F. Agterberg, Nonsymmorphic symmetry and field-driven odd-parity pairing in $\mathrm{Ce}{\mathrm{rh}}_{2}{\mathrm{as}}_{2}$, Phys. Rev. B 105, L020505 (2022).
  * Xie  _et al._ [2020] Y.-M. Xie, B. T. Zhou, and K. T. Law, Spin-orbit-parity-coupled superconductivity in topological monolayer ${\mathrm{WTe}}_{2}$, Phys. Rev. Lett. 125, 107001 (2020).
  * Zhang  _et al._ [2014] X. Zhang, Q. Liu, J.-W. Luo, A. J. Freeman, and A. Zunger, Hidden spin polarization in inversion-symmetric bulk crystals, Nature Physics 10, 387 (2014).
  * Smidman  _et al._ [2017] M. Smidman, M. B. Salamon, H. Q. Yuan, and D. F. Agterberg, Superconductivity and spin–orbit coupling in non-centrosymmetric materials: a review, Reports on Progress in Physics 80, 036501 (2017).
  * Youn  _et al._ [2012] S. J. Youn, M. H. Fischer, S. H. Rhim, M. Sigrist, and D. F. Agterberg, Role of strong spin-orbit coupling in the superconductivity of the hexagonal pnictide srptas, Phys. Rev. B 85, 220505 (2012).
  * Yoshida  _et al._ [2012] T. Yoshida, M. Sigrist, and Y. Yanase, Pair-density wave states through spin-orbit coupling in multilayer superconductors, Phys. Rev. B 86, 134514 (2012).
  * Vafek and Chubukov [2017] O. Vafek and A. V. Chubukov, Hund Interaction, Spin-Orbit Coupling, and the Mechanism of Superconductivity in Strongly Hole-Doped Iron Pnictides, Phys. Rev. Lett. 118, 087003 (2017).
  * Ramires and Sigrist [2016] A. Ramires and M. Sigrist, Identifying detrimental effects for multiorbital superconductivity: Application to ${\mathrm{Sr}}_{2}{\mathrm{RuO}}_{4}$, Phys. Rev. B 94, 104501 (2016).
  * Suh  _et al._ [2020] H. G. Suh, H. Menke, P. M. R. Brydon, C. Timm, A. Ramires, and D. F. Agterberg, Stabilizing even-parity chiral superconductivity in ${\mathrm{sr}}_{2}{\mathrm{ruo}}_{4}$, Phys. Rev. Research 2, 032023 (2020).
  * Brydon  _et al._ [2016] P. M. R. Brydon, L. Wang, M. Weinert, and D. F. Agterberg, Pairing of $j=3/2$ Fermions in Half-Heusler Superconductors, Phys. Rev. Lett. 116, 177001 (2016).
  * Savary  _et al._ [2017] L. Savary, J. Ruhman, J. W. F. Venderbos, L. Fu, and P. A. Lee, Superconductivity in three-dimensional spin-orbit coupled semimetals, Phys. Rev. B 96, 214514 (2017).
  * Kobayashi  _et al._ [2022] S. Kobayashi, A. Bhattacharya, C. Timm, and P. M. R. Brydon, Bogoliubov fermi surfaces from pairing of emergent $j=\frac{3}{2}$ fermions on the pyrochlore lattice, Phys. Rev. B 105, 134507 (2022).
  * Venderbos  _et al._ [2018a] J. W. F. Venderbos, L. Savary, J. Ruhman, P. A. Lee, and L. Fu, Pairing States of Spin-$\frac{3}{2}$ Fermions: Symmetry-Enforced Topological Gap Functions, Phys. Rev. X 8, 011029 (2018a).
  * Kim  _et al._ [2022] D. Kim, T. Sato, S. Kobayashi, and Y. Asano, Spin susceptibility of a $j=3/2$ superconductor, arXiv:2206.14994 (2022).
  * Ramires  _et al._ [2018] A. Ramires, D. F. Agterberg, and M. Sigrist, Tailoring ${T}_{c}$ by symmetry principles: The concept of superconducting fitness, Phys. Rev. B 98, 024501 (2018).
  * Brydon  _et al._ [2018] P. M. R. Brydon, D. F. Agterberg, H. Menke, and C. Timm, Bogoliubov Fermi surfaces: General theory, magnetic order, and topology, Phys. Rev. B 98, 224509 (2018).
  * Abrikosov [1974] A. A. Abrikosov, Calculation of critical indices for zero-gap semiconductors, Sov. Phys. JETP 39, 709 (1974).
  * Fu and Berg [2010] L. Fu and E. Berg, Odd-Parity Topological Superconductors: Theory and Application to ${\mathrm{Cu}}_{x}{\mathrm{Bi}}_{2}{\mathrm{Se}}_{3}$, Phys. Rev. Lett. 105, 097001 (2010).
  * Yanase [2016] Y. Yanase, Nonsymmorphic Weyl superconductivity in UPt3 based on ${E}_{2u}$ representation, Phys. Rev. B 94, 174502 (2016).
  * Agterberg  _et al._ [2017] D. F. Agterberg, T. Shishidou, J. O’Halloran, P. M. R. Brydon, and M. Weinert, Resilient nodeless $d$-wave superconductivity in monolayer fese, Phys. Rev. Lett. 119, 267001 (2017).
  * Ilić  _et al._ [2017] S. Ilić, J. S. Meyer, and M. Houzet, Enhancement of the upper critical field in disordered transition metal dichalcogenide monolayers, Phys. Rev. Lett. 119, 117001 (2017).
  * Denys and Brydon [2021] M. D. E. Denys and P. M. R. Brydon, Origin of the anomalous hall effect in two-band chiral superconductors, Phys. Rev. B 103, 094503 (2021).
  * Link and Herbut [2022] J. M. Link and I. F. Herbut, $p$-wave superconductivity in Luttinger semimetals, Phys. Rev. B 105, 134522 (2022).
  * Venderbos  _et al._ [2018b] J. W. F. Venderbos, L. Savary, J. Ruhman, P. A. Lee, and L. Fu, Pairing States of Spin-$\frac{3}{2}$ Fermions: Symmetry-Enforced Topological Gap Functions, Phys. Rev. X 8, 011029 (2018b).

##  Appendix A General results on the effect of broken symmetry

###  A.1 The Green’s function

In the absence of symmetry-breaking perturbations, the Green’s function is

$G\left(\bm{k},i\omega_{n}\right)=\sum_{a=\pm}\frac{1}{i\omega_{n}-\xi_{a,\bm{k}}}\mathcal{P}_{a,\bm{k}},$
(13)

where $\mathcal{P}_{a,\bm{k}}$ projects onto band $a$ at wavevector
${\bm{k}}$, given explicitly by

$\mathcal{P}_{a,\bm{k}}=\frac{\mathbb{1}+a\hat{\varepsilon}_{\bm{k}}\cdot\vec{\gamma}}{2}=\frac{1}{2}\left[\mathbb{1}+a\frac{\bar{\mathcal{H}}_{0,\bm{k}}}{\left|\vec{\varepsilon}_{\bm{k}}\right|}\right].$
(14)

We now consider the addition of the symmetry-breaking perturbation. Assuming
that $h_{\bm{k}}\ll|\vec{\epsilon}_{\bm{k}}|$, i.e. the mixing between the
unperturbed bands is negligible, we can approximate the Green’s function by

$G\left(\bm{k},i\omega_{n}\right)=\sum_{a,b=\pm}\frac{1}{i\omega_{n}-\xi_{a,b,\bm{k}}}\mathcal{P}_{a,b,\bm{k}},$
(15)

where $\xi_{a,b,{\bm{k}}}=\xi_{a,{\bm{k}}}+b\tilde{g}_{a,{\bm{k}}}h_{\bm{k}}$
and the projection operator is defined

$\mathcal{P}_{a,b,\bm{k}}=\frac{1}{4}\left(\mathbb{1}+a\frac{\bar{\mathcal{H}}_{0,\bm{k}}}{\left|\vec{\varepsilon}_{\bm{k}}\right|}\right)\left(\mathbb{1}+ab\frac{\left\\{\bar{\mathcal{H}}_{0,\bm{k}},\mathcal{H}_{h,\bm{k}}\right\\}}{2\tilde{g}_{a,\bm{k}}h\left|\vec{\varepsilon}_{\bm{k}}\right|}\right)$
(16)

This expression is valid to lowest order in
$h_{\alpha,\bm{k}}/|\vec{\varepsilon}_{\bm{k}}|$, which is equivalent to
ignoring interband corrections.

###  A.2 The linearized gap equation

The pairing potential $\Delta_{0}$ is determined by self-consistent solution
of the gap equation

$\Delta_{0}=-V_{0}\beta\sum_{i\omega_{n},\bm{k}}\text{Tr}\\{\tilde{\Delta}^{\dagger}_{\bm{k}}F\left(\bm{k},i\omega_{n}\right)\\},$
(17)

where $F({\bm{k}},i\omega_{n})$ is the anomalous Green’s function and $V_{0}$
is the pairing interaction. Ignoring interband pairing, just below the
critical temperature the anomalous Green’s function in the presence of
symmetry-breaking perturbations can be approximated by

$F\left(\bm{k},i\omega_{n}\right)=\sum_{a,b,b^{\prime}}\frac{\mathcal{P}_{a,b,\bm{k}}\Delta_{\bm{k}}\mathcal{P}_{a,b^{\prime},-\bm{k}}^{T}}{\left(i\omega_{n}-\xi_{a,b,\bm{k}}\right)\left(i\omega_{n}+\xi_{a,b^{\prime},\bm{k}}\right)}\,.$
(18)

Inserting this into the gap equation we obtained the linearized form

$-\frac{1}{V_{0}}=\beta\sum_{i\omega_{n},\bm{k}}\sum_{a,b,b^{\prime}}\frac{\text{Tr}\\{\tilde{\Delta}_{\bm{k}}^{\dagger}\mathcal{P}_{a,b,\bm{k}}\tilde{\Delta}_{\bm{k}}\mathcal{P}_{a,\bar{b}^{\prime},\bm{k}}\\}}{\left(i\omega_{n}-\xi_{a,b,\bm{k}}\right)\left(i\omega_{n}+\xi_{a,b^{\prime},\bm{k}}\right)}\,.$
(19)

Performing the sums over Matsubara frequency and momentum gives Eq. 10 after
the zero-field result is subtracted.

###  A.3 Spin susceptibility

The static long-wavelength paramagnetic $i$-axis susceptibility can be divided
into an intraband Pauli contribution and an interband van Vleck contribution

$\chi_{ii}=\sum_{a}\\{\chi_{ii,aa}+\chi_{ii,a\bar{a}}\\}\,.$ (20)

In the normal state these terms are given by

$\displaystyle\chi_{ii,aa}^{(N)}$
$\displaystyle=\sum_{\bm{k}}\frac{\mu_{B}^{2}\bar{\chi}_{ii,aa}^{(N)}\left(\bm{k}\right)}{8\left|\vec{\varepsilon}_{\bm{k}}\right|^{2}}\text{Tr}\\{|\left\\{{S}_{i},\bar{\mathcal{H}}_{0,\bm{k}}\\}|^{2}\mathcal{P}_{a,\bm{k}}\right\\}$
(21) $\displaystyle\chi_{ii,a\bar{a}}^{(N)}$
$\displaystyle=\sum_{\bm{k}}\frac{\mu_{B}^{2}\bar{\chi}_{ii,a\bar{a}}^{(N)}\left(\bm{k}\right)}{16\left|\vec{\varepsilon}_{\bm{k}}\right|^{2}}\text{Tr}\\{|[{S}_{i},\bar{\mathcal{H}}_{0,\bm{k}}]|^{2}\\}.$
(22)

Here ${S}_{i}$ is the spin operator in the $i$ direction. In
spin-$\frac{1}{2}$ materials this is given by
${S}_{i}=\tau_{0}\otimes\sigma_{i}$, whereas ${S}_{i}=J_{i}$ for effective
spin-$\frac{3}{2}$ materials. The traces in these expression encode the effect
of the spin-orbital texture of the band electron states; note the appearance
of the effective $g$-factor in the Pauli contribution. This modulates the
Lindhard function

$\displaystyle\bar{\chi}_{ii,aa^{\prime}}^{(N)}\left(\bm{k}\right)$
$\displaystyle=\lim\limits_{\bm{q}\rightarrow\bm{0}}\frac{n_{F}\left(\xi_{a,\bm{k}+\bm{q}}\right)-n_{F}\left(\xi_{a^{\prime},\bm{k}}\right)}{\xi_{a,\bm{k}+\bm{q}}-\xi_{a^{\prime},\bm{k}}}$
(23)

where $n_{F}$ is the Fermi-Dirac function.

The Pauli susceptibility is significantly modified below the critical
temperature due to the opening of the superconducting gap. Moreover, it
develops an anomalous contribution due to the pairing

$\displaystyle\chi_{ii,aa}^{(A)}$
$\displaystyle=\sum_{\bm{k}}\frac{\mu_{B}^{2}\bar{\chi}_{ii,aa}^{(A)}\left(\bm{k}\right)}{256\left|\vec{\varepsilon}_{\bm{k}}\right|^{4}}\text{Tr}\\{(|\\{\\{{S}_{i},\bar{\mathcal{H}}_{0,\bm{k}}\\},\\{\tilde{\Delta}_{\bm{k}},\bar{\mathcal{H}}_{0,\bm{k}}\\}\\}|^{2}-|[\\{{S}_{i},\bar{\mathcal{H}}_{0,\bm{k}}\\},\\{\tilde{\Delta}_{\bm{k}},\bar{\mathcal{H}}_{0,\bm{k}}\\}]|^{2})\mathcal{P}_{a,\bm{k}}\\}.$
(24)

The modified Lindhard function appearing in this expression is

$\displaystyle\bar{\chi}_{ii,aa^{\prime}}^{(A)}\left(\bm{k}\right)$
$\displaystyle=\lim\limits_{\bm{q}\rightarrow\bm{0}}\frac{\left|\Delta_{0,\bm{k}+\bm{q}}\right|\left|\Delta_{0,\bm{k}}\right|}{E_{a,\bm{k}+\bm{q}}E_{a^{\prime},\bm{k}}}\left\\{\frac{n_{F}\left(E_{a,\bm{k}+\bm{q}}\right)-n_{F}\left(E_{a^{\prime},\bm{k}}\right)}{E_{a,\bm{k}+\bm{q}}-E_{a^{\prime},\bm{k}}}-\frac{1-n_{F}\left(E_{a,\bm{k}+\bm{q}}\right)-n_{F}\left(E_{a^{\prime},\bm{k}}\right)}{E_{a,\bm{k}+\bm{q}}+E_{a^{\prime},\bm{k}}}\right\\},$
(25)

where $E_{a,{\bm{k}}}=\sqrt{\xi_{a,{\bm{k}}}^{2}+|\Delta_{a,{\bm{k}}}|^{2}}$.
The Lindhard function appearing in the “normal” contribution is also modified
as

$\displaystyle\bar{\chi}_{ii,aa^{\prime}}^{(N)}\left(\bm{k}\right)$
$\displaystyle=\lim\limits_{\bm{q}\rightarrow\bm{0}}\frac{-\mu_{B}}{2}\left\\{\left[1+\frac{\xi_{a,\bm{k}+\bm{q}}\xi_{a^{\prime},\bm{k}}}{E_{a,\bm{k}+\bm{q}}E_{a^{\prime},\bm{k}}}\right]\frac{n_{F}\left(E_{a,\bm{k}+\bm{q}}\right)-n_{F}\left(E_{a^{\prime},\bm{k}}\right)}{E_{a,\bm{k}+\bm{q}}-E_{a^{\prime},\bm{k}}}\right.$
$\displaystyle\qquad\qquad\qquad\left.+\left[1-\frac{\xi_{a,\bm{k}+\bm{q}}\xi_{a^{\prime},\bm{k}}}{E_{a,\bm{k}+\bm{q}}E_{a^{\prime},\bm{k}}}\right]\frac{1-n_{F}\left(E_{a,\bm{k}+\bm{q}}\right)-n_{F}\left(E_{a^{\prime},\bm{k}}\right)}{E_{a,\bm{k}+\bm{q}}+E_{a^{\prime},\bm{k}}}\right\\},$
(26)

If the superconducting gap has even parity, the anomalous contribution to the
Pauli susceptibility will exactly cancel the normal part as $T\rightarrow 0$.
More generally any such cancellation depends sensitively on the spin-orbital
texture and structure of the gap.

###  A.4 Broken inversion symmetry

In the main text we focus on the effect of a perturbation which breaks time-
reversal symmetry. Our formalism, can also account for inversion symmetry-
breaking perturbations such as a weak antisymmetric spin-orbit coupling. In
the presence of an inversion symmetry breaking perturbation, the critical
temperature is given by Eq. 10, with the field-fitness replaced by

$\tilde{F}_{h,{\bm{k}}}^{(a)}=\frac{\text{Tr}\\{|[\\{\bar{\mathcal{H}}_{0,\bm{k}},\tilde{\Delta}_{\bm{k}}\\},\\{\bar{\mathcal{H}}_{0,\bm{k}},\mathcal{H}_{h,\bm{k}}\\}]|^{2}\mathcal{P}_{a,{\bm{k}}}\\}}{2\text{Tr}\\{|\\{\bar{\mathcal{H}}_{0,\bm{k}},\mathcal{H}_{h,\bm{k}}\\}|^{2}\mathcal{P}_{a,{\bm{k}}}\\}\text{Tr}\\{|\\{\bar{\mathcal{H}}_{0,\bm{k}},\tilde{\Delta}_{\bm{k}}\\}|^{2}\mathcal{P}_{a,{\bm{k}}}\\}}.$
(27)

The key difference compared to Eq. 9 of the main text is that the numerator
involves the commutator of the anticommutators, as opposed to the
anticommutator of the anticommutators. It can be straightforwardly shown that
$\tilde{F}_{h,{\bm{k}}}^{(a)}=0$ for all even-parity superconductors,
reflecting the fact that the intraband potential of these states always pairs
time-reversed partners.

##  Appendix B $J=\frac{3}{2}$ superconductors

###  B.1 Spherical limit

In the spherical limit the Luttinger-Kohn Hamiltonian takes the form

$H=\alpha|{\bm{k}}|^{2}\hat{1}_{4}+\beta({\bm{k}}\cdot{\bf J})^{2}\,.$ (28)

We observe that the Hamiltonian has two doubly-degenerate bands corresponding
to $\langle{\bm{k}}\cdot{\bf J}\rangle=\pm\frac{3}{2}|{\bm{k}}|$ and
$\langle{\bm{k}}\cdot{\bf J}\rangle=\pm\frac{1}{2}|{\bm{k}}|$. We will refer
to these two bands as the spin-$\frac{3}{2}$ and spin-$\frac{1}{2}$ bands,
respectively. A natural choice of basis for these bands is the “helical” basis

$\\{|\tfrac{3}{2}\rangle_{\bm{k}}\,,\,\,|-\tfrac{3}{2}\rangle_{\bm{k}}\\}\,,\qquad\\{|\tfrac{1}{2}\rangle_{\bm{k}}\,,\,\,|-\tfrac{1}{2}\rangle_{\bm{k}}\\}$
(29)

where ${\bm{k}}\cdot{\bf
J}|\sigma\rangle_{\bm{k}}=\sigma|{\bm{k}}||\sigma\rangle_{\bm{k}}$. Note that
the operator ${\bm{k}}\cdot\hat{\bf J}$ is invariant under time-reversal, and
so the “helicity” $\sigma$ is _not_ a pseudospin index. Expressed in terms of
the $z$-spin eigenstates we have

$\displaystyle|\tfrac{3}{2}\rangle_{\bm{k}}$ $\displaystyle=$ $\displaystyle
e^{-3i\phi}\cos^{3}(\tfrac{\theta}{2})|\tfrac{3}{2}\rangle+\sqrt{3}e^{-2i\phi}\cos^{2}(\tfrac{\theta}{2})\sin(\tfrac{\theta}{2})|\tfrac{1}{2}\rangle$
$\displaystyle+\sqrt{3}e^{-i\phi}\cos(\tfrac{\theta}{2})\sin^{2}(\tfrac{\theta}{2})|-\tfrac{1}{2}\rangle+\sin^{3}(\tfrac{\theta}{2})|-\tfrac{3}{2}\rangle$
$\displaystyle|\tfrac{1}{2}\rangle_{\bm{k}}$ $\displaystyle=$
$\displaystyle-\sqrt{3}e^{-3i\phi}\sin(\tfrac{\theta}{2})\cos^{2}(\tfrac{\theta}{2})|\tfrac{3}{2}\rangle$
(31)
$\displaystyle+\frac{1}{2}e^{-2i\phi}[-1+3\cos(\theta)]\cos(\tfrac{\theta}{2})|\tfrac{1}{2}\rangle$
$\displaystyle+\frac{1}{2}e^{-i\phi}[1+3\cos(\theta)]\sin(\tfrac{\theta}{2})|-\tfrac{1}{2}\rangle$
$\displaystyle+\sqrt{3}\sin^{2}(\tfrac{\theta}{2})\cos(\tfrac{\theta}{2})|-\tfrac{3}{2}\rangle$
$\displaystyle|-\tfrac{1}{2}\rangle_{\bm{k}}$ $\displaystyle=$
$\displaystyle\sqrt{3}e^{-3i\phi}\sin^{2}(\tfrac{\theta}{2})\cos(\tfrac{\theta}{2})|\tfrac{3}{2}\rangle$
(32)
$\displaystyle-\frac{1}{2}e^{-2i\phi}[1+3\cos(\theta)]\sin(\tfrac{\theta}{2})|\tfrac{1}{2}\rangle$
$\displaystyle+\frac{1}{2}e^{-i\phi}[-1+3\cos(\theta)]\cos(\tfrac{\theta}{2})|-\tfrac{1}{2}\rangle$
$\displaystyle+\sqrt{3}\sin(\tfrac{\theta}{2})\cos^{2}(\tfrac{\theta}{2})|-\tfrac{3}{2}\rangle$
$\displaystyle|-\tfrac{3}{2}\rangle_{\bm{k}}$ $\displaystyle=$
$\displaystyle-e^{-3i\phi}\sin^{3}(\tfrac{\theta}{2})|\tfrac{3}{2}\rangle+\sqrt{3}e^{-2i\phi}\sin^{2}(\tfrac{\theta}{2})\cos(\tfrac{\theta}{2})|\tfrac{1}{2}\rangle$
$\displaystyle-\sqrt{3}e^{-i\phi}\sin(\tfrac{\theta}{2})\cos^{2}(\tfrac{\theta}{2})|-\tfrac{1}{2}\rangle+\cos^{3}(\tfrac{\theta}{2})|-\tfrac{3}{2}\rangle$

where $(\phi,\theta)$ are the angular coordinates defined by the direction of
${\bm{k}}$. The helicity index reverses under inversion ${\cal I}$ but is
preserved under time-reversal ${\cal T}$:

$\displaystyle{\cal
I}|\pm\tfrac{3}{2}\rangle_{\bm{k}}=|\mp\tfrac{3}{2}\rangle_{\bf-k}\,,\quad{\cal
I}|\pm\tfrac{1}{2}\rangle_{\bm{k}}=|\mp\tfrac{1}{2}\rangle_{\bf-k}\,,$ (34)
$\displaystyle{\cal T}|\pm\tfrac{3}{2}\rangle_{\bm{k}}=\mp
e^{3i\phi}|\pm\tfrac{3}{2}\rangle_{\bf-k}\,,\quad{\cal
T}|\pm\tfrac{1}{2}\rangle_{\bm{k}}=\pm
e^{3i\phi}|\pm\tfrac{1}{2}\rangle_{\bf-k}\,,$ (35)

Because of the momentum-dependent quantization axis, it is useful to define
the momentum-dependent angular-momentum operators

$J^{{\bm{k}}}_{z}=\hat{\bf k}\cdot{\bf J}\,,\quad
J_{\pm}^{\bm{k}}=\left(\hat{\bm{\theta}}_{\bm{k}}\pm
i\hat{\bm{\phi}}_{\bm{k}}\right)\cdot{\bf J}$ (36)

where $\hat{\bm{\theta}}_{\bm{k}}=\cos\phi\cos\theta\hat{\bf
x}+\sin\phi\cos\theta\hat{\bf y}-\sin\theta\hat{\bf z}$ and
$\hat{\bm{\phi}}_{\bm{k}}=-\sin\phi\hat{\bf x}+\cos\phi\hat{\bf y}$ are the
canonical spherical-coordinate unit vectors in the directions perpendicular to
$\hat{\bm{k}}$.

###  B.2 The Zeeman field

We include a Zeeman field in Hamiltonian as

$H_{\text{Z}}=g\mu_{B}{\bf h}\cdot{\bf J}$ (37)

To understand the effect of the spin-orbit coupling we now re-express this in
terms of the momentum-dependent angular momentum operators:

$\displaystyle H_{\text{Z}}^{\bm{k}}$ $\displaystyle=$ $\displaystyle
g\mu_{B}h\left[\hat{\bf h}\cdot\hat{\bm{k}}J_{z}^{\bm{k}}+\hat{\bf
h}\cdot\left(\hat{\bm{\theta}}_{\bm{k}}+i\hat{\bm{\phi}}_{\bm{k}}\right)J_{+}^{\bm{k}}\right.$
(38) $\displaystyle\left.+\hat{\bf
h}\cdot\left(\hat{\bm{\theta}}_{\bm{k}}-i\hat{\bm{\phi}}_{\bm{k}}\right)J_{-}^{\bm{k}}\right]\,.$

Note that for ${\bm{k}}$ in the direction of the applied field, only the first
term is present, and the Zeeman field is diagonal in the helical basis.
Conversely, for ${\bm{k}}$ perpendicular to the applied field, only the second
and third terms are present, and the Zeeman field couples the helical states
with helicity $\sigma$ differing by $\pm 1$.

###  B.3 The pairing states

####  B.3.1 The $p$-wave $A_{1u}$ state

The $p$-wave $A_{1u}$ has

$\tilde{\Delta}({\bm{k}})={\bm{k}}\cdot{\bf J}=|{\bm{k}}|J_{z}^{\bm{k}}$ (39)

Since $\tilde{\Delta}({\bm{k}})$ is proportional to $J_{z}^{\bm{k}}$, and the
helicity is preserved under time-reversal, this pairs states with the same
helicity. This implies that in the spherical limit the $A_{1u}$ state
corresponds to purely intraband pairing.

####  B.3.2 The $p$-wave $A_{2u}$ state

The $p$-wave $A_{2u}$ has the pair potential

$\displaystyle\tilde{\Delta}({\bm{k}})$ $\displaystyle=$
$\displaystyle\sum_{i=x,y,z}k_{i}[J_{i+1}J_{i}J_{i+1}-J_{i+2}J_{i}J_{i+2}]$
(40) $\displaystyle=$ $\displaystyle|{\bm{k}}|\left(\frac{1}{32}[6\cos
2\phi\sin 2\theta+3i(3+\cos 2\theta)\sin\theta\sin
2\phi]J^{\bm{k}}_{+}J^{\bm{k}}_{+}J^{\bm{k}}_{+}\right.$
$\displaystyle+\frac{1}{32}[6\cos 2\phi\sin 2\theta-3i(3+\cos
2\theta)\sin\theta\sin 2\phi]J^{\bm{k}}_{-}J^{\bm{k}}_{-}J^{\bm{k}}_{-}$
$\displaystyle+\frac{1}{8}[4\cos 2\phi\cos 2\theta+i(1+3\cos
2\theta)\cos\theta\sin 2\phi]J^{\bm{k}}_{+}J^{\bm{k}}_{z}J^{\bm{k}}_{+}$
$\displaystyle+\frac{1}{8}[4\cos 2\phi\cos 2\theta-i(1+3\cos
2\theta)\cos\theta\sin 2\phi]J^{\bm{k}}_{-}J^{\bm{k}}_{z}J^{\bm{k}}_{-}$
$\displaystyle-\frac{1}{128}\left[2\cos 2\phi\sin 2\theta+i(1+3\cos
2\theta)\sin\theta\sin
2\phi\right](7J^{\bm{k}}_{+}-20J^{\bm{k}}_{z}J^{\bm{k}}_{+}J^{\bm{k}}_{z})$
$\displaystyle\left.-\frac{1}{128}\left[2\cos 2\phi\sin 2\theta-i(1+3\cos
2\theta)\sin\theta\sin
2\phi\right](7J^{\bm{k}}_{-}-20J^{\bm{k}}_{z}J^{\bm{k}}_{-}J^{\bm{k}}_{z})\right)$

Although very complicated compared to the $A_{1u}$ case, we observe that the
pair potential never pairs electrons with the same helicity, since it always
involves operators which change the helicity. That is, projected into the band
basis, the pair potential only pairs the states $|\sigma\rangle_{\bm{k}}$ and
$|-\sigma\rangle_{-{\bf k}}$. Moreover, we have the following special cases:

  * •

along the crystal axes $(\phi,\theta)=(0,\frac{\pi}{2}),$
$(\frac{\pi}{2},\frac{\pi}{2})$, $(\phi,0)$, only the coefficients of
$J_{\pm}^{\bm{k}}J_{z}^{\bm{k}}J_{\pm}^{\bm{k}}$ are nonzero, i.e. the pairing
potential pairs states with helicity differing by $\pm 2$. This cannot be
satisfied in either band, and so the intraband gap develops a node.

  * •

along the crystal body diagonals
$(\phi,\theta)=(\pm\frac{\pi}{4},\pm\arctan\sqrt{2})$ only the coefficients of
$J_{\pm}^{\bm{k}}J_{\pm}^{\bm{k}}J_{\pm}^{\bm{k}}$ are nonzero, i.e. the
pairing potential pairs states with helicity differing by $\pm 3$, i.e. there
is only pairing in the spin-$\frac{3}{2}$ band. Consequently the gap in the
spin-$\frac{1}{2}$ band has nodes.


Phase Recognition in Contrast-Enhanced CT Scans based on Deep Learning and
Random Sampling
Binh T. Dao1,†, Thang V. Nguyen1,†, Hieu H. Pham1,2,3, and Ha Q. Nguyen1,2
1 Smart Health Center, VinBigData JSC, Hanoi, Vietnam
2 College of Engineering & Computer Science, VinUniversity, Hanoi, Vietnam
3 VinUni-Illinois Smart Health Center, Hanoi, Vietnam

${\dagger}$ These authors contributed equally to this work
Corresponding author<EMAIL_ADDRESS>(Hieu H. Pham)

###### Abstract

Purpose: A fully automated system for interpreting abdominal computed
tomography (CT) scans with multiple phases of contrast enhancement requires an
accurate classification of the phases. Current approaches to classify the CT
phases are commonly based on 3D convolutional neural network (CNN) approaches
with high computational complexity and high latency. This work aims at
developing and validating a precise, fast multi-phase classifier to recognize
three main types of contrast phases in abdominal CT scans.

Methods: We propose in this study a novel method that uses a random sampling
mechanism on top of deep CNNs for the phase recognition of abdominal CT scans
of four different phases: non-contrast, arterial, venous, and others. The CNNs
work as a slice-wise phase prediction, while the random sampling selects input
slices for the CNN models. Afterward, majority voting synthesizes the slice-
wise results of the CNNs, to provide the final prediction at scan level.

Results: Our classifier was trained on 271,426 slices from 830 phase-annotated
CT scans, and when combined with majority voting on 30% of slices randomly
chosen from each scan, achieved a mean F1-score of 92.09% on our internal test
set of 358 scans. The proposed method was also evaluated on 2 external test
sets: CTPAC-CCRCC (N = 242) and LiTS (N = 131), which were annotated by our
experts. Although a drop in performance has been observed, the model
performance remained at a high level of accuracy with a mean F1-score of
76.79% and 86.94% on CTPAC-CCRCC and LiTS datasets, respectively. Our
experimental results also showed that the proposed method significantly
outperformed the state-of-the-art 3D approaches while requiring less
computation time for inference.

Conclusions: In comparison to state-of-the-art classification methods, the
proposed approach shows better accuracy with significantly reduced latency.
Our study demonstrates the potential of a precise, fast multi-phase classifier
based on a 2D deep learning approach combined with a random sampling method
for contrast phase recognition, providing a valuable tool for extracting
multi-phase abdomen studies from low veracity, real-world data.

###### Contents

  1. I. Introduction
  2. II. Proposed Approach
    1. II.A. Overview of approach
    2. II.B. Data collection and annotation
    3. II.C. Model development
      1. II.C.1. Network architecture
      2. II.C.2. Training methodology
  3. III. Experiments and Results
    1. III.A. Experimental setup
    2. III.B. Evaluation metrics
    3. III.C. Experimental results
      1. III.C.1. Model performance on internal test set
      2. III.C.2. Comparison to state-of-the-art
      3. III.C.3. Model performance on external test set
  4. IV. Discussion
    1. IV.A. Key findings
    2. IV.B. Future works
  5. V. Conclusion
  6. VI. Acknowledgements
  7. VII. Conflict of Interest
  8. References
  9. VIII. Supplementary Materials
    1. VIII.A. EfficientNet architecture

##  I. Introduction

Contrast enhancement in computed tomography (CT) scans, especially of the
abdomen, is crucial for successful lesion diagnosis 1, 2. Certain types of
lesions can only be observed on the CT scans taken after the injection of
contrast agents into the blood veins. The contrast enhancement process
generally consists of three main phases 3 as follows.

  * •

Non-Contrast: the CT scan is acquired without injection of any contrast
agents;

  * •

Arterial: the CT scan is acquired 35-40 seconds after the bolus injection,
which can help with identifying the hepatocellular carcinoma (HCC), the focal
nodular hyperplasia (FNH), and the adenoma in the liver;

  * •

Venous: the CT scan is acquired 70-80 seconds after the bolus injection, in
which the liver parenchyma is enhanced through the blood supply by the portal
vein, highlighting hypovascular liver lesions.

Machine learning algorithms over the last decades have achieved a great
success in the interpretation and diagnosis of medical imaging data 4,
including the automatic detection of liver lesions in contrast-enhanced CT
images 5, 6, 7, 8. For instance, a multi-phase analysis of abdominal CT scans
9 was performed to detect cirrhosis and HCC liver. To obtain a robust
performance, however, such an algorithm often requires training from a large-
scale dataset of patient’s preoperative multi-phase CT scans and clinical
features 8. Hence, a reliable method is needed for the collection and
annotation of imaging data of abdominal structures 10. The data mining process
usually starts with accessing and collecting retrospective medical imaging
data through the Picture Archiving and Communication Systems (PACS).
Unfortunately, the current generation of PACS systems does not support the
curation of large-scale, multi-phase contrast-enhanced CT datasets. The key
obstacle is that DICOM tags related to series description (e.g., non-contrast,
arterial, or venous) are manually input, non-standardized, and often
incomplete 11. These limitation leads to the impossibility to automatically
categorize medical image data based solely on their DICOM metatags, as around
15% of all studies were mislabeled due to human factor 12. As a result, these
datasets often rely on physicians for manual re-annotation of CT scans,
typically expensive and time-consuming.

In addition, a typical machine learning-based computer-aided diagnosis (CAD)
system for interpreting abdominal CT scans is often trained images from a
specific phase, or stacks of images from different phases in some fixed order
13. Consequently, in the deployment scenario, it is essential that the system
knows precisely which phase each CT scan (series) belongs to, so that the
right phase scans can be fed to models. This leads to a dire need for a phase
identification module for abdominal CT series.

Several approaches 14, 15 have been proposed to identify multiple phases from
CT scans. For instance, Zhou et al. 14 focuses on volumetric characteristics
of the scan, resizes the entire scan to a 3D block of size 32$\times$
128$\times$128 by interpolation, then feeds them through a 3D CNN. In the
original paper, the proposed model is trained with 43,000 scans, which is a
significant amount of data to be acquired, and this sparks the question of
whether this method would perform well on our smaller dataset. Another work by
Tang et al. 15 suggests using a generative adversarial network on each slice
instead of the whole 3D scan. The authors state that CD-GAN is trained for a
period of approximately 36 hours on an NVIDIA 2080Ti GPU with 11G memory,
showing that this method is too computationally expensive.

Different from the previous approaches 14, 15, we aim to develop a fast,
highly accurate deep learning system for recognizing phases from CT scans.
Specifically, we propose an efficient strategy solely based on a 2D
representation of the slices. The proposed system consists of two main stages:
(1) Random sampling that randomly picks R% of slices from the input CT scan
and inputs to the deep learning model. Here, R% denotes the percentage of
slices selected from the input scan; and (2) Slice-level prediction that
identifies phases of each chosen slice then uses majority voting to conclude
the phase of the given scan. Our experimental results on internal and external
(i.e., CTPAC-CCRCC 16, LiTS 17) datasets showed that the proposed method
significantly outperforms the state-of-the-art 3D approaches while requiring
less computation time for inference.

To summarize, the main contributions of this work are the following:

$\bullet$ We develop and evaluate a novel deep learning system for the
recognition of multi-phase in contrast-enhanced CT scans. The proposed system
exploits a random sampler to reduce the reduce computational time of the input
examples. Majority voting is used to boost the final prediction of the system.
Our extensive experiments show that the proposed approach surpasses previous
state-of-the-art 3D approaches in terms of both accuracy and inference time.
The proposed deep learning system can be easily reused or finetuned, therefore
potential benefits for several applications in clinical settings.

$\bullet$ The imaging dataset used in this study will be shared on our project
website at https://vindr.ai/datasets/abdomen-phases, while the codes will be
published at https://github.com/vinbigdata-medical/abdomen-phases. To the best
of our knowledge, this is the biggest annotated dataset for the recognition of
multi-phase in contrast-enhanced CT scans.

##  II. Proposed Approach

###  II.A. Overview of approach

Our main goal in this study is to develop and evaluate a fast, accurate deep
learning system for the recognition of multi-phase in contrast-enhanced CT
scans. To this end, we first randomly sample $R$% of the slices from the whole
original scan. Each of these chosen slices is then passed through a CNN model,
which was trained to output the phase classification at the slice level.
Finally, the scan-level prediction is predicted by a majority voting of the
results obtained in the previous step. The proposed scheme for the phase
recognition of abdominal CT scans is illustrated in Figure 2.

###  II.B. Data collection and annotation

To develop deep learning algorithms for contrast phase recognition, we built
an internal dataset of abdominal CT scans. The construction of this dataset
was divided into three main steps: (1) data collection, (2) data de-
identification, and (3) data annotation.

Data collection: A total 265 abdominal studies comprising 1,188 CT scans in
the Digital Imaging and Communications in Medicine (DICOM) format were
retrospectively randomly selected from the PACS databases of the Hospital 108
and the Hanoi Medical University Hospital – two major hospitals in Vietnam –
within the period from 2015 until 2020. The ethical clearance of this study
was approved by the Institutional Review Board (IRB) of each hospital before
any data processing steps. The need for obtaining patient consent was waived
because these studies did not impact clinical care.

Data characteristics, including patient demographics and the prevalence of
each contrast-phase class, are summarized in Table 1. A general statistics of
the slice and scan distribution for each class is featured in Figure 3. The
distribution of the CT scanner models and their manufacturers are shown in
Figure 4. There are 6 different levels of slice thickness in the entire
dataset, whose distribution is captured in Table 2. Additionally, the number
of slices per scan ranges from 30 to 2,350 with a mean of 281. The
distribution of the number of slices per scan is illustrated in Figure 5.

Data de-identification: To protect patient’s privacy, all personally
identifiable information associated with the DICOM images has been removed.
Specifically, a Python script was written to remove all DICOM tags of
protected health information (PHI) such as patient’s name, patient’s date of
birth, patient ID, etc. We only retained a limited number of DICOM attributes
that are necessary for processing raw images. The full list of these DICOM
attributes is provided in Supplement LABEL:list-dicom-tag, Table 9.

Data annotation: The dataset was labeled for a total of 4 contrast-phase
classes: (1) Non-contrast, (2) Venous, (3) Arterial, and (4) Others. Here the
Others category refers to all scans that cannot be correctly classified as
either of 3 phases Non-contrast, Venous, and Arterial. They may include scans
of the delay phase or scans that belong to a transitional state between 2
phases. To annotate the imaging data, we designed and built a web-based
labeling framework called VinDr Lab (https://vindr.ai/vindr-lab) 18. Two
radiologists were hired to remotely annotate the data. Once the labeling has
been completed, the labels were exported in Comma-Separated Values (CSV)
format and will be used for training deep learning algorithms.

In total, 265 studies have been annotated. The whole dataset was then divided
into training and validation set by a ratio of 70%/30% accordingly. Since each
study usually contains multiple scans of the same patient, the train-test
split was stratified by the study level to avoid data leakage. As a result,
our training set consists of 271,426 slices from 830 scans (186 studies),
while our validation set contains 121,134 slices from 358 scans (79 studies).

Data Records: To encourage new advances in this research direction, we will
make the dataset freely accessible via our project website at
https://vindr.ai/datasets/abdomen-phases. Specifically, all imaging data and
the corresponding ground truth labels for the training and validation set will
be provided. The images are organized into two folders, one for training and
the other one for validation in which each image has a unique, anonymous
identifier.

###  II.C. Model development

This section describes in detail our model development method. We exploit
state-of-the-art, high-performing deep CNN architectures for the task of
recognizing multi-phase in contrast-enhanced CT scans. We describe our network
architecture choice and training methodology as the following.

####  II.C.1. Network architecture

A set of state-of-the-art deep CNN models has been deployed and evaluated on
the collected dataset, including ResNet-18 19, ResNet-34 19, SEResNet-18,
ResNext-50 20, EfficientNet-B0 21, EfficientNet-B2 21, GhostNet 22, and CD-GAN
23. These deep networks are well-known to be effective for image recognition
tasks. Each network accepts a CT scan as input and predicts the corresponding
contrast phase label. For implementation, we followed the same instructions
and recommendations from the original papers 19, 20, 21, 22, 23. We considered
EfficientNet 21 model as our main network architecture choice due to the high-
level of accuracy and efficiency of this architecture compared to previous
deep CNNs. Details of the EfficientNet 21 architecture is provided in
Supplement VIII.A..

####  II.C.2. Training methodology

In the training stage, all images were fed into the networks with a size of
224$\times$224 pixels. Input images extracted from raw DICOM files were,
first, converted to standard Hounsfield units (HU), using Rescale Slope and
Rescale Intercept from DICOM headers. Afterward, we applied HU window with the
window center of 50 and the window width of 400 to the image. During the
training process, we used Adam optimizer 24 with an initial learning rate of
$10^{-2}$ and cosine annealing scheduler 25 with linear warm up 26. Each
network was trained end-to-end for 15 epochs. To this end, we minimized the
binary cross-entropy loss function between the ground-truth labels and the
predicted label by the network over the training samples. The proposed deep
network was implemented in Python using Pytorch version 1.7.1
(https://pytorch.org/). All experiments were conducted on a Ubuntu 18.04
machine with a single NVIDIA Geforce RTX 2080 Ti with 11GB memory.

##  III. Experiments and Results

###  III.A. Experimental setup

Internal validation: Extensive experiments were conducted to evaluate the
performance of the proposed method. Specifically, we first evaluated the
slice-wise classification performance of trained CNN models (i.e., ResNet-18
19, ResNet-34 19, SEResNet-18, ResNext-50 20, EfficientNet-B0 21,
EfficientNet-B2 21, GhostNet 22, and CD-GAN 23) on the validation set of
121,134 slices. Next, we reported the classification performance of the best
performing network at the scan-wise level by applying the majority voting on
R% of the slices selected from each scan. We experimented with R ranging from
1 to 20 at an interval of 5, then 20 to 100 at an interval of 10. Finally, to
compare the proposed 2D approach with previous 3D state-of-the-art approaches,
we re-implemented two-phase recognition approaches on CT scans including 3DSE
14 and CD-GAN 23. These approaches were trained on the training dataset using
the same hyperparameters setting as described in the original papers 14, 23.
We also measured the average inference time (second) per scan for each
approach and compared it with our proposed 2D method.

External validation: To verify the generalization ability of the proposed deep
learning model, we evaluated it on two external datasets, including LiTS 17
and CPTAC-CCRCC 16. The LiTS 17 dataset contains 131 CT scans in the training
set and 70 CT scans in the test set. It was original developed for the
development of liver segmentation methods. The CPTAC-CCRCC 16 was introduced
by the National Cancer Institute’s Clinical Proteomic Tumor Analysis
Consortium (CPTAC) and developed for investigating the clear cell renal cell
carcinoma (CCRCC). We utilized the imaging data for CCRCC tumors, containing
242 CT scans for our external testing. However, the original aim of these
datasets did not match the purpose of this study. As a result, there were no
target labels for this dataset. Our radiologist team, therefore, classified
scans from these datasets into 4 phase categories. As a result, the LiTS 17
dataset has 8 scans of the Arterial phase, and 123 scans of the Venous phase.
Meanwhile, the CPTAC-CCRCC 16 contains 57, 69, 53, and 63 scans from 4
categories Non-Contrast, Venous, Arterial, and Others, respectively.

###  III.B. Evaluation metrics

We report the classification performance using mean accuracy, macro-average
precision precision, macro-average precision recall, macro-average precision
and macro-average F1-score. These performance indicators are defined as
follows.

Accuracy
$\displaystyle=\frac{\text{TP}+\text{TN}}{\text{TP}+\text{FP}+\text{TN}+\text{FN}},$
Precision $\displaystyle=\frac{\text{TP}}{\text{TP}+\text{FP}},$ Recall
$\displaystyle=\frac{\text{TP}}{\text{TP}+\text{FN}},$ F1
$\displaystyle=2\times\frac{\text{Precision}\times\text{Recall}}{\text{Precision}+\text{Recall}}.$

Here, TP, FP and FN are the number of True Positive, False Positive and False
Negative samples accordingly.

###  III.C. Experimental results

####  III.C.1. Model performance on internal test set

Table 3 summarizes quantitative results for several state-of-the-art CNN
classification models on the internal test set of 121,134 slices. Note that,
while training and benchmarking those models, we chose to fix the input image
size to be 128$\times$128 for the sake of saving computations. It can be seen
that EfficientNet-B2 achieved the best performance with a macro-averaged
recall of 85.92%, a macro-averaged precision of 84.70%, and a macro-averaged
F1-score of 85.26%. This architecture was then selected to conduct all the
remaining experiments.

We further investigated the impact of different image input sizes on the
performance of the picked model, i.e. EfficientNet-B2, as shown in Table 4. We
observed that using the input images with a size of 224$\times$224 for
training gave us the best result: a macro-averaged accuracy of 93.51%; a
macro-averaged recall of 85.46%; a macro-averaged precision of 87.45% and a
macro-averaged F1-score of 86.43%. In addition, training the model with the
224$\times$224 images only took 10 minutes for each epoch instead of 60
minutes when using the input images of size 512$\times$512.

The classification performance of EfficientNet-B2 on our test set is shown in
Table 5. We computed the precision, recall, and F1-score for both the slice
and scan levels along with their 95% confidence interval (CI) using
bootstrapping over 5,000 resamples of the test set. Our proposed model
achieved a mean F1-score of 0.8643 (95% CI 0.8602, 0.8664) for the slice-level
prediction and a mean F1-score of 0.9209 (95% CI 0.9033, 0.9374) for the scan-
level prediction with majority voting on 30% of the total slices. We observed
that the reported performance remained consistent between our three main
classes: Non-Contrast, Venous, and Arterial, while there was a visible gap
between Others and the rest. Additionally, the ROC (Receiver Operating
Characteristic) curves of the proposed model for the four classes are plotted
in Figure 6 along with their corresponding AUC (Area Under the ROC Curve)
scores on the slice-level test set. Unlike F1-score, AUC is a threshold-
independent metric. Nevertheless, it can be seen that the AUC scores reported
in Figure 6 are strongly correlated with the slice-wise F1-scores given in
Table 5.

The effect of using different values of R when performing random sampling with
the majority voting for scan-level prediction is illustrated in Figure 7. It
can be clearly seen that the macro-average F1-score increased as R (the
percentage of slices in each scan to be selected randomly for inference)
approaches 20%, and leveled off as R increased. By applying $R=30$, we
observed a 5.66% increase in macro-averaged F1-score compared to the slice-
wise performance.

####  III.C.2. Comparison to state-of-the-art

To demonstrate the effectiveness of the proposed 2D approach, we compared our
result with recent state-of-the-art methods 14, 23 for the recognition of
multi-phase in contrast-enhanced CT scans. To this end, we reproduced the 3DSE
by Zhou et al. 14 and the CD-GAN 23 by Tang et al. and reported their
performance of these approaches using F1-score on the test set. For a fair
comparison, we applied the same training methodologies and hyper-parameter
settings as reported in the original papers 14, 23. In particular, the input
image size to the model was fixed to 128$\times$128 while comparing to CD-GAN.
The experimental results are provided in Table 6. We found that the proposed
2D approach significantly surpassed the previous state-of-the-art approaches
(an improvement of 6.09% compared to the 3DSE 14 and 3.07% compared to CD-GAN
23), while required less time for inference.

####  III.C.3. Model performance on external test set

We reported in Table 7 the experimental results on two external test sets LiTS
and CPTAC-CCRCC. The averaged F1-scored on the LiTS was 86.94%, while the
averaged F1-score on the CPTAC-CCRCC was 76.79%. We found that the proposed
method suffered from covariate shift, however, it still remains at a high
level of F1-score.

##  IV. Discussion

###  IV.A. Key findings

The phase recognition is important for medical imaging data collection and the
deployment of machine learning models in practice. From the clinical
perspective, a method for fast and precise recognition of CT phases can
effectively aid in diagnosis of abdominal pathologies 27. By training a set of
strong deep CNN models on a large-scale, annotated dataset, we built an
automated system that is able to accurately recognize contrast multi-phases
from CT scans. In particular, we empirically showed that a major improvement
has been achieved, in terms of F1-score and inference time by applying the
proposed random sampling and majority voting. Compared to previous state-of-
the-art 3D approaches, our model showed 30 times improved in inference time
and nearly 6% improvement on F1-score on our dataset.

Although a highly accurate performance has been achieved across three classes:
Non-Contrast, Arterial and Venous, we acknowledge that the proposed method
reveals some limitations. To make a correct classification of contrast phases,
experts often rely on the multiple slices containing arteries, veins, and
parenchyma. However, in our method, slices from each scan are predicted
independently, without incorporating information from other regions. In
addition, since contrast materials are absorbed differently for each
individual, slices of the same regions from 2 different phases, such as the
Arterial and Venous phases, could have similar brightness in the arteries. An
example is demonstrated in Figure 8: there is a clear brightness difference
between the 2 images in the Arterial phase and the top left Arterial image
resembles slices from the Venous phase. For this reason, our slice-level
predictions are prone to errors. Another challenge has related to the nature
of our dataset. Our samples vary in scan range, some abdominal CT scans can
include neck or thighs where contrast material does not pass through, making
these slices indifferentiable across our 4 classes. Therefore, our model is
likely to produce a false prediction on these slices, which contributes
negatively to the performance of the scan level prediction.

The low generalizability of deep learning-based diagnostic systems 28, 29, 30,
31 to datasets and scanners beyond the ones they have been trained with, has
been limiting the use of such methods in real-world clinical settings. We
showed that the proposed deep learning method was successfully generalized to
two different datasets from other hospital sites, for which each with a
different CT scanner.

###  IV.B. Future works

There are several possible mechanisms to improve our current method. The most
promising direction is to eliminate non-affected contrast-enhanced regions of
a scan, such as pelvis areas. This would improve slice-level prediction since
the model is forced to learn and predict images with clearer features. In
addition, due to the performance drop in “Others”, future work includes
applying techniques for reducing the impact of imbalanced data. For example,
weighted BCE losses 32, which directly penalized the probabilistic false
positives, can be used. We also plan to experiment with training and testing
the proposed method on the coronal and sagittal projections of the CT scans,
so that each input image could contain all necessary components: arteries,
veins, and parenchyma, which are used to identify the correct contrast phases.
Moreover, we would conduct additional experiments, incorporating training
procedure refinements 33 such as data augmentation methods to further improve
the generalization of our method. Lastly, it is worth investigating more
sophisticated methods for sampling and synthesizing the slice-level
predictions, such as the multi-instance learning paradigm 34, rather than the
straightforward random sampling and majority voting strategies discussed in
the paper.

##  V. Conclusion

In this study, we developed a 2D deep learning-based approach for the
recognition of contrast phases in abdominal CT scans. We adopted a random
sampling strategy to improve the classification performance and to reduce
inference time. The introduction of a random sampling mechanism helps avoid
training and inferring on 3D data, which are usually much more costly, while
still attaining impressive performances. Extensive experimental results on
both the internal and external datasets have demonstrated that the proposed
approach significantly outperformed previous state-of-the-art 3D approaches.

##  VI. Acknowledgements

This work was supported by VinBigData JSC. We are especially thankful to all
of our collaborators, including radiologists, physicians, and technicians, who
participated in the data collection and labeling process.

##  VII. Conflict of Interest

The authors have no conflict of interest to disclose.

## References

##

  * 1 G. Brancatelli, M. P. Federle, L. Grazioli, A. Blachar, M. S. Peterson, and L. Thaete,  Focal nodular hyperplasia: CT findings with emphasis on multiphasic helical CT in 78 patients,  Radiology 219, 61–68 (2001).
  * 2 Y. L. Bronstein, E. M. Loyer, H. Kaur, H. Choi, C. David, R. A. DuBrow, L. D. Broemeling, K. R. Cleary, and C. Charnsangavej,  Detection of small pancreatic tumors with multiphasic helical CT,  American Journal of Roentgenology 182, 619–623 (2004).
  * 3 R. Smithuis,  CT contrast injection and protocols,  https://radiologyassistant.nl/more/ct-protocols/ct-contrast-injection-and-protocols, 2014\.
  * 4 G. Litjens, T. Kooi, B. E. Bejnordi, A. A. A. Setio, F. Ciompi, M. Ghafoorian, J. A. Van Der Laak, B. Van Ginneken, and C. I. Sánchez,  A survey on deep learning in medical image analysis,  Medical Image Analysis 42, 60–88 (2017).
  * 5 K. Yasaka, H. Akai, O. Abe, and S. Kiryu,  Deep learning with convolutional neural network for differentiation of liver masses at dynamic contrast-enhanced CT: a preliminary study,  Radiology 286, 887–896 (2018).
  * 6 W. Wang, Y. Iwamoto, X. Han, Y.-W. Chen, Q. Chen, D. Liang, L. Lin, H. Hu, and Q. Zhang,  Classification of focal liver lesions using deep learning with fine-tuning,  in Proceedings of the International Conference on Digital Medicine and Image Processing, pages 56–60, 2018.
  * 7 Y. Yoshinobu, Y. Iwamoto, H. Xianhua, L. Lanfen, H. Hongjie, Q. Zhang, and C. Yen-Wei,  Deep learning method for content-based retrieval of focal liver lesions using multiphase contrast-enhanced computer tomography images,  in IEEE International Conference on Consumer Electronics (ICCE), pages 1–4, 2020.
  * 8 R. Gao et al.,  Deep learning for differential diagnosis of malignant hepatic tumors based on multi-phase contrast-enhanced CT and clinical data,  Journal of hematology & oncology 14, 1–7 (2021).
  * 9 A. Nayak, E. B. Kayal, M. Arya, J. Culli, S. Krishan, S. Agarwal, and A. Mehndiratta,  Computer-aided diagnosis of cirrhosis and hepatocellular carcinoma using multi-phase abdomen CT,  International journal of computer assisted radiology and surgery 14, 1341–1352 (2019).
  * 10 S. Park et al.,  Annotated normal CT data of the abdomen for deep learning: Challenges and strategies for implementation,  Diagnostic and interventional imaging 101, 35–44 (2020).
  * 11 H. Harvey and B. Glocker,  A standardised approach for preparing imaging data for machine learning tasks in radiology,  in Artificial Intelligence in Medical Imaging, pages 61–72, 2019\.
  * 12 M. O. Gueld, M. Kohnen, D. Keysers, H. Schubert, B. B. Wein, J. Bredno, and T. M. Lehmann,  Quality of DICOM header information for image categorization,  in Medical Imaging 2002: PACS and integrated medical information systems: design and evaluation, volume 4685, pages 280–287, 2002\.
  * 13 C. Sun, S. Guo, H. Zhang, J. Li, M. Chen, S. Ma, L. Jin, X. Liu, X. Li, and X. Qian,  Automatic segmentation of liver tumors from multiphase contrast-enhanced CT images based on FCNs,  Artificial Intelligence in Medicine 83, 58–66 (2017).
  * 14 B. Zhou, A. P. Harrison, J. Yao, C.-T. Cheng, J. Xiao, C.-H. Liao, and L. Lu,  CT data curation for liver patients: phase recognition in dynamic contrast-enhanced CT,  arXiv:1911.06395 [eess.IV], 2019.
  * 15 Y. Tang et al.,  Contrast phase classification with a generative adversarial network,  in Medical Imaging 2020: Image Processing, volume 11313, page 1131310, 2020.
  * 16 S. Kalayci, F. Petralia, P. Wang, and Z. H. Gümüş,  ProNetView-ccRCC: A web-based portal to interactively explore clear cell renal cell carcinoma proteogenomics networks,  Proteomics 20, 2000043 (2020).
  * 17 P. Bilic et al.,  The liver tumor segmentation benchmark (LITS),  arXiv preprint arXiv:1901.04056 (2019).
  * 18 N. T. Nguyen, P. T. Truong, V. T. Ho, T. V. Nguyen, H. T. Pham, M. T. Nguyen, L. T. Dam, and H. Q. Nguyen,  VinDr Lab: A Data Platform for Medical AI,  https://github.com/vinbigdata-medical/vindr-lab, 2021.
  * 19 K. He, X. Zhang, S. Ren, and J. Sun,  Deep residual learning for image recognition,  in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 770–778, 2016.
  * 20 S. Xie, R. Girshick, P. Dollár, Z. Tu, and K. He,  Aggregated residual transformations for deep neural networks,  in Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pages 1492–1500, 2017.
  * 21 M. Tan and Q. V. Le,  EfficientNet: Rethinking model scaling for convolutional neural networks,  in The 36th International Conference on Machine Learning, pages 6105–6114, 2019.
  * 22 K. Han, Y. Wang, Q. Tian, J. Guo, C. Xu, and C. Xu,  Ghostnet: More features from cheap operations,  in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 1580–1589, 2020.
  * 23 Y. Tang, H. H. Lee, Y. Xu, O. Tang, Y. Chen, D. Gao, S. Han, R. Gao, C. Bermudez, M. R. Savona, R. G. Abramson, Y. Huo, and B. A. Landman,  Contrast Phase Classification with a Generative Adversarial Network,  arXiv:1911.06395 [eess.IV], 2019.
  * 24 D. P. Kingma and J. Ba,  Adam: A method for stochastic optimization,  arXiv preprint arXiv:1412.6980 (2014).
  * 25 I. Loshchilov and F. Hutter,  SGDR: Stochastic gradient descent with warm restarts,  arXiv preprint arXiv:1608.03983 (2016).
  * 26 J. Ma and D. Yarats,  On the adequacy of untuned warmup for adaptive optimization,  arXiv preprint arXiv:1910.04209 7 (2019).
  * 27 K. Guite, L. Hinshaw, and F. Lee,  Computed tomography in abdominal imaging: how to gain maximum diagnostic information at the lowest radiation dose,  in Selected Topics on Computed Tomography, IntechOpen, 2013.
  * 28 R. Therrien and S. Doyle,  Role of training data variability on classifier performance and generalizability,  in Medical Imaging 2018: Digital Pathology, volume 10581, page 1058109, 2018.
  * 29 X. Liang, D. Nguyen, and S. B. Jiang,  Generalizability issues with deep learning models in medicine and their potential solutions: illustrated with cone-beam computed tomography (CBCT) to computed tomography (CT) image conversion,  Machine Learning: Science and Technology 2, 015007 (2020).
  * 30 M. J. Willemink, W. A. Koszek, C. Hardell, J. Wu, D. Fleischmann, H. Harvey, L. R. Folio, R. M. Summers, D. L. Rubin, and M. P. Lungren,  Preparing medical imaging data for machine learning,  Radiology 295, 4–15 (2020).
  * 31 S. A. Nadeem, A. P. Comellas, E. A. Hoffman, and P. K. Saha,  Generalizability of a deep learning airway segmentation algorithm to a blinded low-dose CT dataset,  in Medical Imaging 2021: Image Processing, volume 11596, page 115963I, 2021.
  * 32 Y. Ho and S. Wookey,  The real-world-weight cross-entropy loss function: Modeling the costs of mislabeling,  IEEE Access 8, 4806–4813 (2019).
  * 33 T. He, Z. Zhang, H. Zhang, Z. Zhang, J. Xie, and M. Li,  Bag of tricks for image classification with convolutional neural networks,  in Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 558–567, 2019.
  * 34 Z. Yan, Y. Zhan, Z. Peng, S. Liao, Y. Shinagawa, S. Zhang, D. N. Metaxas, and X. S. Zhou,  Multi-Instance Deep Learning: Discover Discriminative Local Anatomies for Bodypart Recognition,  IEEE Transactions on Medical Imaging 35, 1332–1343 (2016).
  * 35 J. Deng, W. Dong, R. Socher, L.-J. Li, K. Li, and L. Fei-Fei,  Imagenet: A large-scale hierarchical image database,  in IEEE Conference on Computer Vision and Pattern Recognition, pages 248–255, 2009.

##  VIII. Supplementary Materials

###  VIII.A. EfficientNet architecture

This section provides a detailed description of EfficientNet 21, which is our
main network architecture choice in this work. In 2019, Mingxing Tan et al.
proposed EfficientNet 21 that achieve state-of-the-art performance on top-1
accuracy on ImageNet 35 dataset by the time their model architecture was first
released. In addition to the improvement of accuracy, the author also claimed
to have significantly improved the model efficiency. One of the main
contribution of Mingxing Tan et al. 21 was the compound scaling method
demonstrated in Figure 9. Width simply means how wide the network is, in terms
of the number of channels in a Conv layer. Wider networks can capture more
fine-grained features, keeping the model small. However, the networks will
saturate quickly with larger width as in-depth scaling. Depth means how many
layers are stacked into the network. Scaling by depth dimension enables the
model to capture richer and more complex features. Yet, as the network gets
deeper, the performance will be affected by gradient vanishing. Resolution
means the resolution of the input image to the network. The higher the
resolution, the more features are captured during training. Nevertheless, this
does not scale linearly, and quickly diminish with the larger image size. The
proposed compound scaling is introduced to solve these problems. The method is
further illustrated in the formula below.

depth: $d=\alpha^{\phi}$

width: $w=\beta^{\phi}$

resolution: $r=\gamma^{\phi}$

s.t. $\alpha\cdot\beta^{2}\cdot\gamma^{2}\approx 2$

$\alpha\geq 1,\beta\geq 2,\gamma\geq 1$

where $\alpha$, $\beta$, and $\gamma$ are factors that can be determined using
grid search. Here $\phi$ controls how many more resources are available for
model scaling, whereas $\alpha$, $\beta$, and $\gamma$ specify how to
distribute these extra resources to the network width, depth, and resolution.
In addition to the novel scaling method, the authors introduced a baseline
EfficientNet-B0, by doing a Neural Architecture Search that optimized both
accuracy and floating-point operations per second (FLOPS). The architecture of
the baseline is depicted in Table 8. The main building block of the network is
the mobile inverted bottleneck MBConv added with squeeze-and-excitation
optimization. By fixing $\alpha$, $\beta$, $\gamma$ and experiment with
different values of $\phi$, we get EfficientNet B1-B7.