\section*{Abstract}
    Coalescences of neutron stars and white dwarfs are relatively frequent phenomena, outnumbering other types of compact object mergers (neutron stars and black holes without involving white dwarfs) altogether. Such event potentially can produce not only optical, but also an X-ray burst. Transient source
    CDF-S XT2 \cite{2019Natur.568..198X} can be an example of this type of events as suggested by \cite{2019MNRAS.488..259F}. In this note we estimate the rate of these transients in the field of view of X-ray instruments on-board Spectrum-RG satellite. We demonstrate that during four years of the survey program several thousand
    of events related to neutron star --- white dwarf mergers  might appear in the field of view of eROSITA. Collimation of X-ray emission can reduce this number. Smaller, but comparable number of transients is expected in the case of ART-XC telescope. However, due to relatively short duration --- $\lesssim 10^4$~s, --- mostly such transients might be visible just in one scan of telescopes ($\sim 40$~s), and so only a few photons are expected to be detected which makes definite identification without additional information very problematic.

\section{Introduction}

Coalescences of compact objects --- white dwarfs (WDs), neutron stars (NSs), and black holes (BHs), --- are powerful phenomena accompanied by different transients: gamma-ray bursts, gravitational wave emission, type Ia supernova, etc.
Rate of different types of mergers varies from one event in few tens of years per a Milky Way-like galaxy for WD-WD
(some fraction of such mergers --- $\sim15$\%, --- results in SN Ia which occur once in $\sim 200$~yrs per a Milky Way-like galaxy) \cite{2018MNRAS.476.2584M}
down to once in $\lesssim 10^5$~yrs per galaxy for NS-BH \cite{2019arXiv190704218P}.\footnote{Rates of compact object coalescence are often given in units [Gpc$^{-3}$~yr$^{-1}$]. Here we re-calculated them in units [per galaxy per year]. In \cite{2019arXiv190704218P} the authors give the value $\sim 10$~--~30~Gpc$^{-3}$~yr$^{-1}$. Rates of BH-BH and NS-NS coalescence are typically a little bit higher: up to $\sim 100$ Gpc$^{-3}$~yr$^{-1}$ \cite{2018arXiv180605820M, 2019arXiv190704218P} and up to $\sim300$ Gpc$^{-3}$~yr$^{-1}$ \cite{2017ApJ...846..170T}, respectively. }

Recently, NS-WD coalescences attracted much interest. As it was de\-mon\-strat\-ed by  \cite{2017MNRAS.467.3556B},  due to angular momentum losses via
disk winds (in contrast with the models in which angular momentum is carried away by a jet),
a wast majority of such events (with $M_\mathrm{WD}\gtrsim 0.2\, M_\odot$) proceeds in an unstable way on a very short time scale (from minutes to weeks depending on the WD mass, see Fig. 13 in \cite{2017MNRAS.467.3556B}).
This rapid merger is necessarily accompanied by a huge energy release.
Thus, bright transients might appear at different wavelengths, including optical range.\footnote{Also gravitational wave signals are expected. They can be detected by eLISA and similar instruments sensitive to low-frequency waves.}

In \cite{2019MNRAS.488..259F} the authors studied the outcome of NS-WD mergers in some details from the point of view of electromagnetic counterparts (see also a recent e-print \cite{2019arXiv190810866Z}). In particular, they suggested that due to jet formation a high-energy transient can appear. Moreover, it was proposed that the X-ray transient CDF-S XT2 discovered by \cite{2019Natur.568..198X} could be an example of such event.

In \cite{2018A&A...619A..53T} the rate of NS-WD mergers was calculated for different ages of stellar populations. These results clearly demonstrate that this type of coa\-le\-sce\-nce outnumbers joint rate of NS-NS, NS-BH, and BH-BH mergers.
Total (time intergrated) number of  NS-WD coalescence produced by a stellar population is: $(3-7)\times 10^{-5}$ per a solar mass of created stars. This corresponds to the rate $\sim (1-2)\times 10^{-4}$~yr$^{-1}$ in a Milky Way-like galaxy (it is important to underline, that mostly mergers happen within the first Gyr after binaries are formed). Note, that this value is $\lesssim10\%$ of the SN Ia rate.

In \cite{2017MNRAS.467.3556B} the authors also estimated the visibility of NS-WD merging systems in X-rays under assumption that the luminosity is equal to the Eddington value. Numbers are low, and it is not very probable that even a single event will be detected by eROSITA (extended ROentgen Survey with an Imaging Telescope Array) aboard Spectrum-RG satellite. However, as was shown in \cite{2017MNRAS.467.3556B}, this instrument can detect in relatively near-by galaxies many (from several tens to several hundreds) NS-WD systems with low-mass WDs evolving through the stage of stable mass transfer.
Here we explore a different possibility.
We estimate the rate of eROSITA (and ART-XC) detections assuming that NS-WD coalescences can produce bright transients, similar to CDF XT2 \cite{2019Natur.568..198X}, due to jet (or another origin) X-ray emission.

\section{Estimate of the cosmic rate of NS-WD mergers}

In this section we estimate cosmic rate of NS-WD mergers integrated over different redshifts.

Hereafter, we assume flat Lambda-CDM cosmological model with the Planck 2018 cosmological parameters \cite{2018arXiv180706209P}: $H_0 = 67.4$ km s$^{-1}$ Mpc$^{-1}$, $\Omega_m = 0.315$, $\Omega_{\Lambda} = 0.685$. %

We use parameters of cosmic star formation history from \cite{2014ARA&A..52..415M}:
\begin{equation}
    \psi (z) = 0.015 \frac{(1+z)^{2.7}}{1 + [(1+z)/2.9]^{5.6}} \text{ M}_{\odot} \text{ yr}^{-1} \text{ Mpc}^{-3},
\end{equation}
where $z$ is redshift.
This function is the best fit for observational UV and IR data.

It is important to note, that our estimates below are sensitive only to the total star formation rate (SFR) at a given $z$, and do not depend on specific SFR in different types of galaxies, etc.

The rate of NS-WD mergers is given by the equation:
\begin{equation}
    \frac{d\dot{N}_\mathrm{NSWD} (z)} {dz} = \frac{1}{1+z} \int_z^{z_{max}} \rho(z, z') \psi(z') \frac{dV_c(z')}{dz'} \frac{dt(z')}{dz'} dz'.
\end{equation}
Here $t$ is the age of the universe and we can write:
\begin{equation}
    \frac{dt(z)}{dz} = - \frac{1}{H_0 (1+z) \sqrt{(1+z)^3 \Omega_m + \Omega_{\Lambda}}}.
\end{equation}
The maximum considered redshift $z_{max} = 11.247$ corresponds to the age of the universe 400 Myr, where intensive star formation began,
$dV_\text{c}$ is differential comoving volume, and $\rho(z, z') = \rho(t(z) - t(z'))$ is determined by the  delay time function of NS-WD mergers. For the later quantity we use results obtained with the model $\alpha \alpha$ (with kick velocity distribution described in \cite{2017A&A...608A..57V}) from \cite{2018A&A...619A..53T}.

The $\alpha \alpha$-model is based on the classical energy balance of a common-envelope (CE) evolution:
\begin{equation}
    E_{\text{bin}} = \frac{G M_\text{d} M_\text{c}}{\lambda R} = \alpha_{\text{CE}} E_{\text{orbit}},
\end{equation}
where $M_\text{d}$ is the mass of the donor star, $M_\text{c}$ --- the mass of its core, $R$ --- its radius, $\lambda$ -- the structure parameter of its envelope, and $E_{\text{bin}}$ is the binding energy of the envelope. The orbital energy $E_{\text{orbit}}$ is considered as the source of energy to unbind the envelope with an efficiency determined by the parameter $\alpha_{\text{CE}}$.
The model assumes $\alpha_{\text{CE}} \lambda = 2$. Kick velocity distribution is taken from \cite{2017A&A...608A..57V}.

To calculate the differential comoving volume (and to take into account cosmological effects below) we use equations from \cite{1999astro.ph..5116H}:
\begin{equation}
    \frac{dV_c(z)}{dz} = \frac{c}{H_0} \frac{4 \pi D_M^2(z)}{\sqrt{(1+z)^3 \Omega_m + \Omega_{\Lambda}}}.
\end{equation}
Here $\Omega_m$ is the density of non-relativistic matter in units of the critical density, $\Omega_{\Lambda}$ is the density of the cosmological constant in units of the critical density, $D_\text{M}$ is comoving distance.

The total rate of NS-WD mergers intergrated over all redshifts is $\sim$ 843
per day for the whole sky. The maximum of the differential distribution of NS-WD merger rate is at $z = 1.8$.

Transient X-ray emission can be mildly collimated. We assume that jet emission fills solid angle  $\sim1$ steradian. Thus, the number of potentially observable transients is reduced by the factor $(1/4\pi)$ in comparison with the merger rate. The obtained distribution of transients rate is shown in Fig.~1.

\begin{figure}
\includegraphics[scale=0.9]{trans_rate_new.eps}
\caption{The rate of NS-WD transients per day per unit redshift interval for different redshifts. It is assumed that collimated emission covers on a sphere the solid angle equal to one steradian. Total number of transients per day from all redshifts is equal to 67.}
\end{figure}

\section{Log $N$ --- Log $S$ distribution for NS-WD mergers}

In \cite{2019MNRAS.488..259F} the authors obtain peak luminosities of transients $\lesssim 10^{47}$ erg s$^{-1}$. The peak isotropic luminosity of the CDF-S XT2 event was $L_\mathrm{X} \sim 3 \times 10^{45}$ erg s$^{-1}$.
We construct the Log $N$ --- Log $S$ distribution for three characteristic luminosity values (Fig.~2) assuming that the luminosity distribution is a delta-function and that X-ray emission is collimated within a solid angle one steradian. The flux is given by the formula:
\begin{equation}
    F = \frac{L}{4 \pi D_\mathrm{L}^2},
\end{equation}
where $D_\mathrm{L} = (1 + z) D_\mathrm{M}$ is the luminosity distance.
Spatial distribution of mergers is determined by the SFR from \cite{2014ARA&A..52..415M} and delay time between system's birth and NS-WD coalescence from \cite{2018A&A...619A..53T}.

In these estimates we neglect absorption and redshift influence (i.e., shift in spectral energy distribution) because the intrinsic spectrum of transients is unknown.

\begin{figure}
\includegraphics[scale=0.9]{logNlogS.eps}
\caption{All-sky Log $N$ --- Log $S$ distributions for different peak
luminosities of NS-WD transients. S is the flux. The luminosity distribution is assumed to be a delta-function in each case. It is assumed that emission is collimated within one steradian. Flat parts of each curve correspond to $z> z_\mathrm{max}$.  }
\end{figure}

\begin{figure}
\includegraphics[scale=0.9]{L_dist.eps}
\caption{Critical luminosity versus distance.  Above the critical luminosity an NS-WD merger with the spectrum similar to the one of CDF-S XT2 produces more than two counts at eROSITA. The solid line shows results of calculations.  The dashed line is a power-law approximation. }
\end{figure}

\section{Expected number of NS-WD transients in the field of view of eROSITA and ART-XC}

eROSITA  is the primary instrument of the Spectrum-Roentgen-Gamma (SRG) mission successfully launched in July 2019.
The detailed description of the telescope is presented in \cite{merloni2012}.

The field of view (FoV) of eROSITA is 0.833 square degrees.
During 4 years of survey observations it will perform 8 all-sky surveys in soft X-rays within the energy range  $\sim 0.3 - 10\, \mathrm{keV}$, completing one full circle on the sky (one scan) every 4 hours.

In this section we estimate the number of NS-WD transients that fall in the field of view of eROSITA and ART-XC, the second telescope on-board SRG sensitive at the classical X-ray range from $\sim 3-5$~keV up to $\sim30$~keV.

The CDF-S XT2 X-ray outburst lasted for $\approx 20$ ks. However, after $t_\mathrm{break} = 2.3$ ks the flux started to decrease significantly. The $T_{90}$ parameter was estimated to be $\sim 11.1$ ks.
We therefore suggest that such transients last for $\sim 10$ ks.
During this time, telescope will cover a band of sky length $360^{\circ} \times { 10^4} / (4 \times 3600) = 250^{\circ}$.
This gives rate of
\begin{equation}
N = N_{\text{NSWD}} \times \frac {\text{(eROSITA FoV)}} {4 \pi (\frac{180}{\pi})^2 \text{ deg}^2} \times 250 \text{ deg} = 4.3
\end{equation}
events per day falling in the eROSITA FoV.

Similarly, for the ART-XC telescope with FoV $ = 0.3$ square degrees, the rate will be $\sim 1.5$ events per day.

This values are obtained without beaming. If we account for this effect then numbers are correspondingly reduced. Thus, we expect that after 4 years of survey with $\sim 1$ steradian beaming we expect $\gtrsim 500$
events in the field of view of eROSITA, and $\gtrsim 180$
--- for ART-XC.

\section{eROSITA rate of CDF-S XT2-like events}

In this section we model observations of CDF-S XT2-like events with eROSITA.

At the moment, there are no confirmed NS-WD coalescence detections. Accordingly, typical spectral parameters  of these events are unknown. However, in \cite{2019MNRAS.488..259F} the
authors claim that the X-ray transient CDF-S XT2 can be an example of such mergers. Therefore, it is reasonable to consider the spectrum of this event as a typical one in order to estimate how many objects eROSITA can observe and what are typical count rates.

In \cite{2019Natur.568..198X} the authors fit the spectrum of CDF-S XT2 with absorbed power-law with photon index $\Gamma = 1.57$ and hydrogen column density $N_\mathrm{H} = 7.7 \times 10^{21}$~cm$^{-2}$. Galactic contribution is fixed to a column density  $ 8.8 \times 10^{19}$~cm$^{-2}$.

To calculate the number of photons that can be detected by eROSITA from the source with such spectrum, we use the method described in \cite{2019AstL...45..120K}.
In this approach we consider interstellar absorption (data taken from \cite{morrison1983}), dependence of the effective area of the telescope on wavelength (these data is taken from the web-site https://wiki.mpe.mpg.de/eRosita), and redshift influence (shift of the spectral energy distribution). Luminosity (and spectrum) are assumed to be constant during the whole burst duration, as according to observations \cite{2019Natur.568..198X} flux evolves very slowly --- as $\propto t^{-0.14}$, --- during the first several thousand of seconds.

We assume that the dependence of hydrogen column density on intergalactic distance is linear:
\begin{equation}
    N_H = 8.8 \times 10^{19}\text{ cm}^{-2} + 7.7 \times 10^{21}\text{ cm}^{-2} \frac{D_L}{4687 \text{ Mpc}}.
\end{equation}

The resulting number of detected photons is
\begin{equation}
N = \int_{E_1}^{E_2} {  \Delta t \text{ }\frac{f(E) e^{-\sigma N_\mathrm{H}} S_\mathrm{eff}(E) } {4 \pi D_\mathrm{L}^2}} dE,
\end{equation}
where $E$~is the photon energy in keV, $f(E)$ is the spectrum of the source,
$ S_\mathrm{eff}(E)$ is the dependence of the effective area of the telescope on the wavelength, $\Delta t$ is the duration of the observation, $E_1$ and $E_2$ are the boundaries of the range of sensitivity of the telescope.

The host galaxy of CDF-S XT2 is located at $z = 0.738$ (which corresponds to $D_\mathrm{L}=4687$~Mpc), and for the corresponding distance the number of photons detected by eROSITA is more than one hundred for the whole time of the burst.
However, due to eROSITA observing strategy, such transient will be visible just in one scan of telescopes ($\sim 40$~s). In this case only 2 detected photons are expected.
This value is above the expected background, which is $\sim 10^{-2}$ photons/s/pixel \cite{merloni2012}, and so the transient can be marginally detectable.

Similar calculations for ART-XC give $\lesssim  1.6 \times 10^{-2}$ photons/s for the actual response matrix (which corresponds to the effective area $\sim 400$~cm$^{-2}$ \cite{2019ExA....47....1P}). The background for this instrument is currently unknown, but it is expected to be larger than for eROSITA. Thus, we assume that ART-XC will not be able to detect CDF-S XT2-like events, unless they happen at significantly smaller distance.

Taking into account similar technical characteristics of {XMM-Newton} (in comparison with eROSITA) and long operation time in pointing mode, it is quite probable that archival data of this instrument fits better requirements for searching for X-ray transients related to NS-WD mergers.

\section{Discussion}

\subsection{Uncertainties}

Degree of collimation of X-ray emission in transients under study is uncertain. Recently, in \cite{2019arXiv190710523D} the authors provided arguments that early X-ray afterglows after NS-NS coalescence can be unbeamed. Emission is provided by a newborn rapidly spinning magnetar. Such an object can be formed also in a NS-WD merger. Thus, similar arguments are applicable in this case, too. Then, the numbers given in Fig.~1 might be increased by the factor $4\pi$.

In some of our estimates (Sec. 3) we neglected absorption in interstellar and intergalactic medium. As NS-WD mergers typically happen at least hundreds of million years after formation of stars, a binary can travel far from the birth site, so interstellar absorption in the host galaxy can be negligible. In principle, intergalactic contribution can be significant in the eROSITA range of sen\-si\-ti\-vi\-ty (for ART-XС it is less important). We can estimate it as
$\frac{N_\mathrm{H}}{10^{22}\, \mathrm{cm}^{-2}}=\frac{DM}{333 \, \mathrm{pc\, cm}^{-3}}$, see \cite{2013ApJ...768...64H}, where $DM$ is dispersion measure, which for extragalactic sources (for example, Fast Radio Bursts --- FRBs) is about few hundreds for distances $\sim$ Gpc. Still, taking into account unknown spectral parameters (and other uncertainties) we prefer to make simple estimates with a toy model (see Fig. 2) in addition to calculations with parameters of CDF-S XT2 (see Sec. 5) which takes absorption into account.

Uncertainties in SFR at large redshifts (see \cite{2014ARA&A..52..415M})
might not influence our results drastically, as according to \cite{2018A&A...619A..53T}
the delay between binary formation and its coalescence is usually not very large: most of mergers happen after few Gyrs. Thus, the observable events might be dominated by sources born at $z\lesssim2$.
However, in these low redshifts the SFR is not known very precisely, too. And so, it is necessary to check how much this can influence our main results.

We studied the effect of uncertainty in the SFR applying an alternative value  presented in a recent e-print \cite{2019arXiv190710523D}:
 \begin{equation}
     \mathrm{SFR}(z) \approx 0.25 \text{ } e^{-[\ln{((1+z)/3.16)]^2 / 0.524}} \text{ M}_\odot \text{ Mpc}^{-3} \text{ yr}^{-1}.
 \end{equation}
 With this rate we obtained a higher number of NS-WD mergers: 1983
 per day across the sky assuming that emission is collimated within one steradian. This number might be compared with the rate of 843 events per day across the sky calculated above using the SFR given in eq.(1).
 This estimate makes the probability of detection of these transients with eROSITA (and/or ART-XC) much higher.

Uncertainties in the critical WD mass above which  transfer of matter to a NS is unstable are also not very important, as most of WDs have masses $>0.5$ M$_{\odot}$, signi\-fi\-cant\-ly above the recently derived limit $\gtrsim 0.2$ M$_\odot$, see \cite{2017MNRAS.467.3556B}.

The authors of \cite{2018A&A...619A..53T} made calculations only for the solar metallicity.
Of course, especially at large redshifts, metallicity can be different from this value. However, as we noticed above, the main contribution to the total rate is expected from sources born at $z \lesssim 2$,
where metallicity might be not very far from the solar value.
We expect that uncertainties due to this factor might be smaller than those from unknown features of binary evolution (common envelope parameters, etc.) and NS properties (kick velocity distribution, in the first place).

\subsection{NS-WD mergers and FRBs}

Recently in \cite{2018Ap&SS.363..242L} NS-WD mergers were considered as progenitors of FRBs. Indeed, these collisions have several attractive features which make dis\-cus\-sion of this possibility interesting. Below we briefly comment on this proposal on the base of expected number of transients.

Total rate of NS-WD mergers obtained above is $\sim$ 850 per day across the sky (or  $\sim$ 2000 with another SFR).
This is noticeably lower than several thousands per day --- the estimated rate of FRBs if weak events are taken into account (see reviews in \cite{2018PhyU...61..965P, 2019A&ARv..27....4P, 2019arXiv190605878C}).
Possible beaming reduces the number of observable transients related to NS-WD merger further (probably, by an order of magnitude).

As all other catastrophic scenarios NS-WD mergers by themselves cannot explain repeating bursts \cite{2019arXiv190803507T}.
However, a magnetar formed after coalescence potentially can be discussed as a source of subsequent flares. In this case we expect that numerous repeating events differ from the first main burst, which can be accompanied by counterparts in other wavelengths \cite{2019MNRAS.488..259F}.
In particular, NS-WD coalescences can be detectable by future gravitational wave detectors and in optical surveys as transient events.

Thus, we conclude that coalescences of NSs and WDs cannot explain significant fraction of FRBs. However, it is not excluded that some small percentage of these events can be due to NS-WD interaction (for example, those which are not related to sites of recent star formation activity). At the moment it is not clear if the population of FRBs is uniform or not. Drastic increase of statistics of these events due to such instruments as CHIME, ASKAP, MeerKAT, Apertif, etc. expected in near future might help to resolve this issue.

\section{Conclusions}

We estimated the total rate and distribution over redshift of NS-WD coalescences following calculations by \cite{2018A&A...619A..53T} and using starformation history from \cite{2014ARA&A..52..415M}.
The total number of events is $\sim 800-2000$ per day depending on the assumed model of star formation history, and the maximum of the distribution is at $z\approx 1.8$. Basing on this calculations we claim that the rate of NS-WD mergers is too low to explain significant fraction of FRBs.

We demonstrate that taking into account possible collimation of X-ray emission $\gtrsim 500$ events might appear in the field of view of eROSITA instrument aboard SRG satellite during 4 years of its survey program. Re\-spec\-ti\-ve\-ly, $\gtrsim 180$ are expected in the field of view of ART-XC --- the second instrument of the mission. Without collimation these numbers are an order of magnitude higher.

Expected numbers of counts are small because putative transients are expected to be visible just in one scan for tens of seconds, as the rotational period of the satellite is slightly longer than expected duration of a transient.

\section*{Acknowledgements}

We thank drs. Ekaterina Filippova and Galina Lipunova for useful discussions, and the referee for her/his comments.
The authors acknowledge the support from the Program of development of M.V. Lomonosov Moscow State University (Leading Scientific School ``Physics of stars, relativistic objects, and galaxies'').


\section{Introduction}

The quantum effects in low-dimensional frustrated magnets have been extensively studied in recent condensed matter physics.\cite{intfrust,diep} Various exotic quantum phases emerge from the interplay of quantum fluctuation and frustration. From a theoretical point of view, it is remarkable that the ground states of a class of models are written down exactly owing to the frustration. A well-known example is the Majumdar-Ghosh model whose ground state is a dimer state with a spontaneous breakdown of translational symmetry.\cite{mg} The diamond chain\cite{Takano-K-S,tsh,ht2017,kiku2,kiku3} is another example of a model with an exact ground state. Since the diamond chain has an infinite number of local conservation laws, the ground states can be classified analytically. For strong frustration, each ground state consists of a periodic array of spin clusters separated by singlet dimers owing to these local conservation laws. With the decrease of frustration, the length of each spin cluster increases. For weak frustration, various nonmagnetic or ferrimagnetic states are realized depending on the magnitudes of the constituent spins and/or the exchange parameters. \cite{Takano-K-S,tsh,ht2017}

In the real materials with a diamond-chain structure, the lattice distortion is inevitable. In addition to the well-known natural mineral azurite,\cite{kiku2,kiku3} which is a spin-1/2 distorted diamond chain, it has been recently pointed out that
 [Ni$_3$(OH)$_2$(O$_2$C-C$_2$H$_2$-CO$_2$)(H$_2$O)$_4$]$_2$H$_2$O\cite{guillou} can be regarded as a spin-1 distorted
 diamond chain.\cite{kuni,kiku1} Although the exchange parameters are fixed for real materials, even more exotic quantum states emerge in the whole parameter space of distorted diamond chains theoretically.

The ground states of spin-1 diamond chains(S1DC) without distortion have been investigated in Refs. \citen{Takano-K-S} and \citen{ht2017}. In addition to the spin cluster ground states, the nonmagnetic Haldane state and the ferrimagnetic state with spontaneous magnetization $m=1/6$ and 1/3 are found. A series of quantum phase transitions take place between these ground-state phases with different periodicities and with or without a spontaneous translational symmetry breakdown.

In the previous work\cite{hida2019}, we have investigated the effect of two types of lattice distortions that do not preserve the local conservation laws. In the present work, we consider the effect of the bond alternation depicted in Fig. \ref{lattice}. In contrast to the distortions discussed in Ref. \citen{hida2019}, the bond alternation preserves the same local conservation laws as the undistorted case.

This paper is organized as follows.
 {In \S 2}, the model Hamiltonian is presented.
 {In \S 3}, the ground-state phase diagram is determined analytically and numerically.
The last section is devoted to summary and discussion.

\section{Hamiltonian}

 The alternating-bond S1DC is described by the Hamiltonian
\begin{align}
{\mathcal H} = &\sum_{l=1}^{L} \left[(1+\delta)\v{S}_{l}(\v{\tau}^{(1)}_{l}+\v{\tau}^{(2)}_{l}) \right.
\nonumber\\
&+(1-\delta)(\v{\tau}^{(1)}_{l}+\v{\tau}^{(2)}_{l})\v{S}_{l+1}
+ \left. \lambda\v{\tau}^{(1)}_{l}\v{\tau}^{(2)}_{l}\right] ,
\label{hama}
\end{align}
where $\v{S}_{l}$ and $\v{\tau}^{(\alpha)}_{l}\ (\alpha=1,2)$ are spin 1 operators. The number of the unit cells is denoted by $L$, and the total number of sites $\Ns$ is $3L$.
 The Hamiltonian has three different exchange parameters, namely,
$1+\delta$, $1-\delta$, and $\lambda$. Here, $\delta$ represents the strength of the bond alternation, and $\lambda$ controls the frustration as depicted in Fig. \ref{lattice}.
In the case of $\delta = 0$, Eq.~(\ref{hama}) reduces to
the Hamiltonian of S1DC without bond alternation.\cite{Takano-K-S,ht2017}

\begin{figure}[t]
\centerline{\includegraphics[width=6cm]{69444Fig1.eps}}
\caption{Structure of an alternating-bond S1DC. $S=\tau^{(1)}=\tau^{(2)}=1$.
The exchange parameters are denoted by $\lambda$, $1+\delta$, and $1-\delta$,
where $\lambda$ and $\delta$ control the strengths of frustration and bond alternation, respectively.}
\label{lattice}
\end{figure}

The Hamiltonian (\ref{hama}) has a series of local conservation laws.
To see it, we rewrite Eq. (\ref{hama})
 in the form,
\begin{align}
\H &= \sum_{l=1}^{L} \H_l,\label{ham2}\\
\H_l&=(1+\delta)\v{S}_{l}\v{T}_{l}+(1-\delta)\v{T}_{l}\v{S}_{l+1}
+ \frac{\lambda}{2}\left(\v{T}^2_{l}-4\right),
\end{align}
where the composite spin operators $\v{T}_l$ are defined as
\begin{align}
\v{T}_{l} \equiv \v{\tau}^{(1)}_{l}+\v{\tau}^{(2)}_{l}
\quad (l = 1, 2, \cdots ,L).
\end{align}
Then, it is evident that
\begin{align}
[\v{T}_l^2, {\mathcal H}] = 0 \quad (l = 1, 2, \cdots , L).
\end{align}
Thus, we have $L$ conserved quantities $\v{T}_l^2$ for all $l$, even in the presence of bond alternation $\delta$.
By defining the magnitude $T_l$ of the composite spin $\v{T}_l$ by $\v{T}_l^2 = T_l (T_l + 1)$, we have a
 set of good quantum numbers $\{T_l; l=1,2,...,L\}$ where $T_l=$ 0,1, and 2.
The total Hilbert space of the Hamiltonian (\ref{ham2}) consists of
separated subspaces, each of which is specified by
a definite set of $\{T_l\}$, i.e., a sequence of 0, 1, and 2.
 A pair of spins with $T_l=0$ is a singlet dimer.
 A cluster including $n$ successive $T_l\neq 0$ pairs bounded by two dimers is called a cluster-$n$.

 Thus,  a cluster-$n$ is equivalent to the ground state of an alternating-bond antiferromagnetic Heisenberg chain consisting of $2n+1$ effective spins $\v{S}_l$ and $\v{T}_l$.
 We follow the terminology of Ref. \citen{tsh} to call the ground state consisting of an array of cluster-$n$'s separated by dimers as dimer-cluster-$n$ (DC$n$) phase. The phases without dimers are called DC$\infty$ phases.

\section{Ground-State Phase Diagram}

\subsection{Ground states for $\lambda>2$}

For $\lambda>2$, the state with $T_l=2$ is not allowed as in the case of undistorted S1DCs. This can be proven as follows:

Following Ref. \citen{Takano-K-S}, let us consider the ground state of the Hamiltonian of a single diamond $\H_l$
with $T_l=2$. According to the Lieb-Mattis theorem\cite{Lieb-Mattis}, the ground state is a singlet. By elementary manipulations, we can show that the following state $\ket{G}_l$ is the ground state of $\H_l$.
\begin{align}
\ket{G}_l&=\sqrt{\frac{1}{5}}\Big[\ket{-1}_{S_{l+1}}\ket{-1}_{S_l}\ket{2}_{T_{l}}
+\ket{1}_{S_{l+1}}\ket{1}_{S_l}\ket{-2}_{T_{l}}\nonumber\\
&+\sqrt{\frac{{2}}{{3}}}\ket{0}_{S_{l+1}}\ket{0}_{S_l}\ket{0}_{T_{l}}\nonumber\\
&-\frac{1}{\sqrt{2}}\left(\ket{-1}_{S_{l+1}}\ket{0}_{S_l}\ket{1}_{T_{l}}
+\ket{0}_{S_{l+1}}\ket{-1}_{S_l}\ket{1}_{T_{l}}\right.\nonumber\\
&\left.+\ket{0}_{S_{l+1}}\ket{1}_{S_l}\ket{-1}_{T_{l}}+\ket{1}_{S_{l+1}}\ket{0}_{S_l}\ket{-1}_{T_{l}}\right)
\nonumber\\
&+\frac{1}{\sqrt{6}}\ket{1}_{S_{l+1}}\ket{-1}_{S_l}\ket{0}_{T_{l}}+\frac{1}{\sqrt{6}}\ket{-1}_{S_{l+1}}\ket{1}_{S_l}\ket{0}_{T_{l}}\Big],
\end{align}
where $\ket{S_l^z}_{S_l}$ and $\ket{T_l^z}_{T_l}$ are the eigenstates of $S_l^z$ and $T_l^z$, respectively. In this state, we can show that
\begin{align}
\v{S}_l\v{T}_l&=\v{S}_{l+1}\v{T}_l=-3.\label{eq:st}
\end{align}
Hence, we find the ground-state energy $E_l$ of $\H_l$ with $T_l=2$ as
\begin{align}
E_l(T_l=2)&=-6 +\lambda.
\end{align}
Obviously, $E_l(T_l=0)=-2\lambda$. Since these are independent of $\delta$, the conclusion that the state with $T_l=2$ is not allowed in the ground state for $\lambda >2$ remains valid even for $\delta \neq 0$ following the proof of Lemma 1 of Ref. \citen{Takano-K-S}. Hence, the ground-state phase diagram for $\lambda >2$ is identical to that of the mixed diamond chain with $(S,\tau)=(1,1/2)$ discussed in Ref. \citen{hts_alt_mdc} as shown in Fig. \ref{phase}(a). The main features are summarized as follows:

The phase diagram consists of DC$n$ phases with finite $n$ and DC$\infty$ phases corresponding to the Haldane and dimer phases.  The DC$\infty$ phase is the Haldane phase or the dimer phase according to whether $\delta < \deltac\simeq 0.2598$ or $\delta > \deltac$, since the Hamiltonian (\ref{ham2}) reduces to the alternating-bond spin-1 Heisenberg chain.\cite{affleck1985,Affleck-Haldane1987,Kato-Tanaka1994,Yamamoto1994,Totsuka1995,Kitazawa-Nomura1997}  For $\delta=\deltac$ and $\delta \simeq 1$ the analytical results suggest the presence of an infinite series of phase transitions between DC$n$$(n=0, ...,\infty)$ phases with varying $\lambda$. Hence, we speculate that they should take place in the entire range of $\deltac \leq \delta \leq 1$. For $\delta=0$, only the DC$n$ phases for $n \leq 3$ and Haldane phase are realized.\cite{ht2017} For small $\lambda$, the ground state is the DC$\infty$ state that corresponds to the Haldane state for $\delta < \deltac$ and to the dimer state for $\delta > \deltac$.
\begin{figure}
\centerline{\includegraphics[width=7cm]{69444Fig2a.eps}}
\centerline{\includegraphics[width=7cm]{69444Fig2b.eps}}
\caption{(a) Overall ground-state phase diagram. The region $\lambda >2$ is reproduced from Ref. \citen{hts_alt_mdc}. Small open circles are the phase boundaries between DC$n$ and DC$(n-1)$ phases. For $n \geq 4$, these phase boundaries are only shown by solid lines to avoid complications. Other phase boundaries are shown by the large open circles and solid lines. The DC0-DC1 phase boundary calculated exactly in Ref. \citen{hts_alt_mdc} is shown by the thick dotted line. The open square is the accumulation point of infinite series of critical points at $\delta=\deltac$.  The filled squares are triple points. (b) The enlarged phase diagram in the small $\lambda$ region. The thick dotted line is the approximate relation (\ref{eq:decouple}) for $ \delta \simeq  1$}
\label{phase}
\end{figure}
The physical quantities characterizing these phases are low temperature magnetic susceptibility and residual entropy.\cite{hts_alt_mdc} The behaviors of these quantities are also the same as those  described in Sect. 5.3 of Ref. \citen{hts_alt_mdc}.

\subsection{Ground-state phase diagram for $\lambda <2$}
\subsubsection{Possible candidates of ground states}
\label{sec_cand}
As in the case of $\delta=0$,\cite{ht2017} we employ the finite size DMRG method to calculate the ground-state energies of the Hamiltonian (\ref{ham2}) with open boundary condition  up to $L=13$ for all possible configurations $\{T_l\}$. For $\lambda <2$, we find that the possible configurations that give the lowest energies are (a) $\forall l, T_l=1$, (b) $\forall l, (T_{2l},T_{2l+1})=(1,2)$ or (2,1), (c) $\forall l, T_{l}=2$, and the configurations that differ from them only locally up to the system size $L=13$. From this observation, we choose the above three configurations (a)-(c) as most plausible candidates of the configurations of $T_l$'s in the DC$\infty$ ground states.
We will also examine some of the states with longer periodicities numerically.
 For the case (a), the ground state is in the Haldane phase or the dimer phase according to whether $\delta < \deltac$ or $\delta > \deltac$\cite{affleck1985,Affleck-Haldane1987,Kato-Tanaka1994,Yamamoto1994,Totsuka1995,Kitazawa-Nomura1997}.
 Thus, we expect the following four different phases in the DC$\infty$ phase.

\begin{enumerate}
\item Haldane phase (H phase) :
$\{T_l\}=\{\dot{1}\}$ and $\delta < \deltac$.
\item Dimer phase (D phase) :
$\{T_l\}=\{\dot{1}\}$ and $\delta > \deltac$.
\item Ferrimagnetic phase with $M=\Ms/6$ (F$_{1/6}$ phase) :
$\{T_l\}= \{\dot{1}\dot{2}\}.$\label{f6}
\item Ferrimagnetic phase with $M=\Ms/3$ (F$_{1/3}$ phase) :
$\{T_l\}=\{\dot{2}\}$.
\end{enumerate}
Here, $\{\dot{T_1}\cdots\dot{T}_l\}$ denotes the configuration consisting of a periodic array of the sequence $T_1\cdots T_l$ over the whole chain. The values of the spontaneous magnetizations in the ferrimagnetic phases can be identified from the Lieb-Mattis theorem\cite{Lieb-Mattis} .

\subsubsection{Approximation for $\delta \simeq 1$}

We decompose the Hamiltonian (\ref{ham2}) into two parts
\begin{align}
\H&=\H_0+\H_1,
\end{align}
where
\begin{align}
\H_0&=\sum_{l=1}^L \H_{0l},\ \
\H_1=\sum_{l=1}^L \H_{1l},\\
\H_{0l}&=(1+\delta)\v{S}_l\v{T}_l+\frac{\lambda}{2}[T_l(T_l+1)-4], \\
\H_{1l}&=(1-\delta)\v{S}_{l+1}\v{T}_l.
\end{align}
We denote the ground state of $\H_{0l}$ with total spin $S^{\rm tot}_l$, where $\v{S}^{\rm tot}_l\equiv \v{S}_l+\v{T}_l$, as $\kket{S^{\rm tot}_l,S^{{\rm tot} z}_l}_{T_lS_{l}}$, and corresponding eigenvalue as $E_{0,l}(S^{\rm tot}_l)$. Depending on the values of $T_l (=1, 2)$, these are given as follows:
\begin{enumerate}
\item $T_l=2$:

The ground state is a triplet state that has total spin $S^{\rm tot}_l=1$. They are expressed in terms of the states $\ket{{S}^z_l}_{S_l}$ and $\ket{{T}^z_l}_{T_{l}}$ as.
\begin{align}
\kket{1,1}_{T_lS_{l}}&=\sqrt{\frac{3}{5}}\left(\ket{-1}_{S_l}\ket{2}_{T_{l}}
-\frac{1}{\sqrt{2}}\ket{0}_{S_l}\ket{1}_{T_{l}}\right.\nonumber\\
&+\left.\frac{1}{\sqrt{6}}\ket{1}_{S_l}\ket{0}_{T_{l}}\right),\label{eq:kket11}\\
\kket{1,0}_{T_lS_{l}}&=\sqrt{\frac{3}{10}}\Big(\ket{-1}_{S_l}\ket{1}_{T_{l}}
-\frac{2}{\sqrt{3}}\ket{0}_{S_l}\ket{0}_{T_{l}}\nonumber\\
&+\ket{1}_{S_l}\ket{-1}_{T_{l}}\Big),\label{eq:kket10}\\
\kket{1,-1}_{T_lS_{l}}&=\sqrt{\frac{3}{5}}\left(\ket{1}_{S_l}\ket{-2}_{T_{l}}
-\frac{1}{\sqrt{2}}\ket{0}_{S_l}\ket{-1}_{T_{l}}\right.\nonumber\\
&+\left.\frac{1}{\sqrt{6}}\ket{-1}_{S_l}\ket{0}_{T_{l}}\right),\label{eq:kket1m1}
\end{align}
and the corresponding energy is given by $E_{0l}(1)=-3(1+\delta)+\lambda$.

\item $T_l=1$

The ground state is a singlet state $\kket{0,0}_{T_lS_{l}}$ with a total spin $S^{\rm tot}_l=0$ and the corresponding energy is given by $E_{0l}(0)=-2(1+\delta)-\lambda$. The wave function is not presented since it is unnecessary for the following calculations.
\end{enumerate}

For $\delta=1$, the candidates of the ground states found in section \ref{sec_cand} can be constructed from the above local states as follows:

The $\{\dot{2}\}$ phase that corresponds to the ferrimagnetic ground state with magnetization $1/3$ per spin is given by
\begin{align}
\ket{F_{1/3}}=\prod_{l=1}^{L}\kket{1,1}_{T_lS_{l}},
\end{align}
if the spontaneous magnetization is taken along the $z$-axis. The energy of this phase is given by
\begin{align}
E_{1/3,0}/L=-3(1+\delta)+\lambda.
\end{align}
The $\{\dot{1}\}$ phase that corresponds to the nonmagnetic dimer ground state is given by
\begin{align}
\ket{D}=\prod_{l=1}^{L}\kket{0,0}_{T_lS_{l}}.
\end{align}
The energy of this phase is given by
\begin{align}
E_{D,0}/L=-2(1+\delta)-\lambda.
\end{align}

The $\{\dot{1}\dot{2}\}$ phase that corresponds to the ferrimagnetic ground state with magnetization $1/6$ per spin is given by
\begin{align}
\ket{F_{1/6}}&=\prod_{l=1}^{L/2}\kket{0,0}_{T_{2l-1}S_{2l-1}}\kket{1,1}_{T_{2l}S_{2l}},
\end{align}
or
\begin{align}
\ket{F_{1/6}}&=\prod_{l=1}^{L/2}\kket{1,1}_{T_{2l-1}S_{2l-1}}\kket{0,0}_{T_{2l}S_{2l}},
\end{align}
if the spontaneous magnetization is taken along the $z$-axis. The energy of this phase is given by
\begin{align}
E_{1/6,0}/L=-\frac{5}{2}(1+\delta).
\end{align}

Next, we calculate the energy of each ground state up to the first order with respect to the perturbation $\H_1$.
Among the terms in $\H_{1l}$, ${T}^{\pm}_l{S}^{\mp}_{l+1}$ has no diagonal elements within the ground state subspace of $\H_0$. Hence, we only consider the term ${T}^{z}_l{S}^{z}_{l+1}$.
By elementary manipulations, we find
\begin{align}
\bbra{1,1}_{T_lS_{l}}T^z_l\kket{1,1}_{T_lS_{l}}&
=\frac{3}{2},\\
\bbra{1,1}_{T_{l+1}S_{l+1}}S^z_{l+1}\kket{1,1}_{T_{l+1}S_{l+1}}&
=-\frac{1}{2}.
\end{align}
Hence, we have
\begin{align}
E_{1/3,1}&=\bra{F_{1/3}}\H_1\ket{F_{1/3}}=-\frac{3L}{4}(1-\delta).
\end{align}
Since $\ket{D}$ is nonmagnetic, expectation values of $T^z_l$ and $S^z_{l+1}$ vanish. Hence, we have
\begin{align}
E_{{\rm D},1}&=\bra{D}\H_1\ket{D}=0.
\end{align}

In $\ket{F_{1/6}}$, one of the groups of spins ($\v{S}_l$,$\v{T}_l$) or ($\v{S}_{l+1}$,$\v{T}_{l+1}$) connected by $\H_{1l}$ is always nonmagnetic. Hence, the expectation value of $T^z_lS^z_{l+1}$ vanishes. Thus, we find
\begin{align}
E_{1/6,1}&=\bra{F_{1/6}}\H_1\ket{F_{1/6}}=0.
\end{align}
Then, the ground state energies of the whole diamond chain in these three phases are given by
\begin{align}
E_{1/3}/L&=-3(1+\delta)-\frac{3}{4}(1-\delta)+\lambda,\\
E_{1/6}/L&=-\frac{5}{2}(1+\delta),\\
E_{\rm D}/L&=-2(1+\delta)-\lambda,
\end{align}
up to the first order in $\H_1$. The state $\ket{F_{1/6}}$ always has higher energy than the states $\ket{F_{1/3}}$ and $\ket{D}$ within the present approximation. Hence, it is not realized for $\delta \simeq 1$. Thus, the phase transition between the phases $\ket{F_{1/3}}$ and $\ket{D}$ takes place at
\begin{align}
\lambda_{\rm DF_{1/3}}(\delta)=1-\frac{1}{8}(1-\delta).
\label{eq:decouple}
\end{align}

\subsubsection{Numerical calculation}

Following the argument in Ref. \citen{ht2017}, for general values of $\delta$, the phase boundaries are determined from the ground-state energy per spin $\epsgtilde(\{T_l\},\delta)$ of an infinite-length Heisenberg chain with bond alternation $\delta$,
\begin{align}
\tilde{\H}=\sum_{l=1}^L \left[(1+\delta)\v{S}_l\v{T}_l+(1-\delta)\v{T}_l\v{S}_{l+1}\right],
\end{align}
with the corresponding configuration of $\{T_l\}$ in the following way:
\begin{enumerate}
\item Haldane-F$_{1/6}$ phase boundary for $\delta < \deltac$:
\begin{align}
 {\lambda}_{\rm HF_{1/6}}(\delta)&=2\epsgtilde(\{\dot{1}\},\delta) -2\epsgtilde(\{\dot{1}\dot{2}\},\delta).
\end{align}
\item Dimer-F$_{1/6}$ phase boundary for $\delta > \deltac$:
\begin{align}
 {\lambda}_{\rm DF_{1/6}}(\delta)
&=2\epsgtilde(\{\dot{1}\},\delta) -2\epsgtilde(\{\dot{1}\dot{2}\},\delta).
\end{align}
\item F$_{1/6}$-F$_{1/3}$ phase boundary:
\begin{align}
{\lambda}_{\rm F_{1/6}F_{1/3}}(\delta)
=2\epsgtilde(\{\dot{1}\dot{2}\},\delta)-2\epsgtilde(\{\dot{2}\},\delta).
\end{align}
\item D-F$_{1/3}$ direct transition phase boundary:
\begin{align}
 {\lambda}_{\rm DF_{1/3}}(\delta)
=\epsgtilde(\{\dot{1}\},\delta)-\epsgtilde(\{\dot{2}\},\delta).
\end{align}
\end{enumerate}
The energy $\epsgtilde(\{T_l\},\delta)$ can be calculated from the increment of the ground-state energy in each step of the iDMRG calculation.\cite{iDMRG,DMRG_rev} The phase diagram is shown in Fig. \ref{phase}. The number of states $\chi$ kept in each iDMRG step ranged from 240 to 480. Convergence with respect to $\chi$ is confirmed.

There are two triple points for $\lambda <2$. Among them, the D-F$_{1/3}$-F$_{1/6}$ triple point is $(\lambda,\delta)\simeq(1.03261,0.42495)$ and the D-H-F$_{1/6}$ triple point is  $(\lambda,\delta)=(\lambda_{\rm c},\deltac)$ where $\lambda_{\rm c}\simeq 1.0832$.

So far, we restricted ourselves to the candidate ground states with configurations (a)-(c). In addition to these most plausible candidates, we also examined the states with $\{\dot{1}...\dot{1}\dot{2}\}$ structures that have longer periodicities. We find these states are not ground states for almost all values of $\delta$ and $\lambda$. Nevertheless, in the close neighborhood of the point  $(\lambda_{\rm c},\deltac)$, we obtained numerical results suggesting the presence of the F$_{1/9}$ and  partial ferrimagnetic ground states\cite{pferri_s_h,furuya2014} in a tiny region of width $\Delta\lambda \sim 10^{-3}$. However, the complete analysis of this region requires extensive high precision numerical calculations for all possible configurations of $\{T_l\}$. Hence, we leave the full investigation of this region for future studies.

\section{Summary and Discussion}

 The ground-state phases of the alternating-bond S1DC are investigated. Owing to the local conservation laws, the ground states are rigorously constructed, once the ground states of corresponding Heisenberg chains are known. For $\lambda >2$, the local state $T_l=2$ is not allowed. Hence, the ground state is identical to that of the (1,1/2) alternating-bond mixed diamond chain discussed in Ref. \citen{hts_alt_mdc}. Each ground-state phase is described as a DC$n$ phase consisting of a uniform array of cluster-$n$'s that are equivalent to the ground states of spin-1 alternating-bond Heisenberg chains of length $2n+1$. The nonmagnetic DC$\infty$ phases that correspond to Haldane or dimer phases also appear for intermediate $\lambda$.

 For small $\lambda$, the ground states are the DC$\infty$ states. Especially for $\lambda <2$, the local state with $T_l=2$ is allowed. Hence, ferrimagnetic phases with $m=1/6$ and 1/3 are also realized in addition to the nonmagnetic Haldane or dimer phases. The ferrimagnetic phase with $m=1/6$ is accompanied by the spontaneous translational symmetry breakdown.

For $\delta \simeq 1$, we have carried out the perturbation analysis with respect to $1-\delta$. In this regime, only the dimer phase and ferrimagnetic phase with $m=1/3$ are realized and no spontaneous translational symmetry breakdown takes place.

In the close neighborhood of the point  $(\lambda_{\rm c},\deltac)$, we obtained numerical results suggesting the presence of a narrow partial ferrimagnetic phase. This might be related to the presence of the adjacent Haldane-dimer critical line with a Luttinger-liquid ground state along the line $\delta=\deltac$ and $\lambda > \lambda_{\rm c}$, since the partial ferrimagnetic phase can be regarded as a spontaneously magnetized Luttinger liquid\cite{furuya2014}.

Thus, we found that the introduction of the bond alternation to the S1DC gives rise to a rich variety of quantum phases and phase transitions. The experimental realization of the S1DC has been reported in Refs. \citen{kuni} and \citen{kiku1}, although the distortion pattern is different. We hope the rich phase diagram predicted in the present paper stimulates the synthesis of other S1DC materials in the near future.

\acknowledgments

We thank K. Takano for stimulating comments and discussion. The numerical diagonalization program is based on the package TITPACK ver.2 coded by H. Nishimori. A part of the numerical computation in this work has been carried out using the facilities of the Supercomputer Center, Institute for Solid State Physics, University of Tokyo and Yukawa Institute Computer Facility at Kyoto University.


\section{Introduction}

A hallmark of unconventional superconductivity is a nodal pairing
state where the excitation gap vanishes at points or lines in
momentum space~\cite{Mineev1999}.  Recently, however, a third type of
node has been proposed: extended Bogoliubov Fermi surfaces (BFSs)
where the excitation gap vanishes at a surface in momentum space
\cite{Agterberg2017, Brydon2018b}.  In clean, inversion-symmetric
(even-parity) superconductors that spontaneously break time-reversal
symmetry (TRS), all nodes are generically expected to be BFSs.
Crucial for the appearance of BFSs is that the superconductivity
involves more than one band: Specifically, the pairing between
electrons in different bands generates a pseudomagnetic field, which
``inflates'' point and line nodes of the intraband pairing potential
into BFSs. These nodal surfaces are robust against perturbations that
preserve particle-hole and inversion symmetries, which can be formulated
in terms of a $\mathbb{Z}_2$ topological
invariant~\cite{Kobayashi2014, Zhao2016, Agterberg2017, Bzdusek2017}.

A natural setting for the appearance of BFSs is in systems where a
multiband structure arises from the presence of discrete low-energy
electronic degrees of freedom apart from spin, e.g., atomic-orbital or
sublattice indices. This permits the construction of novel
``internally anisotropic'' pairing states where the Cooper-pair
wave function has nontrivial dependence upon the orbital or sublattice
indices~\cite{Nomoto2016b, Cheung2018, Brydon2018b}. Crucially, for the
appearance of BFSs, internally anisotropic pairing states are
typically characterized by both intraband and interband pairing
potentials~\cite{Brydon2018b}. Such pairing states have been proposed
for a wide variety of multiband systems of current interest, such as
iron-based superconductors~\cite{Gao2010, Nicholson2012, Ong2016,
  Vafek2016, AgterbergFeSe2017, Setty2019, Hu2019},
Cu$_x$Bi$_2$Se$_3$~\cite{Fu2014}, half-Heusler
compounds~\cite{Brydon2016, Agterberg2017, Savary2017, Yang2017,
  Roy2017, Timm2017, Boettcher2018a, Boettcher2018b, Kim2018,
  Brydon2018b}, the antiperovskite Sr$_{3-x}$SnO~\cite{Kawakami2018},
Sr$_2$RuO$_4$~\cite{Taylor2012}, UPt$_3$~\cite{Nomoto2016a,
  Yanase2016}, transition-metal dichalcogenides~\cite{Oiwa2018,
  Moeckli2018}, and twisted bilayer graphene~\cite{Guo2018, Su2018,
  FWu2019}. This long list of materials---some of which are believed
to support a time-reversal-symmetry-breaking (TRSB) state---is
encouraging for the existence of~BFSs.

Although BFSs are robust against symmetry-preserving perturbations,
this topological protection does not guarantee the existence of such a
state.  Instead, it is necessary to consider the thermodynamic
stability. Since a TRSB combination of two nodal pairing states
eliminates all nodes that are not common to both states, it is expected
to be energetically favored over time-reversal-symmetric
combinations~\cite{Sigrist1991}. This argument does not hold if the
resulting TRSB state possesses a BFS, however, as this implies a
nonzero density of states (DOS) at the Fermi energy, which, at first
glance, is unfavorable compared to the line nodes generic for
time-reversal-symmetric states.  It was argued in Ref.~\cite{Agterberg2017}
that a TRSB state with a BFS could, nevertheless, be energetically
favorable in the presence of sufficiently strong spin-orbit coupling
(SOC). This analysis was restricted to temperatures close to $T_c$,
however, and so did not account for the effect of the expected large
residual DOS at low temperatures. Moreover, although the TRSB state
becomes more stable with increasing SOC, the size of the BFS
decreases, as shown below.  It is, thus, unclear if BFSs can be realized
in a limit where they have a detectable effect on the electronic
structure~\cite{Lapp2019}.  Another interesting question raised by the
analysis in Refs.~\cite{Agterberg2017, Brydon2018b} is what happens at SOC
strengths insufficient for a stable TRSB state.

In this paper, we use mean-field theory to study the appearance of
BFSs in a paradigmatic model of a multiband system with strong SOC,
specifically the Luttinger-Kohn Hamiltonian of $j=3/2$ fermions in a
cubic material~\cite{LK56}. The $j=3/2$ degree of freedom naturally
leads to a multiband system and to internally anisotropic
pairing. Assuming pairing in a \textit{s}-wave $J=2$
channel~\cite{Brydon2016}, we construct the superconducting phase
diagram as a function of the SOC strength and temperature.  We focus
on a particular set of pairing states belonging to the irreducible
representation (irrep) $T_{2g}$, which is expected to provide a
typical picture for pairing in a higher-dimensional representation. At
vanishing SOC, a fully gapped time-reversal-symmetric superconducting
phase is realized as was predicted in Ref.~\cite{Ho1999}.  For nonzero
SOC, we obtain a rich phase diagram, which, in particular, contains a
sizable region with TRSB superconductivity. The largest BFSs that we
find lead to a residual zero-temperature DOS at the Fermi energy of
approximately 20\% of the normal-state DOS, which should leave clear
signatures in thermodynamic measurements. We also verify the existence
of a subdominant magnetic order parameter which is induced by the TRSB
superconductivity.

Our paper is organized as follows: In Sec.~\ref{sec:mft}, we
introduce our microscopic model and outline the mean-field theory,
including a discussion of previously known limits of vanishing and
strong SOC. We present the mean-field phase diagram in
Sec.~\ref{sec:trsb} and study the effect of cubic anisotropy of the
SOC. A key feature of the phase diagram is the first-order transition
into a time-reversal-symmetric superconducting state at intermediate
SOC strength, which we explain in terms of a simplified model.  This
is followed in Sec.~\ref{sec:trsbstate} by a detailed study of the
TRSB state and the induced magnetic order parameter. We summarize our
results and draw additional conclusions in Sec.~\ref{sec:conclusion}.

\section{Model and mean-field theory
\label{sec:mft}}

Our starting point is the Luttinger-Kohn Hamiltonian for $j=3/2$
fermions in a cubic material~\cite{LK56},
\begin{equation}
  \label{eq:LK}
  h(\bm{k}) = (\alpha \abs{\bm{k}}^2 - \mu)\, \openone_4 + \beta \sum_i k_i^2 J_i^2 + \gamma \sum_{i \neq j} k_i k_j J_i J_j ,
\end{equation}
where $i = x,y,z$ and $i+1 = y$ if $i = x$, etc., and $J_i$ are the
$4\times 4$ matrix representations of the angular momentum operators
$j = 3/2$. The $j=3/2$ fermions can arise due to the strong atomic
SOC, e.g., of spins $s=1/2$ and orbital angular momenta $l=1$ for
\textit{p} orbitals.  In addition to the spin-independent dispersion
coefficient $\alpha$ and the chemical potential $\mu$, the Hamiltonian
in Eq.~\eqref{eq:LK} includes the symmetry-allowed SOC terms
proportional to $\beta$ and $\gamma$. The Hamiltonian has doubly
degenerate eigenvalues given by
\begin{align}
  \label{eq:EV}
  \epsilon_{\bm{k},\pm}
  & = \biggl( \alpha + \frac{5}{4}\, \beta \biggr) \abs{\bm{k}}^2 - \mu \nonumber \\
  & \quad{} \pm \beta \sqrt{\sum_i \biggl[ k_i^4 + \biggl( \frac{3\gamma^2}{\beta^2} - 1 \biggr)
    k_i^2 k_{i+1}^2 \biggr]} .
\end{align}
Note that SOC lifts the fourfold degeneracy of the $j=3/2$ manifold
away from the $\Gamma$ point. Due to the presence of time-reversal and
inversion symmetries, the bands remain doubly degenerate so that the
states in each band can be labeled by a pseudospin-$1/2$
index~\cite{Brydon2018b}.

The description in terms of an effective spin $j=3/2$ permits Cooper
pairs with total angular momentum $J = 0$ (singlet) and $J = 1$
(triplet), but also $J = 2$ (quintet) and $J = 3$
(septet)~\cite{Boettcher2016, Brydon2016, Savary2017, Roy2017,
  Yang2017, Venderbos2018, Boettcher2018a, Boettcher2018b, Kim2018,
  Sim2019}. Similar to singlets and triplets, the quintet and septet
pairings correspond to even- and odd-parity orbital wave functions,
respectively. In particular, this allows for a broader variety of
\textit{s}-wave pairing states: Besides the usual singlet, there are
five additional quintet states with on-site pairing.

Restricting ourselves to such local pairing states, the pairing
interaction has the general form
\begin{equation}
  H_{\text{pair}} = \sum_{j} \sum_l \sum_{l_i\in l}
    V_l\, b^\dagger_{l_i,j}b_{l_i,j} ,
\end{equation}
where $b^\dagger_{l_i,j}$ creates a Cooper pair at site $j$ in channel
$l_i$ belonging to the irrep $l$~\cite{Brydon2016}. There are three
irreps in the cubic $O_h$ point group which support \textit{s}-wave
pairing: The singlet state belongs to the one-dimensional $A_{1g}$
irrep, whereas the five quintet states are distributed into the
two-dimensional $E_g$ and the three-dimensional $T_{2g}$ irreps.

Within the standard mean-field treatment, the interaction is decoupled
to obtain the effective single-particle Hamiltonian,
\begin{equation}
  \label{eq:H_MF}
  H_{\mathrm{MF}} = \sum_{\bm{k}} \biggl(
  \frac{1}{2}\, \Psi_{\bm{k}}^\dagger \mathcal{H}(\bm{k}) \Psi_{\bm{k}}
  + \sum_l \frac{\tr[\Delta_l \Delta_l^\dagger]}{V_l} \biggr),
\end{equation}
with the Bogoliubov--de~Gennes (BdG) Hamiltonian,
\begin{equation}
  \label{eq:H_BdG}
  \mathcal{H}(\bm{k}) =
  \begin{pmatrix}
    h(\bm{k})      & \Delta        \\
    \Delta^\dagger & -h^T(-\bm{k}) \\
  \end{pmatrix} ,
\end{equation}
and the Nambu spinors
$\Psi_{\bm{k}} = (c_{\bm{k}}, c_{-\bm{k}}^\dagger)^T$ with
$c_{\bm{k}} = (c_{\bm{k},3/2},\allowbreak c_{\bm{k},1/2},\allowbreak
c_{\bm{k},-1/2},\allowbreak c_{\bm{k},-3/2})^T$, where
$c_{\bm{k},\sigma}$ is the annihilation operator for a fermion with
momentum $\bm{k}$ and spin~$\sigma$.

In this paper, we focus on pairing states in the $T_{2g}$ irrep, where
a general pairing state can be written as
\begin{equation}
  \Delta = \Delta_0\,
  ( l_{yz} \eta_{yz}
  + l_{xz} \eta_{xz}
  + l_{xy} \eta_{xy}) ,
\end{equation}
with the amplitude $\Delta_0$, the three-component order parameter
$\bm{l} = (l_{yz},\allowbreak l_{xz},\allowbreak l_{xy})$, and the gap
matrices
$\eta_{\alpha\beta} = (J_\alpha J_\beta + J_\beta J_\alpha) U_T /
\sqrt{3}$, where
\begin{equation}
  U_T =
  \begin{pmatrix}
    0  & 0 & 0  & 1 \\
    0  & 0 & -1 & 0 \\
    0  & 1 & 0  & 0 \\
    -1 & 0 & 0  & 0 \\
  \end{pmatrix}
\end{equation}
is the unitary part of the time-reversal operator.  From the
fourth-order expansion of the corresponding Landau free energy, four
possible ground states are known: $\bm{l} = (1,0,0)$, $(1,1,1)$,
$(1,\ii,0)$, and $(1,\omega,\omega^2)$ with $\omega = \ee^{2\pi\ii/3}$
(as well as symmetry-related vectors) \cite{Sigrist1991}.  The states
$(1,0,0)$ and $(1,1,1)$ are time-reversal symmetric, whereas the
chiral state $(1,\ii,0)$ and the cyclic state $(1,\omega,\omega^2)$
break TRS and, therefore, support BFSs~\cite{Brydon2018b}.  In the
following, however, we focus on the submanifold of $T_{2g}$ states
spanned by the $\bm{l} = (1,0,0)$ and $(1,\ii,0)$ states by adopting
the mean-field ansatz:
\begin{equation}
  \Delta = \Delta_{yz} \eta_{yz} + \ii \Delta_{xz} \eta_{xz}
  , \label{eq:meanfield}
\end{equation}
with two \emph{real} variational parameters $\Delta_{xz}$ and
$\Delta_{xz}$. If one of the parameters is zero, we obtain the
TRS-preserving $\bm{l} = (1,0,0)$ state. On the other hand, a TRSB
state is realized if both parameters are nonzero; in particular, the
case of $\Delta_{yz}=\Delta_{xz}$ corresponds to $\bm{l} = (1,\ii,0)$.
Although this restricted ansatz is artificial for a cubic system, we
are motivated by the observation that the $\eta_{xz}$ and $\eta_{yz}$
pairing potentials are the only \textit{s}-wave quintet states in our
cubic model which are also degenerate in hexagonal and tetragonal
crystals.  For example, a chiral \textit{d}-wave state with the same
symmetry is believed to be realized in tetragonal
URu$_2$Si$_2$~\cite{Kasahara2007}. We, therefore, expect our conclusions
to be applicable to any TRSB superconductor with two degenerate
pairing potentials.  The $(1,\ii,0)$ state has an (inflated)
equatorial line node, which should lead to a higher free energy
compared to a state with only (inflated) point nodes.  By considering
the likely less favorable pairing state, we, at worst, underestimate
the stability of the BFSs.  In fact, performing the same analysis for
the pair of $E_g$ states does not result in qualitative changes in the
phase diagram.  The pairing state in the spherically symmetric limit
has been considered in Refs.~\cite{Boettcher2016, Boettcher2018a,
  Herbut2019}.

\subsection{Free energy}

In a weak-coupling approach, the leading pairing instability can be
obtained by direct minimization of the Helmholtz free energy with
respect to the mean fields.  From the BdG Hamiltonian \eqref{eq:H_MF},
we obtain the Helmholtz free energy
\begin{equation}
\label{eq:Helmholtz}
  F = \sum_{\bm{k}} \frac{\tr[\Delta \Delta^\dagger]}{V_0}
  - 2 k_B T \sum_{\bm{k},\nu} \ln\biggl[ 2 \cosh\biggl(\frac{E_{\bm{k},\nu}}
    {2 k_B T} \biggr) \biggr] ,
\end{equation}
where $V_0$ is the attractive pairing interaction in the $T_{2g}$
channel and $E_{\bm{k},\nu}$ are the positive eigenvalues of
$\mathcal{H}(\bm{k})$ in Eq.~\eqref{eq:H_BdG}.  Inserting the
mean-field ansatz from Eq.~\eqref{eq:meanfield}, we numerically
minimize the Helmholtz free energy to obtain the self-consistent
values of $\Delta_{xz}$ and $\Delta_{yz}$.

To compare with our numerical calculation and previous
results~\cite{Agterberg2017}, we also use the complementary approach
of expanding the free energy in the pairing potential to obtain the
Ginzburg-Landau (GL) free energy~\cite{Dahl2007},
\begin{subequations}
\label{eq:GL}
\begin{equation}
  \label{eq:GL-F}
  F = \sum_{\bm{k}} \frac{\tr[\Delta \Delta^\dagger]}{V_0}
  + k_B T \sum_{\bm{k},\,\ii \omega_n}
    \sum_{l=1}^\infty \frac{1}{l} \tr[ ( G_0 \Sigma )^l ] ,
\end{equation}
with
\begin{equation}
  \label{eq:G0-Sigma}
  G_0 =
  \begin{pmatrix}
    G(\bm{k},\ii\omega_n) & 0 \\
    0 & \tilde{G}(\bm{k},\ii\omega_n) \\
  \end{pmatrix}
  \quad\text{and}\quad
  \Sigma =
  \begin{pmatrix}
    0 & \Delta \\
    \Delta^\dagger & 0 \\
  \end{pmatrix}
  ,
\end{equation}
\end{subequations}
where $G$ and $\tilde{G}$ are the particlelike and holelike Green's
functions of the normal-state Hamiltonian $h(\bm{k})$ and
$\ii \omega_n = \ii (2n+1)\pi k_B T$ are the fermionic Matsubara
frequencies.  For this choice of $\Sigma$, all terms with odd $l$
vanish.  The GL free energy can be evaluated analytically, see
Appendix~\ref{sec:GL} for an example calculation and the necessary
approximations.

\subsection{Known limits}

Previous work has revealed the behavior of the model in the limiting
cases of vanishing and strong SOC \cite{Ho1999, Brydon2016,
  Agterberg2017, Brydon2018b, Venderbos2018}. We summarize the results
in the following.

\subsubsection{Vanishing spin-orbit coupling}

The case of vanishing SOC was studied by \citet{Ho1999} in the context
of pairing in fermionic cold atomic gases.  They found that for
\textit{s}-wave quintet pairing, a TRS-preserving state is
energetically favored compared to a TRSB state. To understand this
limit, we first note that the vanishing SOC implies that the
eigenvalues of the normal-state Hamiltonian are fourfold
degenerate. As such, the pairing potential and the normal-state
Hamiltonian can be simultaneously diagonalized by a
momentum-independent spin rotation.  The resulting eigenvalues are
identical to the case of a \textit{s}-wave singlet gap and so the gap
is uniform across the Fermi surface.  For TRSB pairing states, two of
the diagonal entries of the diagonalized pairing potential vanish,
indicating that two of the four degenerate Fermi surfaces remain
ungapped in the superconducting state. On the other hand, a
TRS-preserving state opens a gap on all the Fermi surfaces and is,
thus, energetically favorable.  In real materials, a nonzero SOC is
always present, which lifts the fourfold degeneracy.  We, nevertheless,
expect that for sufficiently weak SOC, the time-reversal-symmetric
state proposed by \citet{Ho1999} persists.

\subsubsection{Strong spin-orbit coupling}
\label{sss:largeSOC}

In the limit where the SOC-induced splitting of the bands is much
larger than the pairing potential, an effective single-band model can
be used for the states close to the Fermi
energy~\cite{Brydon2016}. Specifically, we write the effective BdG
Hamiltonians for the two bands labeled by $\pm$ in the pseudospin
basis as
\begin{equation} \label{eq:Heff}
  \mathcal{H}_{\text{eff},\pm}(\bm{k}) = \begin{pmatrix}
    \epsilon_{\bm{k},\pm}s_0 + \delta H_{\bm{k},\pm} & \pm
    \psi^{\text{intra}}_{\bm{k}}is_y \\
    \mp \psi^{\text{intra}\,\ast}_{\bm{k}}is_y & -\epsilon_{\bm
      k,\pm}s_0 - \delta H^T_{-\bm{k}, \pm}
    \end{pmatrix} ,
\end{equation}
where $s_\mu$ are the Pauli matrices in the pseudospin space. The
effective Hamiltonian describes intraband pseudospin-singlet pairing
with potential,
\begin{equation}
  \label{eq:psieff}
  \psi_{\bm{k}}^{\text{intra}} = \frac{\sqrt{3}\gamma}{2}\,
  \frac{\Delta_{yz} k_y k_z + i\Delta_{xz} k_x k_z}{\sqrt{\sum_i \bigl[ \beta^2 k_i^4
      + (3\gamma^2 - \beta^2) k_i^2 k_{i+1}^2 \bigr]}} .
\end{equation}
The interplay of the quintet pairing with the normal-state spin-orbit
texture gives the intraband potential a \textit{d}-wave form factor,
reflecting the $J=2$ total angular momentum of the Cooper pairs and
imposes a sign difference between the bands. The nodal structure
of the intraband potential favors a TRSB combination of the quintet
states as this gaps out nonintersecting line nodes thereby
enhancing the average gap magnitude and, thus, lowering the free
energy~\cite{Sigrist1991}.  Since the $\eta_{yz}$ pairing potential
leads to line nodes on the $k_y=0$ and $k_z=0$ planes, whereas the
$\eta_{xz}$ state has line nodes on the $k_x=0$ and $k_z=0$ planes,
the $\bm{l} = (1,i,0)$ state is characterized by point nodes along the
$k_z$ axis and a line node on the $k_z=0$ plane.

The diagonal blocks of the effective BdG Hamiltonian in
Eq.~\eqref{eq:Heff} obtain a correction term $\delta H_{\bm{k},\pm}$
from including the effect of interband pairing to second order in
perturbation theory~\cite{Agterberg2017, Brydon2018b,
  Venderbos2018}. This correction has the general form
\begin{equation}
\delta H_{\bm{k},\pm}=  \gamma_{\bm{k},\pm} s_0 + \bm{h}_{\bm{k},\pm}
  \cdot \bm{s} ,
\end{equation}
where $\gamma_{\bm{k},\pm}$ renormalizes the band dispersion and is
always nonzero in the presence of interband pairing, whereas
$\bm{h}_{\bm{k},\pm}$ describes an effective pseudomagnetic field that
is only present for TRSB states. The two contributions can be written
as
\begin{align}
    \gamma_{\bm{k},\pm} &= \frac{1}{2 (\epsilon_{\bm{k},+} - \epsilon_{\bm{k},-})}\,
  \tr[\mathcal{P}_{\bm{k},\pm} \Delta\Delta^\dagger \mathcal{P}_{\bm{k},\pm}] ,\\
  \label{eq:pseudomagn}
  \bm{h}_{\bm{k},\pm} &= \frac{1}{2 (\epsilon_{\bm{k},+} - \epsilon_{\bm{k},-})}\,
  \tr[\bm{s} \mathcal{P}_{\bm{k},\pm} \Delta\Delta^\dagger \mathcal{P}_{\bm{k},\pm}] ,
\end{align}
where $\mathcal{P}_{\bm{k},\pm}$ are projection operators on the
normal-state Hilbert spaces of the $\pm$ bands.  The pseudomagnetic
field is crucial for the appearance of BFSs as can be seen from the
dispersion in the effective low-energy model,
\begin{equation}
  \label{eq:Eeff}
  E_{a,b,\pm} = a \abs{\bm{h}_{\bm{k},\pm}}
    + b \sqrt{[\epsilon_{\bm{k},\pm} + \gamma_{\bm{k},\pm}]^2
    + \abs{\psi^{\text{intra}}_{\bm{k}}}^2} ,
\end{equation}
where $a$ and $b$ are independently chosen to be $\pm 1$, giving four
bands. In the absence of the pseudomagnetic field, a node occurs where
the square root vanishes, but the pseudomagnetic field is generally
nonzero at these momenta. This lifts the pseudospin degeneracy by
shifting the pseudospin-up and pseudospin-down bands in opposite
directions and leads to the formation of BFSs \cite{Agterberg2017,
  Brydon2018b}.  Although this increases the free energy of the TRSB
state, for sufficiently small $\abs{\bm{h}_{\bm{k},\pm}}$, it should
not cause a transition to a TRS-preserving phase, since the energy
difference between the lowest TRSB and TRS-preserving states is
generically finite. In particular, from Eq.~\eqref{eq:pseudomagn} we
expect that a TRSB state with BFSs is stable for
$\abs{\Delta_{yz}}, \abs{\Delta_{xz}} \ll
\abs{\epsilon_{\bm{k},+}-\epsilon_{\bm{k},-}}$.

\section{Phase diagram}
\label{sec:trsb}

\begin{figure}[tb]
  \centering
  \includegraphics{fig01.pdf}
  \caption{Representative dispersion relations for the spherically
    symmetric model described by Eq.~\eqref{eq:hspherical}. (a)
    Without SOC, the bands are fourfold degenerate, and there is a
    single Fermi surface with wave vector
    $k_{F,0}=\sqrt{\mu/\alpha}$. Note that we take $\alpha>0$ so that
    the band has positive effective mass.  (b) Turning on the SOC
    lifts the fourfold degeneracy, yielding two doubly degenerate
    quadratically dispersing bands with positive effective mass. There
    are now two Fermi surfaces with wave vectors
    $k_{F,\pm} = \sqrt{\mu/(\alpha + 5\beta/4 \pm \beta)}$. (c) For
    $\beta < -4\alpha/9$, the effective mass of one of the bands
    becomes negative, and there is only a single Fermi surface.}
  \label{fig:cartoon}
\end{figure}

We start by considering the case of spherically symmetric SOC, i.e.,
$\beta = \gamma$, and later generalize to the case of cubic
anisotropy.  In the spherical limit, the normal-state Hamiltonian
simplifies to
\begin{equation}
  h(\bm{k}) = (\alpha \abs{\bm{k}}^2 - \mu)\, \openone_4 + \beta (\bm{k}
  \cdot \bm{J})^2 . \label{eq:hspherical}
\end{equation}
Representative examples of the normal-state band structure are shown
in Fig.~\ref{fig:cartoon}.

In Fig.~\ref{fig:spher}, we present the phase diagram as a function of
temperature and SOC strength.  Figures~\ref{fig:insets}(a)--\ref{fig:insets}(f) show
the band structure around the Fermi surface in the $[100]$ direction
where we anticipate the appearance of nodes from the projected gap in
Eq.~\eqref{eq:psieff}. Any gaps in the spectrum along this direction
at nonzero SOC strength are, therefore, due entirely to the interband
pairing potential.  To obtain comparable results over a wide range of
SOC strengths, we fix the critical temperature $T_c$ and vary the
attractive interaction $V_0$ such that the second-order coefficient of
the GL free energy vanishes at the chosen $T_c$. This eliminates
effects due to the changing DOS at the Fermi energy as the SOC is
varied.

\begin{figure*}[tb]
  \centering
  \includegraphics{fig02.pdf}
  \caption{Phase diagram for the $T_{2g}$ pairing states given by
    Eq.~\eqref{eq:meanfield} on the SOC-temperature plane.  The color
    code indicates the gap magnitude
    $\sqrt{\Delta_{xz}^2+\Delta_{yz}^2}$ where brighter colors mean
    larger gaps and white means no superconductivity.  The horizontal
    line at $T/T_c = 1$ denotes the critical temperature $T_c$
    predicted by GL theory.  Lines of first-order (second-order) phase
    transitions are indicated in red (orange).  The blue dot in both
    panels indicates the point of TRSB from GL theory, the red dot in
    the panel on the right denotes the onset of the first-order phase transition
    estimated by GL theory. The left panel is a zoom of the box in the
    right panel.  The SOC strength $\beta$ is plotted as an effective
    spin-orbit energy $\beta k_F^2/(k_B T_c)$ where
    $k_F^2 = \mu/(\alpha + 5\beta/4)$.}
  \label{fig:spher}
\end{figure*}

\begin{figure*}[tb]
  \includegraphics{fig03.pdf}
  \caption{(a)--(f) Band structure in the vicinity of the Fermi energy
    for parameter sets indicated by the corresponding labels in
    Fig.~\ref{fig:spher} along the $[100]$ direction where we expect
    nodes in a nodal state.  The Fermi wave vectors are given by
    $k_{F,\pm} = \sqrt{\mu/(\alpha + 5\beta/4 \pm \beta)}$.}
  \label{fig:insets}
\end{figure*}

Starting at $\beta = 0$, we find the fully gapped TRS-preserving state
(``nodeless TRS'') predicted by \citet{Ho1999}.  Switching on the SOC,
we observe that the gap just below the critical temperature has nodes
(``nodal TRS''), but the nodeless TRS state is recovered at lower
temperatures. The nodal behavior arises as the SOC lifts the fourfold
degeneracy of the bands, making a distinction between inter- and
intraband pairings possible. Close to the critical temperature, the
strength of the pairing potential is much smaller than the band
splitting so that the gap at the Fermi surface is controlled by the
nodal intraband pairing potential in Eq.~\eqref{eq:psieff}.  However,
as the pairing potential grows upon lowering the temperature, the
interband potential shifts the nodes away from the Fermi surfaces at
$k_{F,\pm} = \sqrt{\mu/(\alpha + 5\beta/4 \pm \beta)}$ as seen in the
band structure at point~(f) in Fig.~\ref{fig:insets}.  At a critical
value of the pairing potential, the nodes meet and annihilate, marking
the recovery of the nodeless TRS phase.

A further increase in SOC leads to an enhancement of the critical
temperature over the one anticipated from the second-order coefficient
of the GL free energy, implying a first-order transition between the
normal and the superconducting states.  The presence of the first-order
phase transition is confirmed by computing the position of the
tricritical point from GL theory, i.e., the point where the
fourth-order coefficient turns negative.  We find very good agreement
between the numerical calculation and our GL theory (cf.\ the red dot
in the right panel of Fig.~\ref{fig:spher}). In this region, the
magnitude of the pairing potential
$\sqrt{\Delta_{xz}^2+\Delta_{yz}^2}$ is much larger than expected from
BCS theory, and the very large interband pairing potentials ensure a
full gap as shown by points~(d) and~(e) in Fig.~\ref{fig:spher}. We,
hence, refer to this state as the ``large-gap'' phase, in contrast to
the other, ``small-gap'' phases. The origin of the first-order phase
transition is discussed in Sec.~\ref{sec:toy-model}.

\begin{figure}
  \centering
  \includegraphics{fig04.pdf}
  \caption{Sketch of the pairing amplitudes $\abs{\Delta_{xz}}$ and
    $\abs{\Delta_{yz}}$ and of the gap structure in the nodal TRS, TRSB
    $C_2$, and TRSB $C_4$ phases, see Fig.~\ref{fig:spher}. The TRSB
    $C_2$ state breaks both TRS and $C_4$ symmetry.}
  \label{fig:phase-cartoon}
\end{figure}

Upon increasing the SOC strength beyond
$\beta k_F^2 \approx -8.4 k_B T_c$, there is an abrupt drop in the
magnitude of the gap, and the nodeless TRS phase gives way to a nodal
state.  Close to $T_c$, this state marked by ``nodal TRS'' has a gap
that is well approximated by Eq.~\eqref{eq:psieff} and exhibits line
nodes [point~(e) in Figs.~\ref{fig:spher}
and~\ref{fig:insets}]. Further below $T_c$, we enter a phase which
breaks TRS but where the two gap parameters $\Delta_{yz}$ and
$\Delta_{xz}$ have unequal magnitude.  We label this the ``TRSB
$C_2$'' state because the unequal gap magnitudes yield a spectrum with
only $C_2$ rotational symmetry about the $z$ axis.  The magnitudes of
$\Delta_{yz}$ and $\Delta_{xz}$ converge as the SOC is increased, thus,
realizing the ``TRSB $C_4$'' state where the spectrum has $C_4$
rotational symmetry about the $z$ axis.  This is the
$\bm{l} = (1,i,0)$ state and is consistent with predictions of the
strong-SOC limit discussed in Sec.~\ref{sss:largeSOC}.  The
intermediate TRSB $C_2$ state can be visualized as a continuous
rotation of the vector $\bm{l}$ from $(1,0,0)$ to $(1,i,0)$, see
Fig.~\ref{fig:phase-cartoon}.  The boundary of the TRSB $C_4$
phase shows reentrant behavior, but it is realized at all temperatures
for $\beta k_F^2 \approx -9.7 k_B T_c$.  Both the TRSB $C_4$ and the
$C_2$ phases display BFSs.

The critical value of the SOC strength for which the TRSB state
becomes stable just below $T_c$ is estimated from an expansion of the
GL free energy to fourth order at
$\beta k_F^2 \approx -8.957 k_B T_c$.  This estimate is shown as the
blue dot in both panels of Fig.~\ref{fig:spher} and is in excellent
agreement with the mean-field theory.  A previous
analysis~\cite{Agterberg2017} had estimated this critical strength to
be $\beta k_F^2 \approx -11.572 k_B T_c$ (expressed in our units) and,
therefore, overestimated it by about $30\%$. The disagreement stems
from the approximate treatment of the band splitting in
Ref.~\cite{Agterberg2017}. Nevertheless, we confirm that the TRSB
state is realized at moderate values of the SOC strength.

\subsection{Effects of cubic anisotropy}
\label{sec:cubic}

Cubic anisotropy is introduced in our model by setting
$\gamma \neq \beta$ in the Luttinger-Kohn Hamiltonian. In
Fig.~\ref{fig:F4}, we show the pairing state realized just below the
critical temperature as a function of $\beta$ and the cubic anisotropy
parameter $\gamma-\beta$. Note that the transition into the large-gap
phase is of first order and the critical temperature, therefore,
exceeds the temperature at which the second-order coefficient in the
GL expansion changes sign. As can be seen, there is a pronounced
asymmetry between the cases of $\abs{\gamma}>\abs{\beta}$ and
$\abs{\gamma}<\abs{\beta}$: The region of first-order transitions into
the large-gap phase is suppressed for $\abs{\gamma}>\abs{\beta}$ and
disappears entirely for sufficiently strong $\gamma$, and the TRSB
state occurs at smaller values of the SOC strength
$\abs{\beta}$. These trends are reversed for $\gamma>\beta$.

\begin{figure}[tb]
  \includegraphics{fig05.pdf}
  \caption{Phase diagram just below the critical temperature as a
    function of SOC strength and cubic anisotropy.  ``$(1,0,0)$
    small'' and ``$(1,0,0)$ large'' refer to the TRS-preserving phase,
    whereas ``$(1,i,0)$'' specifies the TRSB $C_4$ phase with BFSs.}
  \label{fig:F4}
\end{figure}

In Fig.~\ref{fig:cubic}, we show temperature-dependent phase diagrams
along two lines $\gamma = 2\beta$ and $\gamma = \beta/2$.  Along the
cut $\gamma = 2\beta$, there is no first-order phase transition.  The
change in gap magnitude along the nodeless to nodal transition is
steep but not abrupt.  This transition is accompanied by the
disappearance of nodes.  The intermediate $C_2$ phase is also heavily
suppressed.  Along the other cut $\gamma = \beta/2$, we do not recover
the small-gap phase within the boundaries of the graph.  Therefore, we
also do not observe the point of TRSB as predicted by GL theory.

\begin{figure}[t!b]
  \includegraphics{fig06.pdf}\\[2ex]
  \includegraphics{fig07.pdf}
  \caption{Phase diagrams on the SOC-temperature plane along two
    lines (a) $\gamma = 2\beta$ and (b) $\gamma = \beta/2$.  The
    color code represents the gap magnitude where brighter colors
    mean larger gaps and white means no superconductivity.  Lines of
    first-order (second-order) phase transitions are indicated in red
    (orange).  For $\gamma=2\beta$, panel (a), there is no first-order
    phase transition at $T_c$.  Below $T_c$, we find a large-gap phase,
    but the transition to it is not of first order.  For
    $\gamma = \beta/2$, panel (b), the large-gap phase occurs.  The
    critical temperature is strongly enhanced, and larger gaps are
    found.  In both panels, the blue dot at $T_c$ is the point of TRSB,
    and the red dot at $T_c$ is the tricritical point as predicted by
    the GL free energy.}
  \label{fig:cubic}
\end{figure}

The phase diagram in Fig.~\ref{fig:F4} can be understood by looking at
the expression for the effective intraband pairing
Eq.~\eqref{eq:psieff}.  We find that the magnitude of the intraband
pairing is proportional to $\gamma$, i.e., larger (smaller) $\gamma$
means stronger (weaker) intraband pairing compared to interband
pairing. The existence of the large-gap phase depends on the ratio
between intra- and interband pairings as we will discuss in the next
section.

\subsection{Origin of the first-order transition}
\label{sec:toy-model}

The first-order phase transition into the large-gap phase shown in
Figs.~\ref{fig:spher}, \ref{fig:F4}, and \ref{fig:cubic} is one of the
most surprising features of the phase diagram of our model. The
inclusion of cubic anisotropy reveals that it is not generic, however,
but rather depends upon the balance between the two spin-orbit terms.
In this section, we show that the first-order transition is controlled
by the relative strengths of the intra- and interband pairing
potentials, which, in turn, depends on the SOC strengths as noted
above.

The first-order transition can be understood based on a simplified
model with two bands in which we fix the ratio of the inter- and
intraband pairing potentials. In this model, the normal-state bands
have the dispersions
$\xi_{\bm{k},\pm} = (1\pm \delta)\epsilon_{\bm{k}} - \mu$, where
$\delta$ parametrizes the band splitting and the precise form of
$\epsilon_{\bm{k}}$ is unimportant. The splitting parameter $\delta$
plays a role analogous to the SOC strength in the full model where
the band splitting is characterized by differing effective masses of
the Luttinger-Kohn bands as illustrated in Fig.~\ref{fig:cartoon}.

Since the interband and intraband pairing potentials are obtained by
projecting $\Delta$ from Eq.~\eqref{eq:meanfield} into the band basis,
the relative strength of the interband and intraband pairing
potentials is determined by details of the normal-state band
structure.  To represent this aspect, we write the pairing potential
in the band basis as
\begin{equation}
  \label{eq:psitoy}
  \Delta = \eta
  \begin{pmatrix}
    r            & \sqrt{1-r^2} \\
    \sqrt{1-r^2} & -r           \\
  \end{pmatrix},
\end{equation}
where $\eta$ is the magnitude of the pairing potential, and the
coefficient $r$ controls the relative strength of intra- and interband
pairings: $r = 0$ corresponds to pure interband pairing and $r = 1$ to
pure intraband pairing. The intraband pairing has opposite signs in
each band, in agreement with Eq.~\eqref{eq:psieff}. Since the
first-order transition only occurs into a TRS-preserving state, in the
following we assume that $r$ and $\eta$ are real.  Note that for the
pairing potential in Eq.~\eqref{eq:psitoy}, the ratio between intra-
and interband pairings is momentum independent. In contrast, in the
full model, this quantity varies across the Fermi surface.  We can
nevertheless define this ratio for the full model in terms of the
Fermi-surface average,
\begin{equation} \label{eq:avr}
  r^2 = \frac{1}{\Delta_{xz}^2+\Delta_{yz}^2}
  \int \frac{\diff\Omega}{4\pi}\:\abs{\psi_{\bm{k}}^{\text{intra}}}^2 .
\end{equation}

\begin{figure}[tb]
  \includegraphics{fig08.pdf}
  \caption{Phase diagram of the simple model as a function of the
    pairing ratio $r$ and the band splitting $\delta$.  In this model,
    first-order phase transitions are only possible for $x > x_c$
    because only then there is a region where $F_4 < 0$.}
  \label{fig:toy}
\end{figure}

The GL expansion of the free energy of the simple model gives a Taylor
series in the parameter $\eta$,
\begin{equation} \label{eq:GLexpand}
  F = F_2\eta^2 + F_4\eta^4 + \mathcal{O}(\eta^6) ,
\end{equation}
where expressions for the coefficients $F_2$ and $F_4$ can be obtained
from Eqs.~\eqref{eq:GL-F} and~\eqref{eq:G0-Sigma}. A negative sign of
the fourth-order coefficient $F_4$ indicates that the transition into
the superconducting state is of first order. We show the variation of
the sign of this coefficient as a function of the parameters $r$ and
$\delta$ in Fig.~\ref{fig:toy}. For sufficiently small intraband
pairing strength $r$, we find that $F_4$ is positive at small band
splitting $\delta$, becomes negative for increasing $\delta$, and
finally returns to a positive value. Assuming that higher-order terms
in the GL expansion can be ignored, this indicates that the phase
transition becomes discontinuous beyond a critical band splitting, but
a continuous transition is recovered as the band splitting is further
increased.

The conclusions for the simple model are broadly in agreement with the
phase diagrams for the full model in Figs.~\ref{fig:spher}
and~\ref{fig:F4}.  Equation \eqref{eq:avr} gives $r=1/\sqrt{5}$ for
the full model in the spherical limit.  According to the simple model,
the phase transition at this value of $r$ becomes discontinuous at
$\abs{x} = \abs{\delta\mu/(k_B T_c)} \approx 2.460$, which is in very
good agreement with the location of the tricritcal point for the full
model at $\abs{x} \approx 2.594$ (red dot in Fig.~\ref{fig:spher})
where the effective band splitting is given by
$\delta = \beta/(\alpha + 5\beta/4)$. The simple model also explains
the asymmetric effect of the cubic anisotropy seen in
Fig.~\ref{fig:F4}: For $\abs{\gamma}>\abs{\beta}$, the intraband
pairing potential is enhanced, which, in turn, increases the value of
$r$ and, thus, suppresses the first-order transition. Conversely,
$\abs{\gamma}<\abs{\beta}$ reduces the intraband pairing potential and,
thus, $r$ and favors the first-order transition.

The simple model and our full results agree in showing that a
second-order transition is recovered at sufficiently large values of
the band splitting $\delta$. The reappearance of the second-order
transition in the full model, however, does not occur with a
tricritical point, but rather with a discontinuous jump in the minimum
of the free energy from a large value of the gap magnitude
$(\Delta_{xz}^2 + \Delta_{yz}^2)^{1/2}$ (large-gap phase) to a minimum
at a small value of the gap magnitude (small-gap phase).  Properly
capturing this behavior in the simple model would require extending
the GL expansion in Eq.~\eqref{eq:GLexpand} to, at least, eighth order
in~$\eta$.

\section{Properties of the time-reversal-symmetry-breaking state
\label{sec:trsbstate}}

We now investigate features of the TRSB $C_4$ state. We choose the
parameter set labeled by (c) in Figs.~\ref{fig:spher}
and~\ref{fig:insets}.  In this case, we can set
$\Delta_{xz}=\Delta_{yz}=\Delta_0$, and so, the pairing potential is
$\Delta = \Delta_0\, (\eta_{yz} + \ii \eta_{xz})$.

\begin{figure}[tb]
  \includegraphics{fig09.pdf}
  \caption{BFSs for the parameters labeled by (c) in
    Figs.~\ref{fig:spher} and~\ref{fig:insets} (heavy black lines).
    The colored lines denote the normal-state Fermi surfaces of the
    $+$ and $-$ bands. $k_\perp$ is the radial component of the
    momentum on the $k_x k_y$-plane.}
  \label{fig:fs}
\end{figure}

\subsection{Bogoliubov Fermi surfaces}

First, we map out the BFSs by searching for vanishing energy
eigenvalues.  Thanks to rotational symmetry around the $z$ axis and
inversion symmetry, we can restrict ourselves to the first octant. The
resulting nodal surfaces are shown in Fig.~\ref{fig:fs}.  In the TRSB
$C_4$ state, the magnitude of the pseudomagnetic field in
Eq.~\eqref{eq:pseudomagn} is
\begin{equation}
  \abs{\bm{h}(\bm{k})} = \frac{4 \abs{\Delta_0}^2}{(\epsilon_{\bm{k},+}
    - \epsilon_{\bm{k},-})^2}\,
  \beta\, \sqrt{\abs{\bm{k}}^4 - 3 (k_x^2 + k_y^2) k_z^2} .
\end{equation}
The size of the BFSs scales with the magnitude of the pseudomagnetic
field. Since this field is inversely proportional to the band
splitting squared, which grows as $\abs{\bm{k}}^2$, the inner BFS is
larger than the outer one, see the inset of Fig.~\ref{fig:fs}.  The
pseudomagnetic field has the largest magnitude close to the boundary with
the TRSB $C_2$ state since this corresponds to the smallest band
splitting for which the TRSB $C_4$ state is stable.  Here, the BFSs
have the largest volume and are, therefore, clearly distinguishable from
line and point nodes.

The existence of BFSs leads to a nonzero DOS at zero energy, which is
not expected for clean superconductors.  We compute the DOS
numerically from the mean-field dispersion and analytically using
the low-energy dispersion from Eq.~\eqref{eq:Eeff}.  In the absence of
cubic anisotropy, close to the Fermi surface, $\bm{h}(\bm{k})$,
$\gamma_{\bm{k}}$, and $\psi_{\bm{k},\pm}$ only depend upon the polar
angle $\theta$.  The DOS in the $\pm$ band is, thus,
\begin{equation} \label{eq:DOS}
  \rho_\pm(E) = \mathcal{N}_{0,\pm}
    \sum_{a,b} \int_0^{\pi}
  \frac{\bigl| E - a \abs{\bm{h}_\pm(\theta)} \bigr|\, \sin\theta \diff\theta}
    {\sqrt{(E - a \abs{\bm{h}_\pm(\theta)})^2 - \abs{\psi_\pm(\theta)}^2}} ,
\end{equation}
where we have assumed the normal-state DOS $\mathcal{N}_{0,\pm}$ to be
constant in the range of the superconducting gap.  Evaluating
Eq.~\eqref{eq:DOS}, we find excellent agreement with the numerical
results, as shown in Fig.~\ref{fig:dos}. In particular, we clearly see
a large residual DOS at zero energy in the superconducting gap of up
to 20\% of the normal-state DOS. The flat DOS at zero energy results
from the lifting of the pseudospin degeneracy by the pseudomagnetic
field $\bm{h}$. This shifts the DOS for each pseudospin species,
leading to the scaling
$\rho(E) \propto (\abs{E + \abs{\bm{h}}} + \abs{E - \abs{\bm{h}}})/2$ instead
of $\rho(E) \propto \abs{E}$, as would be the case for line nodes.
This gives a constant DOS for $-\abs{\bm{h}} < E < \abs{\bm{h}}$, as
previously reported in Ref.~\cite{Lapp2019}.  The effect of the
pseudomagnetic field is also seen in the splitting of the coherence
peaks: In the absence of the pseudomagnetic field, we expect a single
coherence peak at $\abs{E}=\Delta_0$. Upon adding the pseudomagnetic
field, it is split into four coherence peaks at
$\Delta_0 + \abs{\bm{h}_{\pm}(\theta = \pi/4)}$ and
$\Delta_0 - \abs{\bm{h}_{\pm}(\theta = \pi/4)}$, where $\theta = \pi/4$ is
the angle of maximum gap. Since the pseudomagnetic field has different
magnitude at the two Fermi surfaces, these two peaks are, in turn,
weakly split.

A residual DOS in an unconventional superconductor can also arise due
to the presence of impurities~\cite{Mineev1999}.  We can estimate the
required concentration of impurities to achieve a zero-energy residual
DOS as large as 20\% of the normal-state DOS within the
self-consistent Born approximation.  Using the exact results for the
polar phase of $^3$He, which also has an equatorial line node, we
estimate that the required concentration of impurities would
approximately result in a 40\% suppression of $T_c$ compared to the
clean limit.  It should be possible to rule out the effect of
impurities by considering the residual DOS as a function of $T_c$ for
different samples.  We also emphasize that the splitting of the
coherence peaks seen in Fig.~\ref{fig:dos} cannot be explained by
impurity effects.

\begin{figure}[tb]
  \includegraphics{fig10.pdf}
  \caption{DOS in the superconducting state for the parameters labeled
    by (c) in Figs.~\ref{fig:spher} and~\ref{fig:insets}, based on a
    full two-band calculation (black curve) and on a low-energy
    single-band approximation (green curve). The results of the two
    approaches agree very well.  The residual DOS at zero energy is as
    large as 20\% of the normal-state DOS at the Fermi energy
    $\mathcal{N}_0=\sqrt{\mu}/2(\alpha + 5\beta/4)^{3/2}$.}
  \label{fig:dos}
\end{figure}

\subsection{Induced magnetic order parameter}

\begin{figure}[tb]
  \includegraphics{fig11.pdf}
  \caption{(a) Gap parameter and (b) induced magnetic order parameter
    as functions of temperature in the TRSB phase and with the SOC
    labeled by (c) in Figs.~\ref{fig:spher} and~\ref{fig:insets}.}
  \label{fig:SI}
\end{figure}

As pointed out in Ref.~\cite{Brydon2018b}, the pseudomagnetic field can be
interpreted as manifesting a subdominant secondary magnetic order
parameter, which is induced by the superconductivity. This subdominant
order is related to the time-reversal odd part of the gap product,
\begin{equation}
  \Delta \Delta^\dagger - U_T\Delta^\ast \Delta^TU_T^\dagger =
  \frac{4}{3}\, \Delta_0^2 (7 J_z - 4 J_z^3)
  \mathrel{\equiv} 2\Delta_0^2 \mathcal{J}_z .
\label{eq:subdominant}
\end{equation}
In Fig.~\ref{fig:SI}, we show the expectation value of $\mathcal{J}_z$
together with the superconducting gap as functions of temperature. The
superconductivity and magnetism appear together but their temperature
dependence close to the critical temperature is notably different:
whereas the gap magnitude scales as $\Delta_0 \sim \abs{T-T_c}^{1/2}$, the
expectation value of $\mathcal{J}_z$ scales as
$\<\mathcal{J}_z> \sim \abs{T-T_c}$. This linear temperature dependence
close to $T_c$ reflects its relation to the gap product in
Eq.~\eqref{eq:subdominant}.

The finite expectation value of $\mathcal{J}_z$ generically leads to a
finite pseudomagnetic field in Eq.~\eqref{eq:pseudomagn} and, thus, to a
momentum-dependent spin polarization. To understand the interplay
between magnetism and superconductivity, we include a magnetic order
parameter $m_z$ in the channel that couples to superconductivity in
the GL expansion.  To that end, following~\cite{Dahl2007}, we redefine:
\begin{equation}
\label{eq:Sigma_new}
  \Sigma =
  \begin{pmatrix}
    \mathcal{M}_z & \Delta \\
    \Delta^\dagger & -\mathcal{M}_z^T \\
  \end{pmatrix} ,
\end{equation}
in Eq.\ \eqref{eq:G0-Sigma}, where
$\mathcal{M}_z = m_z \mathcal{J}_z$.  The lowest-order coupling
between the superconducting and magnetic order parameters occurs at
third order and has the form
\begin{equation}
  \ii F_3m_z (\Delta_{xz}\Delta_{yz}^* - \Delta_{xz}^* \Delta_{yz} ),
\end{equation}
which clearly indicates that the TRSB superconducting state induces
the magnetism.  The lengthy expression for the coefficient $F_3$ is
presented in Appendix~\ref{sec:F3}. In particular, we must introduce a
cutoff $\Lambda$ of the attractive pairing interaction to account for
particle-hole asymmetry in the normal state.  In the limit where the
band splitting and cutoff are much larger than $k_BT_c$ (i.e., the
conditions under which the TRSB state is stable), the coefficient
simplifies to
\begin{equation}
  \label{eq:F3}
  F_3 = \frac{\mathcal{N}_0}{\mu}\,
  \frac{48}{5}
  \biggl[1-  \frac{\ln\frac{2\Lambda e^\gamma}{\pi k_BT_c}}
    {3(1-\tilde{\beta}^2)}
  -\frac{1}{4} \ln
    \biggl(1 + \frac{\Lambda^2}{\tilde{\beta}^2\mu^2}\biggr)\biggr]
    ,
\end{equation}
where $\mathcal{N}_0 = \sqrt{\mu}/2\, (\alpha + 5\beta/4)^{3/2}$ is
the normal-state DOS at the Fermi energy,
$\tilde\beta = \beta/(\alpha + 5\beta/4)$, $\tilde{\beta}\mu$ is the
band splitting, and $\gamma$ is the Euler-Mascheroni constant. To
understand this result, we note that the magnetic order paramater
$m_z$ couples to $\mathcal{J}_z$, which is not diagonal in the band
basis but has both interband and intraband components. The intraband
component directly couples to the pseudomagnetic field generated by
the interband pairing potentials and gives a cutoff-independent
contribution to $F_3$. On the other hand, the interband components of
the magnetic order couple to both the intraband and the interband pairing
potentials and give the cutoff-dependent contribution, see Appendix
\ref{sec:F3} for details.  These two contributions have opposite signs
and the contribution from the interband component is likely dominant
when $\Lambda \gg k_BT_c$.

It is interesting to compare our results to the more familiar case of
coupling between ferromagnetic and superconducting order parameters in
a single-band TRSB superconductor~\cite{Mineev1999}. A similar GL
expansion of the free energy in that case also gives a third-order
coupling term with coefficient proportional to $\mathcal{N}_0/\mu$,
which implies that the magnetization in the superconducting state is
on the order of $\Delta_0^2/\mu^2$ and is, hence, expected to be weak.
This property is thought to be generic for TRSB superconductors
\cite{Cross1975,Mermin1980}. In the present case, it can be understood
as being due to the fact that the $j = 1/2$ and $j = -3/2$
quasiparticles do not participate in the pairing.  The spin of these
unpaired quasiparticles then compensates the polarization of the
Cooper pairs as is the case for a spin-$1/2$ superconductor where
only the up spin is paired and the unpaired down spin compensates the
polarization~\cite{Mineev1999}. The presence of a BFS, therefore, does
not imply a strong magnetization of the superconductor.

\section{Summary and conclusions}
\label{sec:conclusion}

In this paper, we have used BCS mean-field theory to study the
evolution of the quintet superconducting state in the paradigmatic
Luttinger-Kohn model as a function of the SOC strength. We find a rich
phase diagram in the SOC-temperature plane.  For weak SOC, a
time-reversal-symmetric superconducting state is realized. Upon
increasing the SOC strength, the transition into the superconducting
state becomes first order. The origin of the first-order transition is
the competition between inter- and intraband pairings, which is
controlled by the cubic anisotropy of the SOC; for sufficiently
anisotropic SOC, the first-order transition can be completely
suppressed.  Upon further increasing the SOC strength, first a
second-order transition is recovered, and finally a TRSB pairing state
is stabilized.  At low temperatures, the TRSB state displays reentrant
behavior as well as a first-order transition into the TRS-preserving
state.

The TRSB state exhibits BFSs and a residual DOS at the Fermi energy,
which can be as large as 20\% of the normal-state DOS.  The TRSB
pairing state induces a subdominant magnetic order parameter, which we
find to be small even if the residual DOS is sizable, consistent with
the general result that TRSB superconductors have weak intrinsic
magnetization.

Our analysis establishes that a pairing state with BFSs can be
thermodynamically stable, even when the residual DOS at the Fermi
energy due to the BFSs is a sizable fraction of the normal-state
DOS. This result is encouraging for experimental searches for BFSs as
it shows that the residual DOS due to the BFSs can be of detectable
magnitude. Since the size of the BFSs is controlled by the ratio of
the interband pairing potential to the band splitting, materials where
this ratio is as large as possible are the best candidates. This
suggests that heavy-fermion superconductors are promising. It is,
therefore, intriguing that a residual DOS has been observed in
URu$_2$Si$_2$~\cite{Kasahara2007} and UTe$_2$~\cite{Ran2019,
  Metz2019}.

\begin{acknowledgments}
  The authors acknowledge useful discussions with D.~F.  Agterberg and
  M. Vojta. H.\,M. and P.\,M.\,R.\,B. were supported by the Marsden
  Fund Council from Government funding, managed by Royal Society Te
  Ap\=arangi. H.\,M. acknowledges the hospitality of the Technische
  Universit\"at Dresden. C.\,T.  acknowledges the hospitality of the
  University of Otago and financial support by the Deutsche
  Forschungsgemeinschaft through the Collaborative Research Center
  SFB~1143, Project A04, the Research Training Group GRK~1621, and the
  Cluster of Excellence on Complexity and Topology in Quantum Matter
  ct.qmat (EXC~2147).
\end{acknowledgments}