Low energy inelastic response in the superconducting phases of $PrOs_4Sb_{12}$

TL;DR

This paper investigates the low energy inelastic response in the superconducting phases of PrOs4Sb12, providing a microscopic model, Ginzburg-Landau theory, and Raman spectroscopy techniques to identify gap structures and symmetry-breaking phenomena.

Contribution

It introduces a microscopic derivation of the Ginzburg-Landau expansion and proposes Raman spectroscopy methods to distinguish gap symmetries in PrOs4Sb12's superconducting phases.

Findings

Identification of low energy spectral weight differences in A and B phases.

Proposal of Raman spectroscopy geometries to determine gap symmetries.

Potential to locate nodes on the Fermi surface based on spectral features.

Abstract

Recent AC susceptibility and polar Kerr effect measurements in the skutterudite superconductor PrOs4Sb12 (POS) uncovered the nature of the superconducting double transition from a high temperature, high field, time reversal symmetric phase (or the A phase) to a low temperature, low field, time reversal symmetry broken phase (or the B phase). Starting from a microscopic model, we derive a Ginzburg-Landau expansion relevant to POS that describes this entrance into the time reversal symmetry broken phase along the temperature axis. We also provide a study of the low energy inelastic (Raman) response in both the A and B phases of POS, and seek additional signatures which could help reveal the exact form of the gap functions previously proposed in these phases. By appropriately manipulating the incoming and scattered light geometries, along with additional subtraction procedures and suitable assumptions, we show that one can access the various irreducible representations contained in the point group describing POS. We demonstrate how to use this technique on example order parameters proposed in POS. Depending on whether there exist nodes along the c-axis, we find additional low energy spectral weight within the superconducting gap in the Eg geometry, a feature that could pin point the location of nodes on the Fermi surface.

Tables

Table 1. Table 1: Character table for the T h subscript 𝑇 ℎ T_{h} point group. In the basis function column we have defined r s 2 = x 2 + y 2 superscript subscript 𝑟 𝑠 2 superscript 𝑥 2 superscript 𝑦 2 r_{s}^{2}=x^{2}+y^{2} and r a 2 = x 2 − y 2 superscript subscript 𝑟 𝑎 2 superscript 𝑥 2 superscript 𝑦 2 r_{a}^{2}=x^{2}-y^{2} .

$T_h$	E	$3C_2$	$8C_3$	$i$	$3i{C}_2$	$8i{C}_3$	Basis
$A_g$	1	1	1	1	1	1	$r_s^2+z^2$
$E_g$	2	2	-1	2	2	-1	($r_a^2$,$2z^2-r_s^2$)
$T_g$	3	-1	0	3	-1	0	$(xz,yz,xy)$
$A_u$	1	1	1	-1	-1	-1	-
$E_u$	2	2	-1	-2	-2	1	-
$T_u$	3	-1	0	-3	1	0	$(x,y,z)$

Equations

$${ϵ}_1(\vec{k})$$

$${ϵ}_2(\vec{k})$$

$${ϵ}_3(\vec{k})$$

$$A_g(\vec{k})$$

$${\vec{E}}_g(\vec{k})$$

$${\vec{T}}_u(\vec{k})$$

$$S(\bar{c},c)$$

$${\mathcal{F}}_s$$

$$D_s$$

$$D_d$$

$$\frac{{\partial }^2\sigma }{\partial \omega \partial \Omega }$$

$$S_{\gamma \gamma }(\vec{q},\omega )$$

$${\chi }_{\gamma \gamma }(\omega )={\int }_0^{\beta }d\tau e^{-i{\omega }_m\tau }⟨T_{\tau }{\tilde{\rho }}_{\gamma }(\tau ),{\tilde{\rho }}_{\gamma }(0)⟩{|}_{i{\omega }_m\to \omega +i\delta }$$

$$\widetilde{{\rho }_{\gamma }}=\sum _{\vec{k},\sigma }\sum _{n,m}{\gamma }_{n,m}(\vec{k})c_{n,\sigma }^{†}(\vec{k})c_{m,\sigma }(\vec{k}),$$

$${\gamma }_n(\vec{k})$$

$${\gamma }_n(\vec{k})=\sum _{\mu }{\gamma }_n^{\mu }{\mathcal{I}}_n^{\mu }(\vec{k})$$

$${\chi }_{\gamma \gamma }^{(n)}(i\omega )$$

$${\lambda }_n(\vec{k},i\omega )=\frac{{\Delta }_n(\vec{k}{)}^2}{E_n(\vec{k}{)}^2}\left(\frac{1}{2E_n(\vec{k})+i\omega }+\frac{1}{2E_n(\vec{k})-i\omega }\right),$$

$${ϵ}_{E_g}(\vec{k})$$

$${\vec{E}}_g(\vec{k})$$

$${\gamma }_n(\vec{k})$$

$$R_{\Gamma }^{\mu \nu }$$

$$R_{E_g}^{\mu \nu }=\left(\begin{array}{ccc}3k_x^2-k_z^2& 0& -2k_x{k}_z\\ 0& 3k_y^2-k_z^2& -2k_y{k}_z\\ -2k_x{k}_z& -2k_y{k}_z& -k_x^2-k_y^2+6k_z^2\end{array}\right).$$

$${\gamma }_n(\vec{k}{)}_{\pm }$$

$${\gamma }_n(\vec{k}{)}_{\pm }=\gamma (\vec{k}{)}_{A_g}\pm \gamma (\vec{k}{)}_{T_g(1)},$$

$${\gamma }_n(\vec{k}{)}_{\pm }=\gamma (\vec{k}{)}_{E_g(2)}\pm i\gamma (\vec{k}{)}_{T_g(1)},$$

$${ϵ}_{A_g}(\vec{k})=c_x+c_y+c_z+c_x{c}_y+c_y{c}_z+c_x{c}_z-\mu ,$$

$$R_{A_g}^{\mu \nu }=\left(\begin{array}{ccc}-C_{x;yz}& s_x{s}_y& s_x{s}_z\\ s_x{s}_y& -C_{y;zx}& s_y{s}_z\\ s_x{s}_z& s_y{s}_z& -C_{z;xy},\end{array}\right).$$

$$S(\bar{c},c)$$

$$S(\bar{\psi },\psi )$$

$$\mathcal{Z}=\int \mathcal{D}(\psi ,\bar{\psi })e^{-S(\psi ,\bar{\psi })}.$$

$$\pi g{e}^{gab}=\int d\varphi d\bar{\varphi }{e}^{a\varphi +b\bar{\varphi }}{e}^{\frac{-|\varphi {|}^2}{g}}$$

$$\mathcal{Z}$$

Full text

I Introduction

Ever since its original discovery early last decade Bauer et al. (2002); Maple et al. (2002), the heavy fermion skutterudite superconductor $\mathrm{PrOs_{4}Sb_{12}}$ has met every measure used to quantify unconventionality in superconductors Maple et al. (2006); Aoki et al. (2007). Two distinct specific heat jumps have been observed Vollmer et al. (2003) at temperatures close to $T_{c}\sim 1.85K$ and $T_{c}\sim 1.7K$ indicating two different transitions to the superconducting state in the absence of a magnetic field. This conclusion has been backed by additional thermal expansion Oeschler et al. (2003) and transport Izawa et al. (2003) measurements which report identical jump-like features supporting a double transition to the superconducting state. The power-law behavior of the Nuclear Magnetic Resonance (NMR) spin-lattice relaxation rate Tou et al. (2005, 2011), Nuclear Quadrupole Resonance (NQR) Katayama et al. (2007), penetration depth Chia et al. (2003) and specific heat Vollmer et al. (2003); Bauer et al. (2002); Maple et al. (2002) seem to suggest a point-nodal superconducting gap-like behavior, although other thermal conductivity, muon spin resonance ($\mu$SR) and NQR studies seem to support a fully gapped Fermi surface Seyfarth et al. (2006); MacLaughlin et al. (2002); Kotegawa et al. (2003). The presence of a non-zero Kerr angle Aoki et al. (2003); Levenson-Falk et al. (2016) points toward a chiral time reversal symmetry broken (TRSB) superconducting ground state, with a more recent measurement Levenson-Falk et al. (2016) reporting that the TRSB state is specific to the low temperature, low magnetic field superconducting phase (B phase) and absent in the high temperature, high magnetic field phase (A phase). Additionally, Knight shift measurements Higemoto et al. (2007) have shown signs of spin triplet pairing and, along with TRSB, has motivated a recent proposal Kozii et al. (2016) suggesting that POS could prove to be a promising candidate for a chiral superconductor hosting three-dimensional Majorana fermions. In the normal state of POS, low temperature specific heat measurements Aoki et al. (2002) in the presence of a magnetic field have identified a broad magnetic field induced ordered phase characterized by an enhancement of the 4$f$ magnetic moment.

A fundamental question concerning the character of the superconducting ground state is the pairing mechanism and symmetry of the order parameter defined on the Fermi surface. Several theoretical efforts have already been put forth to classify and probe the pairing mechanism and symmetry in POS Sergienko and Curnoe (2004); Maki et al. (2003, 2004); Goryo (2003); Miyake et al. (2003); Ichioka et al. (2003). The structure of the gap functions in the spin singlet pairing channel have been considered Goryo (2003) using a Ginzburg-Landau (GL) phenomenological approach in both the A and B phases. In the absence of any spin orbit coupling, it was concluded that in the A phase a strongly anisotropic $s$-wave is favored; using the basis functions for the irreducible representations of the $T_{h}$ point group (see the Character table in Table 1), a single order parameter component of the forms $\Delta(\vec{k})=s+k_{x}^{4}+k_{y}^{4}+k_{z}^{4}$ or $\Delta(\vec{k})=s+k_{x}^{2}k_{y}^{2}+k_{y}^{2}k_{z}^{2}+k_{z}^{2}k_{x}^{2}$ with an $A_{g}$ symmetry were argued to be competent ($s$ is a scalar number). In the B phase, a TRSB $s+id$ state with a gap structure of the form $\Delta(\vec{k})=(s+k_{x}^{2}k_{y}^{2}+k_{y}^{2}k_{z}^{2}+k_{z}^{2}k_{x}^{2})$ + $i(k_{x}^{2}-k_{y}^{2})$, with combinations of the one dimensional $A_{g}$ and two dimensional $E_{g}$ representations was argued to be the ground state. In the presence of spin orbit coupling, an $f$-wave pairing state with point nodes along all the three axes in momentum space was suggested in Ichioka et al. (2003), an novel $s+g-$ wave pairing was found in ref Maki et al. (2003), and, using a microscopic approach with quadrupole fluctuation mediated pairing, Miyake and coworkers Miyake et al. (2003) found a chiral $p_{x}+ip_{y}$ pairing order parameter. Several other order parameters were proposed Sergienko and Curnoe (2004) based on a rigorous symmetry analysis in which the authors, following the work of Volovik and Gor’kov Volovik and Gor kov (1985), minimized a Landau free energy functional that was invariant under the point group, gauge and time reversal symmetries, and studied the ensuing gap nodal structures.

Electronic Raman scattering in superconductors Klein and Dierker (1984); Devereaux and Hackl (2007); Devereaux and Einzel (1995); Devereaux et al. (1994, 1996) has proved to be an indispensable tool for determining the structure and symmetry of the gap in unconventional, anisotropic superconductors. Light is incident with a certain initial electric polarization and energy ($\hat{e}^{i},\omega_{I}$) and scattered with a final polarization and energy ($\hat{e}^{s},\omega_{S}$). The energy difference corresponds to the energy needed to break a Cooper pair or to excite a low energy in-gap mode. The Raman scattering vertex, in general, can have contributions from all the irreducible representations (IRRP), $\Gamma_{i}$, of the point group of the crystal; however, a proper choice of the polarization geometries of the incoming and outgoing light beams can single out contributions from a single IRRP to the total scattering cross section. As a result of this selective Brillouin zone averaging, the quasiparticle energy gaps are probed only in certain directions in the Brillouin zone. Thus this technique is sensitive to the position of gap nodes, and therefore provides indirect information about the symmetry of the gap function and its anisotropies. If there are several candidate pairing symmetries with different location of nodes on the Fermi surface, this method can (at least) help narrow down the various possible candidates.

In this work, starting from a microscopic model, we derive a GL expansion relevant to POS that shows its entry into the TRSB phase at low temperatures. We also estimate the GL coefficients using their relationships to the various microscopic parameters. As an experimental means to probe the proposed pairing symmetries, we examine the low energy inelastic (Raman) response in both the A and B phases of POS; in particular, we seek to qualitatively understand its dispersive features for low energy and zero momentum transfer. By appropriate manipulation of the incoming and scattered light geometries, as well as other subtraction procedures routinely applied Devereaux and Einzel (1995) in this method, we demonstrate that one can probe the various irreducible representations contained in the $T_{h}$ point group describing POS. For the purposes of illustration, we stick to the spin singlet, even parity forms of the gap functions proposed by Goryo Goryo (2003) in the A and B phases, as was discussed in the preceeding paragraphs. Depending on the existence of nodes along the $c-$ axis, we find an enhancement (nodal) or suppression (gapped) of the low-energy spectral weight in the $E_{g}(1)$ and $E_{g}(2)$ geometries (here (1), (2) etc. denote the individual elements of the respective multi-component IRRP-). This constitutes a defining property in electronic inelastic light scattering that could help pin-point the exact location of nodes on the fermi surface and, thereby, narrow down the candidate pairing symmetries proposed for POS.

II Lattice and Fermi surface

The lattice structure of POS is shown in Fig 1(left). The relevant space group of interest is the $Im\bar{3}$ (No. 204) with a tetrahedral $T_{h}$ point group. The lattice structure consists of Pr atoms (shown in magenta) forming a body centered cubic lattice and the Sb atoms (shown in gray balls) forming icosahedral cages (yellow structures) encapsulating the Pr atoms. The Os atoms (not shown for purposes of clarity) form a cube and are intercalated within the body centered cubic structure. The icosahedral cage is shown in more detail in Fig 1 (right). Each cage consists of three rectangular planes (shown in red, green and blue) of Sb atoms lying perpendicular to each other, along the three coordinate planes, with a common intersection point at the origin. The symmetry group describing of the Pr atoms, by themselves, is the cubic $O_{h}$ group, while that of the Sb atoms is the tetrahedral $T_{h}$ group. Since the cubic $O_{h}$ is a larger symmetry compared to the tetrahedral $T_{h}$ (all symmetry elements in $T_{h}$ are also present in $O_{h}$), the symmetry of the entire lattice is determined by the smaller $T_{h}$ point group. One should, therefore, use the irreducible representations of $T_{h}$ (for example, see refs Dresselhaus et al. (2007); Cardona and Peter (2005)) to construct invariants to describe POS (this is made possible because the symmetry axes belonging to the Sb icosahedral cages coincide with that of the Pr body centered cubic lattice ).

Several quantum oscillation measurements Sugawara et al. (2002, 2009, 2008) have been performed to map out the Fermi surface structure of POS. The Brillouin zone is a rhombicdodecahedron (just as in the case of a body centered cubic lattice), with its origin at the $\Gamma$ point and the center of the rhombic face labelled as the $N$ point (for more details about the other high symmetry points and axes, the reader can refer to Bradley and Cracknell (2010)). All quantum oscillation measurements observe two hole pockets centered around the $\Gamma$ point $-$ an inner spherical pocket and an outer rounded cubic pocket. There is another larger multiply connected pocket, centered again around the $\Gamma$ point, and touching the face of the boundary rhombus at the $N$ point. The Brillouin zone and the Fermi surface of POS is shown in fig 2.

Currently, we do not know of any exhaustive study attempting to write down a simple tight-binding model that captures all the essential kinematics such as the band structure, Fermi surface and orbital character of POS. However, to study the electronic Raman response of a material, we only need the underlying symmetry properties of POS. Hence, for the purposes of this paper, we will only be interested in the symmetry characters of the individual bands and not other microscopic details. To achieve this, and at the same time maintain analytical and numerical tractability of the problem, we will work in the band basis and ignore any interband effects. This is a reasonable assumption as the low energy interband response can be neglected for widely spaced bands. To construct these dispersions, we will rely on the basis functions appearing in the character table presented in Table 1. For the three Fermi surfaces shown in Fig 2, we simply construct invariants (i.e functions that transform as a singlet $A_{g}$) with respect to symmetry operations of the $T_{h}$ point group. We have chosen the invariants such that the Fermi surfaces match the experimentally measured contours in refs Sugawara et al. (2002, 2009, 2008) through quantum oscillations. These invariants are given as

[TABLE]

where $(\epsilon_{1},\epsilon_{2},\epsilon_{3})=(-6.88,-6.7,3)eV$, $\mu=-1eV$, and the basis functions are given by

[TABLE]

Here, we have defined $c_{x}=\cos(k_{x})$, $s_{\frac{x}{2}}=\sin(k_{x}/2)$ and so on. It is easy to check that Eqs (2) transform according to their respective irreducible representations (albeit in a lattice version) as specified in the character table in Table 1. It should also be noted that in the last equation for $\vec{T}_{u}(\vec{k})$, we have used higher order basis functions that are not present in Table 1. This state of affairs arises because writing the dot product invariants with the lowest order basis functions of the $T_{u}$ representation yields a function equivalent to $\epsilon_{1}(\vec{k})$. The three dispersions in Eqs. 1 provide the Fermi surface contours in Fig 2.

III Time reversal symmetry breaking

In the following analysis, we derive a generalized GL expansion to understand time reversal symmetry breaking in POS. On the theoretical sideShiina (2004); Shiina and Aoki (2004); Shiina et al. (2004), the effect of magnetic and quadrupole moments on the phase diagram was outlined. Motivated by the evidence of such magnetic moments Kohgi et al. (2003); Kuwahara et al. (2004); Maple et al. (2002); Bauer et al. (2002); Tayama et al. (2003); Aoki et al. (2002), magnetic exchange correlations, quadrupole moments Kohgi et al. (2003); Tayama et al. (2003); Aoki et al. (2002); Rotundu et al. (2004) and quadrupole exchange correlations, and taking into account the crystal symmetries appropriate for quartic interactions generated by such correlations, we start with the following form of the action written in momentum space

[TABLE]

Here $\bar{c}_{\alpha\sigma}(\textbf{k})$ and $c_{\beta\sigma}(\textbf{k})$ are the creation and annihilation Grassmann numbers that follow Grassmannian algebra for electrons with Bloch momentum $\vec{k}$ and Matsubara frequency $k_{n}$ (collectively denoted by k), spin $\sigma$ and orbital indices $\alpha,\beta$, $\epsilon_{\alpha\beta}(\vec{k})$ are the orbital matrix elements, $g$ is the interaction strength which is taken to be a constant, $\beta$ is the inverse temperature, and $D(\vec{k})\equiv\phi_{s}(\vec{k})+\phi_{d}(\vec{k})$ is the total form factor which is a mixture of the $s-$ wave and $d-$wave form factors. Here we have defined $\phi_{s}(\vec{k})=\Delta_{s}\Phi_{s}(\vec{k})$ and $\phi_{d}(\vec{k})=\Delta_{d}\Phi_{d}(\vec{k})$, where $\Delta_{s}$ and $\Delta_{d}$ are complex scalar amplitudes, along with $\Phi_{s}(\vec{k})=(k_{x}^{2}k_{y}^{2}+k_{y}^{2}k_{z}^{2}+k_{z}^{2}k_{x}^{2})$ and $\Phi_{d}(\vec{k})=(k_{x}^{2}-k_{y}^{2})$. The spin structure in the interaction term is chosen since we anticipate a BCS-like singlet superconducting order parameter. At this point, from an experimental perspective, it is not fully clear if the gaps on the different orbitals need to be distinct or not; hence from now on, for simplicity, we will choose the total form factor $D(\vec{k})$ to be independent of the orbital index $\alpha$. Such an assumption will allow us to work in a basis where the kinetic energy and the order parameter matrices appearing in the Hamiltonian are both diagonal because they commute and hence, will simplify the problem to a collection of independent bands with the same superconducting gap on each band. Note that this is not always true in generic multiband superconductors (for example, the Iron- based superconductors) where the gap and kinetic energy matrices do not commute. Changing to the band basis with energy eigenvalues $\epsilon_{\alpha}(\vec{k})$ by using the unitary transformation, $c_{\alpha\sigma}(\textbf{k})=U_{\alpha\beta}(\vec{k})\psi_{\beta\sigma}(\textbf{k})$, performing a Hubbard-Stratonovich transformation to decouple the quartic terms in the action, and integrating out the fermionic degrees of freedom results in a free energy in the superconducting state of the form (see Appendix and Ref. Wang and Fu (2017))

[TABLE]

with the detailed form of the coefficients given in the Appendix. For the case of the $s-$ and $d-$wave symmetries chosen in our paper for POS, $\alpha^{\prime}$ and $g_{\pm}$ are small compared to the rest of the coefficients due to negligible overlap between odd powers of the $s-$ and $d-$ wave form factors, and hence can be neglected in the lowest-order approximation. We will later include these terms to see their effects on the phase diagram. If we define the relative phase between $\Delta_{s}$ and $\Delta_{d}$ to be $\eta$, minimization of the free energy with respect to $\eta$ fixes $\eta=\pi/2$. Minimizing with respect to $|\Delta_{s}|$ and $|\Delta_{d}|$ gives us four sets of solutions of the form, $\left(|\Delta_{s}|,|\Delta_{d}|\right)=(0,0)$, $\left(|\Delta_{s}|,|\Delta_{d}|\right)=(0,\sqrt{\frac{-(\alpha-\delta\alpha)}{2\beta_{d}}})$, $\left(|\Delta_{s}|,|\Delta_{d}|\right)=(\sqrt{\frac{-(\alpha+\delta\alpha)}{2\beta_{s}}},0)$, $\left(|\Delta_{s}|,|\Delta_{d}|\right)=(D_{s},D_{d})$, where

[TABLE]

Here we have defined $\alpha=(\alpha_{s}+\alpha_{d})/2$ and $\delta\alpha=(\alpha_{s}-\alpha_{d})/2$. The last pair of solutions where the relative phase between $\Delta_{s}$ and $\Delta_{d}$ is fixed by the $\beta_{2}$ term at $\pm\pi/2$ gives rise to the TRSB phase. Such a phase has the lowest free energy but is obviously a stable solution only when $(|\Delta_{s}|,|\Delta_{d}|)$ are both real. Hence, the conditions $\frac{\alpha(\beta_{1}-\beta_{d})-\delta\alpha(\beta_{1}+\beta_{d})}{2(\beta_{s}\beta_{d}-\beta_{1}^{2})}>0$ and $\frac{\alpha(\beta_{s}-\beta_{1})-\delta\alpha(\beta_{s}+\beta_{1})}{2(\beta_{1}^{2}-\beta_{s}\beta_{d})}>0$ define the TRSB phase. The $s-$wave and $d-$wave pairings are stable when $-(\alpha+\delta\alpha)$ and $-(\alpha-\delta\alpha)$ are positive respectively. These phases are shown in Fig 3(a-c) for $\beta_{s}=\beta_{d}$ (Fig 3(a)), $\beta_{s}>\beta_{d}$ (Fig 3(b)) and $\beta_{s}<\beta_{d}$ (Fig 3(c)). Given the candidate pairing symmetries we have chosen for POS, viz. $\Phi_{s}(\vec{k})$ and $\Phi_{d}(\vec{k})$, we find that POS belongs to the case where $\beta_{s}<\beta_{d}$.

The presence of linear terms with the coefficients $\alpha^{\prime}$ and $g_{\pm}$ reflects the fact that the $s$-wave and $d$-wave order parameters are not distinguished by symmetry. Then a transition between $s$-wave and $d$-wave states does not require further breaking any symmetries – in fact, the $s$ and $d$-wave components are always mixed even in the time-reversal invariant phase. Thus when $\Delta_{s}\sim\Delta_{d}$, the two order parameters coexist in two competing ways, either preserving or breaking time-reversal symmetry. Since the linear coupling terms are already present at the quadratic level, while the terms responsible for time-reversal symmetry breaking enter at the quartic level, the mixing of the two orders breaks time-reversal only at low-temperatures when order parameters have larger expectation values. Treating the $\alpha^{\prime}$ and $g_{\pm}$ terms perturbatively, we obtain a schematic phase diagram shown in Fig 3(d). A more detailed calculation for obtaining this phase diagram is found in Ref Maiti and Chubukov (2013).

IV Inelastic (Raman) scattering

The scattered light intensity is written in terms of the differential scattering cross section Klein and Dierker (1984); Devereaux and Hackl (2007); Devereaux and Einzel (1995); Devereaux et al. (1994, 1996) as

[TABLE]

where $n_{B}(\omega)$ is the Bose-Einstein distribution function, $\omega_{S}$ and $\omega_{I}$ are the frequencies of the scattered and incident light respectively and $r_{0}=\frac{e^{2}}{mc^{2}}$. The imaginary part of the inelastic response $Im\chi_{\gamma\gamma}(\vec{q},\omega)$ (the subscript $\gamma$ denotes that the fluctuations are weighted by a vertex function) is related to the generalized structure factor $S_{\gamma\gamma}$ through the fluctuation-dissipation theorem. The inelastic response, in the long wavelength limit, is sensitive to effective density fluctuations $\left(\chi_{\gamma\gamma}(\omega)\equiv\chi_{\gamma\gamma}(0,\omega)\right)$ given by

[TABLE]

where the effective (weighted) density is written as

[TABLE]

$n,m$ denote band indices, $\gamma_{n,m}$ is the vertex, and $c_{n,\sigma}^{\dagger}(\vec{k})$ is the electron creation operator for band $n$, momentum $\vec{k}$ and spin $\sigma$. If we are interested only in the low energy response, we can neglect the off-diagonal vertices by substituting $\gamma_{n,m}\rightarrow\delta_{mn}\gamma_{n}$, and thereby ignore interband transitions to write the vertex in the form

[TABLE]

Here we have defined the electron mass $m$, $|n,\vec{k}\rangle$ as the Bloch state with band index $n$ and crystal momentum $\vec{k}$, $\hat{e}^{i,s}$ denotes the polarization directions of the incident and scattered light respectively, $p$ is the momentum operator and $\epsilon_{n}(\vec{k})$ are the Bloch energies. The vertex function on the $n$th band can be broken down into various contributions from the IRRPs of the point group as

[TABLE]

where the index $\mu$ denotes the contributions from the different point group irreducible representations and the functions $\mathscr{I}^{\mu}_{n}(\vec{k})$ are the corresponding basis functions of the $\mu$’th IRRP in the $n^{\rm th}$ band. The Raman response from the $n^{\rm th}$ band can then be written as (in the independent band approximation)

[TABLE]

with $\lambda_{n}$ the Tsuneto function given byDevereaux and Hackl (2007); Devereaux and Einzel (1995); Devereaux et al. (1994, 1996)

[TABLE]

where $\Delta_{n}(\vec{k})$ is the gap function on the $n^{\rm th}$ band and $E_{n}(\vec{k})^{2}$ are the quasiparticle energies and equals $\epsilon_{n}(\vec{k})^{2}+\Delta_{n}(\vec{k})^{2}$. As we assumed earlier, for simplicity, we will take the gap function to be independent of the orbital or band indices and denote it by $\Delta(\vec{k})$. We will use the expression in Eq. (11) to evaluate the response functions in the following sections. We would like to point out that in deriving the expressions in Eq. (11), we have ignored the role of vertex corrections through final state interactions, and long range Coulomb interactions. A recent work Maiti et al. (2017) studied these effects in a generic multi-band system and concluded that vertex corrections remove the $2\Delta_{0}$ singularity and create a broad peak at higher energies. These conclusions should not invalidate our low energy calculations, but must only modify the exact position of the peaks. The authors also find that in the $q\rightarrow 0$ limit, long range Coulomb interactions are irrelavant to the Raman response for all polarization geometries. Additional effects arising from subleading interaction channels due to vertex corrections leading to collective in-gap ($\omega<2\Delta_{0}$) Leggett and Bardasis-Schrieffer modes will not be considered in this study. We have also assumed that the light scattering occurs within a region where superconductivity is uniform. In the presence of domain walls separating different superconducting states, the responses must be averaged over the domains, and isolating the pairing symmetry is not straightforward.

IV.1 Light Polarization and Symmetry:

Square Lattice

As a warm up, we recall the arguments for the case of probing a $d$-wave gap for a square lattice (e.g. the Cuprates)Devereaux and Einzel (1995). We start with a dispersion of the form $\epsilon(\vec{k})=-2t(\cos(k_{x})+\cos(k_{y}))+4t^{\prime}\cos(k_{x})\cos(k_{y})-\mu$ where we choose $(t,t^{\prime},\mu)=(1,0.65,-1)eV$. This gives us a Fermi surface contour shown in the right panel of fig 4. The crucial feature that makes inelastic light scattering of Bloch electrons powerful, is the ability of the experimentalist to control the vertex function $\gamma_{n}(\vec{k})$ by manipulating the incoming and scattered light polarization geometries $\hat{e}^{i}$ and $\hat{e}^{s}$. For a generic multiband system (one can think of a single band system as a multiband system with infinite energy separation between bands), as long as the bands are substantially (compared to the frequency of incoming/scattered light) far apart near the Fermi energy, one can approximate the equation for the vertex function in Eq (IV) as a band curvature Devereaux and Einzel (1995); Devereaux et al. (1996). Note that this approximation breaks down when the bands touch each other at or near the fermi level and could yield qualitatively different results Setty and Hu (2014). Keeping in mind that there are three distinct irreducible representations, $-A_{1g},B_{1g},B_{2g}-$ that have been accessed in the square lattice point group, each IRRP (selected through its appropriate polarization geometry) can selectively average different regions of the Fermi surface, yielding distinct responses in the superconducting state. The basis functions for the square lattice IRRPs are given by $B_{1g}\sim\cos(k_{x})-\cos(k_{y})$, $B_{2g}\sim\sin(k_{x})\sin(k_{y})$ and $A_{1g}\sim\cos(k_{x})+\cos(k_{y})$. The $B_{2g}$ geometry can be accessed, for example, by setting $\hat{e}^{i}=(1,0),\hat{e}^{s}=(0,1)^{T}$ or vice-versa, while the $B_{1g}$ representation can be obtained by rotating the $B_{2g}$ polarization vectors by $\pi/4$. To access the $A_{1g}$ representation, one needs to use the left and right circularly polarized light for the incoming/outgoing polarization vectors. The low energy inelastic response in these geometries for the square lattice is plotted in Fig 4(left) for a $d-$wave gap function of the form $\cos(k_{x})-\cos(k_{y})$. The low energy curves behave very anisotropically for different geometries. For example, in the $A_{1g}$ and $B_{2g}$ cases, they behave linearly for small $\omega$ while for the $B_{1g}$ case, the response is super-linear ($\sim\omega^{3}$) Devereaux and Einzel (1995); Devereaux and Hackl (2007). To understand this behavior one has to look at the regions of the Fermi surface that are being sampled$-$ in the $B_{2g}$ geometry, the momentum average is focussed along the $k_{x}=\pm k_{y}$ directions where the $d-$wave nodes are located, while the $B_{1g}$ geometry samples the anti-nodal regions (along the $k_{x}$ or $k_{y}$ axes) of the Fermi surface. Thus the low energy response in the $B_{1g}$ grows slower ($\sim\omega^{3}$) relative to the $B_{2g}$ case ($\sim\omega$). In the $A_{1g}$ scenario, there is sampling of both the nodal and anti-nodal regions and, thus, the low energy response is dominated by the nodes ($\sim\omega$). In the following paragraphs, we will follow a similar line of argument for the $T_{h}$ point group.

IV.2 Light polarization and symmetry$-T_{h}$ point group

In this section, we will highlight the connections between light polarizations that can be experimentally manipulated, and the symmetry of the Raman vertex in the Brillouin zone for the $T_{h}$ point group. For simplicity, we will illustrate this connection using dispersions obtained by individual IRRPs for the calculation of the Raman vertex. The total vertex for the band structures (used above) can then be obtained by a simple linear combination of the vertices of these representations. In order to keep our equations tractable in this example, we will work in the continuum limit; this can be easily generalized to the lattice case with the inclusion of trigonometric functions in momenta instead of polynomials. Let us begin with a dispersion obtained from the $\vec{E}_{g}$ representation appearing as a term in Eqs. 1. In the continuum limit, this is given by

[TABLE]

with $\mu$ the chemical potential. We can evaluate the Raman vertex for this dispersion using the band effective mass approximation for the $n^{\rm t}$ band asDevereaux et al. (1996)

[TABLE]

where, $\epsilon_{\Gamma}(\vec{k})$ is the band dispersion formed from the irreducible representation $\Gamma$ and $\mu,\nu$ are the coordinate indices. For these dispersions, the Raman tensor $R^{\mu\nu}$ is

[TABLE]

An arbitrary choice of the incoming and scattered polarization vectors yields a combination of $A_{g}$, $\vec{E}_{g}$ and $\vec{T}_{g}$ representations. We can obtain a pure representation by choosing $\hat{e}^{i}$ and $\hat{e}^{s}$ such that only the representation of interest contributes to the scattering cross section. This is not always straightforward, as we will see, and one will need to resort to subtraction procedures. The simplest cases to access are the $E_{g}(1)$ and $\vec{T}_{g}$ representations. Choosing $\hat{e}^{i}=(1,-1,0)^{T}$ and $\hat{e}^{s}=(1,1,0)$, it is easy to verify that $\sum_{\mu\nu}\hat{e}_{\mu}^{i}R^{\mu\nu}_{\Gamma}\hat{e}^{s}_{\nu}=3(k_{x}^{2}-k_{y}^{2})$, which transforms as the first component of the $\vec{E}_{g}$ representation. Similarly, choosing $\hat{e}^{i}=(0,0,1)^{T}$ and $\hat{e}^{s}=(1,0,0)$, or, $\hat{e}^{i}=(0,0,1)^{T}$ and $\hat{e}^{s}=(0,1,0)$, picks up the $T_{g}(1)$ or $T_{g}(2)$ representations respectively (the $T_{g}(3)$ component is zero here). To obtain the $E_{g}(2)$ and $A_{g}$ representations, we need to choose two different sets of incoming and scattered polarizations. As our first choice, ($c_{+}$), we select $\hat{e}^{i}=(\sqrt{\alpha},i\sqrt{\alpha},\sqrt{\beta})^{T}$ and $\hat{e}^{s}=(\sqrt{\alpha},-i\sqrt{\alpha},\sqrt{\beta})$, and for our second choice ($c_{-}$), we select $\hat{e}^{i}=(\sqrt{\alpha},i\sqrt{\alpha},-\sqrt{\beta})^{T}$ and $\hat{e}^{s}=(\sqrt{\alpha},-i\sqrt{\alpha},-\sqrt{\beta})$. Here, $\alpha$ and $\beta$ are constants that must be determined. These two choices $c_{\pm}$ for the polarization give us vertices

[TABLE]

At this point, the constants $\alpha$ and $\beta$ can be related to each other in two useful ways. We consider first the case where the coefficients of $(k_{x}^{2}+k_{y}^{2})$ and $k_{z}^{2}$ are equal, i.e, when $(3\alpha-\beta)=(6\beta-2\alpha)$ or $5\alpha=7\beta$. With this condition, the vertices $\gamma_{n}(\vec{k})_{\pm}$ in Eq (IV.2) are the sum and difference of the $A_{g}$ ($k_{x}^{2}+k_{y}^{2}+k_{z}^{2}$) and $T_{g}(1)$ ($k_{x}k_{z}$) representations alone. Thus, we have

[TABLE]

where $\gamma(\vec{k})_{A_{g}}$ and $\gamma(\vec{k})_{T_{g}(1)}$ are vertices that individually transform as $A_{g}$ and $T_{g}(1)$ representations respectively. Summing the responses from $\gamma_{n}(\vec{k})_{+}$ and $\gamma_{n}(\vec{k})_{-}$ using Eq. ( 11), and noting that we know how to obtain a pure $T_{g}(1)$ contribution, we can calculate the response from a pure $A_{g}$ representation by a mere subtraction of the two resulting responses. For the pure $E_{g}(2)$ response, we equate the coefficients of the $(k_{x}^{2}+k_{y}^{2})$ and $k_{z}^{2}$: $6\beta-2\alpha=-2(3\alpha-\beta)$ or $\alpha=-\beta$. Substituting this condition into Eq. (IV.2), we obtain the vertices as

[TABLE]

where $\gamma(\vec{k})_{E_{g}(2)}$ is the vertex that transforms as the second component of the $E_{g}$ representation. Again using Eq. (11), for either $\gamma_{n}(\vec{k})_{\pm}$, and equipped with knowledge of the $T_{g}(1)$ response, we can isolate the response in the $E_{g}(2)$ channel by subtracting the two resulting responses.

As an example of a dispersion that can be obtained through an $A_{g}$ representation in Table 1, we choose a dispersion of the form

[TABLE]

which is a generalization of the function $A_{g}(\vec{k})$ used in eq 2. The Raman tensor $R^{\mu\nu}$ for this dispersion is given as

[TABLE]

where we have defined $C_{i;jk}=-c_{i}(1+c_{j}+c_{k})$. To obtain the components of a pure $\vec{T}_{g}$ representation, we can follow arguments analogous to what we did in the previous paragraph. For the third component of the $\vec{T}_{g}$ representation, we need to choose $\hat{e}^{i}=(0,1,0)^{T}$ and $\hat{e}^{s}=(1,0,0)$, and similarly for the first and second components. To obtain the response in the $A_{g}$ channel, we choose $c_{+}$ to be $\hat{e}^{i}=(1,i,1)^{T}$ and $\hat{e}^{s}=(1,-i,1)$ and select $c_{-}$ to be $\hat{e}^{i}=(1,i,-1)^{T}$ and $\hat{e}^{s}=(1,-i,-1)$. The resulting vertices corresponding to these two choices are the sum and difference of the $A_{g}$ and $T_{g}(1)$ symmetries; the pure $A_{g}$ response can then be extracted by using eq 11 for these vertices, along with our knowledge of the pure $T_{g}(1)$ response. One can repeat this analysis for a dispersion $\epsilon_{T_{u}}(\vec{k})$ (of the form appearing in eq 1) obtained from the basis functions of the $\vec{T}_{u}$ representation in Table 1. Fig 5 shows regions of the outer hole pocket that are sampled by the various representations studied in this section. We will use this figure to understand the frequency dependence of the low energy inelastic response in the following section.

At this point, it is important to note certain circumstances where this procedure of extracting pure irreducible representations breaks down. First, in order to isolate the cross section from a pure irreducible representation, we note that the subtraction process requires some basic assumptions about the band structure and superconducting gaps. Such assumptions have been made in the context of the cuprates as well Devereaux and Hackl (2007); Devereaux and Einzel (1995); Devereaux et al. (1994, 1996). In our analysis, we used simple band structures created from invariants using basis functions in Table 1; but once a more accurate tight-binding band structure calculation emerges, our subtraction procedure can be used to analyse inelastic Raman data with greater confidence. More complications could arise when there are multiple bands crossing the Fermi level (as the case is here) due to the fact that a pure IRRP in one band need not correspond to the same IRRP in the other. In such cases, one could still think of indirect ways to avoid mixing of different representations. For example, on could rely on the fact that, in general, the gaps on different bands have different magnitudes. With an approximate knowledge of the values of the gaps on each band (from experiments like ARPES), it is possible to isolate the response from an individual band by keeping only the response in the relevant energy window and subtracting the contributions not in that window; this could be done for each polarization geometry separately. Although this procedure is less straighforward, one could avoid mixing of different IRRPs to a certain extent. However, this procedure becomes ambiguous when the gaps are similar in magnitude and/or have different form factors on different bands, and could have limited applicability. In such cases, one has to be satisfied with studying the responses with multiple contributing IRRPs (all be it with knowledge of what IRRPs form the combination). One could still make reasonably accurate comparisons with the corresponding theoretical predictions of the responses in mixed geometries Second, in the presence of vertex corrections, there are additional contributions that need to be incorporated into the inelastic response Devereaux et al. (1996); Maiti et al. (2017) which might, in general, make it difficult to isolate pure representations. Third, the effective mass approximation used in our analysis is valid only when the frequency of the incoming and scattered light is either small or large in comparison to the inter-band spacings Devereaux et al. (1996); Setty and Hu (2014).

This may not always be accurate Setty and Hu (2014) and one must be careful in not overestimating the vertex in certain regions of the Fermi surface. Due to the above outlined reasons, our calculations should not be considered any more than qualitative.

IV.3 Low energy Raman response

Having outlined the procedure to isolate the response from select irreducible representations above, we now demonstrate its use in determining the exact form of the pairing symmetry in POS. We confine ourselves to the spin singlet, even parity forms of the gap functions proposed in Goryo’s work Goryo (2003). Similar arguments can be carried out in a straightforward manner for other pairing forms. We choose the temperature and broadening parameter in the response function to be 2 meV and 3 meV respectively. Fig 6 (right) shows a color plot of the $A_{g}$ gap function, $\Delta(\vec{k})=\Delta_{0}\Phi_{s}(\vec{k})$, proposed in the A phase of POS. The brighter (green) regions show the location of the six nodes on each face of the cube. Fig 6 (left) shows the low energy response for such a gap function in channels belonging to different irreducible representations. A suppression of spectral weight for small $\omega$ is observed in the $T_{g}(3)$ ($xy$) geometry, with the $E_{g}(1)$ and $E_{g}(2)$ geometries showing enhancement for small $\omega$. To understand this result, we see from Fig 5 that the Fermi surface averaging in the $T_{g}(3)$ channel occurs mainly along the $x=\pm y$ directions on the $x-y$ plane. It is exactly along these directions where the chosen $A_{g}$ gap function makes the Fermi surface fully gapped, thus, suppressing low energy quasiparticle response. The $E_{g}(1)$ and $E_{g}(2)$ geometries, on the other hand, sample at least two of the gap nodes (see Fig 5) and excite a low energy quasi-particle response. To confirm this picture, we repeated the calculation for a four-dip state given by $\Delta(\vec{k})=\Delta_{0}(k_{x}^{2}k_{y}^{2}+k_{y}^{2}k_{z}^{2}+ak_{z}^{2})$, which has nodes only along the $k_{x}$ and $k_{y}$ axes. In this case, in addition to the $T_{g}(3)$ geometry, a small suppression of quasiparticle spectral weight occurs in the $\Gamma_{E_{g}}(2)$ geometry as well. The suppression is only small due to the fact that the $E_{g}(2)$ geometry $-$inspite of its extended lobe-like feature along the $k_{z}$ axis$-$ has a small vertex contribution coming from the equatorial plane as well (see fig 5 bottom right).

We can follow a similar analysis in the B phase. For the sake of convenience, we choose a gap function, $\Delta(\vec{k})=\left(\Phi_{s}(\vec{k})+i\Phi_{d}(\vec{k})\right)\equiv\Delta_{s0}(k_{x}^{2}k_{y}^{2}+k_{y}^{2}k_{z}^{2}+k_{z}^{2}k_{x}^{2})+i\Delta_{d0}(k_{x}^{2}-k_{y}^{2})$, such that there are two nodes along the $k_{z}$ axis (such a choice is possible provided the basis functions are also transformed appropriately). Given that the intensity of the Raman cross section depends only on the absolute value of the gap function, it is insensitive to the total phase of the order parameter but is sensitive to relative phases of the individual components. Fig (7) (right) shows the gap nodal structure in the B phase $-$ the Fermi surface now has nodes along the $k_{z}$ direction and is gapped everywhere else. Following arguments presented in the previous paragraph, we see that in the $T_{g}(3)$ and $E_{g}(1)$ geometries, which represents sample regions primarily along the $k_{x}=\pm k_{y}$ and $k_{x}=0,k_{y}=0$ directions respectively, the low energy response is reduced (see Fig (7) left panel, and Fig 8 which compares the $E_{g}(1)$ response for different pairing forms). However, in the $E_{g}(2)$ geometry, the gap function along the $k_{z}$ axis is projected; as a result, one should expect an increased low energy response. In the $A_{g}$ geometry, all the high symmetry regions are averaged over and results in a response intermediate in curvature to the other geometries. The four-dip state in the TRSB has been argued Goryo (2003) to be highly accidental and hence we will not discuss this state. However, if indeed, this state were to be realized, one would obtain a low energy increase of the intensity in the $\Gamma_{E_{g}}(1)$ geometry as well.

At this juncture, we would like to emphasize what we briefly pointed out earlier regarding the response in the $E_{g}(2)$ geometry. For the gap with four nodes, this geometry exhibits an intermediate behavior between enhancement and suppression. This is because there is also an equatorial contribution, although small, that cannot be entirely neglected compared to the axial contribution (see Fig 5). One should, therefore, in the process of trying to extract nodal behavior in gap functions, interpret the $E_{g}(2)$ response with care.

V Final Remarks

We derived a Ginzburg-Landau theory starting from a microscopic model that describes the entry of the skutterudite superconductor $\mathrm{PrOs_{4}Sb_{12}}$ (POS) into a time reversal symmetry broken (TRSB) phase as a function of temperature, in accordance with recent AC susceptibility and polar Kerr effect measurements Levenson-Falk et al. (2016). Using the expansion, we calculated the GL coefficients using their relationships to the various microscopic parameters and determined the shape of the contour bounding the TRSB phase. As an experimental means to probe the proposed pairing symmetries, we examined the low energy inelastic (Raman) response in both the A and B phases of POS. We provided a qualitative understanding of the low energy and zero momentum transfer response for various light polarization geometries. By appropriate manipulation of the incoming and scattered light geometries, as well as other subtraction procedures, we demonstrated that one can access the various irreducible representations contained in the $T_{h}$ point group describing POS. Depending on the existence of nodes along the $c-$ axis, we found an enhancement (nodal) or suppression (gapped) of the low energy spectral weight, based on what regions of the Fermi surface were being sampled. Inelastic light scattering could, thus, help pin-point the exact location of nodes on the fermi surface and, thereby, narrow down the candidate pairing symmetries proposed in POS.

Acknowledgments: CS and PWP acknowledge support from Center for Emergent Superconductivity, a DOE Energy Frontier Research Center, Grant No. DE-AC0298CH1088. YW is supported by the Gordon and Betty Moore Foundations EPiQS Initiative through Grant No. GBMF4305.

VI Appendix: Landau-Ginzburg Analysis

Here we outline some details about calculations involved in obtaining the phase diagrams appearing in Fig 3 of the main text. We will start by repeating the following form of the action written in momentum space:

[TABLE]

Here $\bar{c}_{\alpha\sigma}(\textbf{k})$ and $c_{\beta\sigma}(\textbf{k})$ are Grassmann numbers that follow Grassmannian algebra for electrons with Bloch momentum $\vec{k}$ and Matsubara frequency $k_{n}$ (collectively denoted by k), spin $\sigma$ and orbital indices $\alpha,\beta$, $\epsilon_{\alpha\beta}(\vec{k})$ are the orbital matrix elements, $\beta$ is the inverse temperature (not to be confused with the index $\beta$), $g$ is the interaction strength which is taken to be a constant, and $\tilde{V}(\vec{k},\vec{q})$ is the four body interaction giving rise to superconductivity. Changing to the band basis with energy eigenvalues $\epsilon_{\alpha}(\vec{k})$ by using the unitary transformation, $c_{\alpha\sigma}(\textbf{k})=U_{\alpha\beta}(\vec{k})\psi_{\beta\sigma}(\textbf{k})$, gives us the action

[TABLE]

where $V(\vec{k},\vec{q})$ is the resulting interaction related to $\tilde{V}(\vec{k},\vec{q})$ through the band matrix elements. As we are interested in zero-momentum pairing, we will assume that $\textbf{k}^{\prime}=-\textbf{k}$. After making the shift $\textbf{q}\rightarrow\textbf{q}-\textbf{k}$, we obtain the standard BCS action for independent bands

[TABLE]

where the interaction has been factorized in terms of $D(\vec{k})\equiv\phi_{s}(\vec{k})+\phi_{d}(\vec{k})$, which is the total form factor and is a mixture of the $s-$ wave and $d-$wave form factors. The total partition function can be written as

[TABLE]

We can now decouple the quartic term appearing in the partition function using the Hubbard-Stratonovich transformation for fermions given by

[TABLE]

where $a$ and $b$ are two commuting elements of a Grassmann algebra, defined by the numbers $\psi,\bar{\psi}$ which follow the anti-commutation relations $\{\psi_{i},\psi_{j}\}=\{\psi_{i},\bar{\psi}_{j}\}=\{\bar{\psi_{i}},\bar{\psi_{j}}\}=0$ with $i$ and $j$ sets of electron quantum numbers. In our case, we have defined $a\equiv\bar{\psi}_{\alpha\uparrow}(\textbf{q})\bar{\psi}_{\alpha\downarrow}(-\textbf{q})D(\vec{q})$ and $b\equiv\psi_{\alpha\downarrow}(-\textbf{k})\psi_{\alpha\uparrow}(\textbf{k})D(\vec{k})$ and it is easily seen that $a$ and $b$ commute. Using this, we can rewrite the partition function as

[TABLE]

which can be further simplified using matrix notation as

[TABLE]

Here we have defined

[TABLE]

along with $\hat{\Psi}_{\textbf{k}\alpha}=\left(\bar{\psi}_{\alpha\uparrow}(\textbf{k}),\psi_{\alpha\downarrow}(-\textbf{k}).\right)$. Performing the Fermionic functional integral, we obtain

[TABLE]

In order to perform an expansion in powers of the superconducting order parameter, we write the total Green function inverse as a sum of non-interacting and superconducting self energy parts; that is, we write $\mathscr{G}_{\textbf{k}\alpha}^{-1}=\hat{\mathscr{G}}_{0\textbf{k}\alpha}^{-1}+\hat{\Sigma}$, where,

[TABLE]

and

[TABLE]

We now perform an expansion keeping in mind that $Log\left[\hat{\mathscr{G}}_{\textbf{k}\alpha}^{-1}\right]=Log\left[\hat{\mathscr{G}}_{0\textbf{k}\alpha}^{-1}+\hat{\Sigma}\right]=Log\left[\hat{\mathscr{G}}_{0\textbf{k}\alpha}^{-1}\right]+\sum_{n=1}^{\infty}\frac{(-1)^{n+1}}{n}\left(\hat{\mathscr{G}}_{0\textbf{k}\alpha}\hat{\Sigma}\right)^{n}$. All the powers odd in $n$ vanish because of the Matsubara sums and we are only left with even powers in $n$. We will only be interested in second and fourth order powers of $\Sigma$ in the Landau-Ginzburg expansion. With this, the superconducting contributions to the free energy become

[TABLE]

We decompose the pairing terms into two channels: the $s-$ wave and $d-$ wave symmetries according to the point group symmetries of the lattice. For the second order contribution, we have

[TABLE]

and the fourth order contribution is given by

[TABLE]

Note that here $\phi_{s}$ and $\phi_{d}$ have the respective momentum dependences from the point group of the crystal. Thus, we define $\Delta_{s}$ and $\Delta_{d}$ by $\phi_{s,d}\equiv\Delta_{s,d}\Phi_{s,d}(\vec{k})$ with the functions $\Phi_{s,d}(\vec{k})$ containing all the momentum dependence. We now evaluate these contributions to the free energy by expanding into the quadratic and quartic powers of $\Delta_{s}$ and $\Delta_{d}$ to give

[TABLE]

with the following definitions of the coefficients (with $g_{\pm}\equiv(g_{s}\pm g_{d})/2$)

[TABLE]

Given that $\alpha^{\prime}$ and $g_{\pm}$ are small compared to the other coefficients due to negligible overlap between odd powers of $\Phi_{s}(\vec{k})$ and $\Phi_{d}(\vec{k})$, and for simplicity in minimizing the free energy, it is easier to study the problem in the absence of linear couplings of $\Delta_{s}$ and $\Delta_{d}$, i.e we ignore terms in the free energy proportional to $(\Delta_{s}^{*}\Delta_{d}+\Delta_{s}\Delta_{d}^{*})$. Defining a phase difference between the $d-$wave and $s-$wave gaps to be equal to $\eta$, the free energy is minimized when $\eta=\pi/2$. Minimization of $\mathscr{F}_{s}$ with respect to $|\Delta_{s}|$ and $|\Delta_{d}|$ yields two sets of equations for $\Delta_{s}$ and $\Delta_{d}$ given by (defining $\alpha_{s,d}=\alpha\pm\delta\alpha$)

[TABLE]

Solving these equations gives us four combinations of solutions $\left(|\Delta_{s}|,|\Delta_{d}|\right)=(0,0)$, $\left(|\Delta_{s}|,|\Delta_{d}|\right)=(0,\sqrt{\frac{-(\alpha-\delta\alpha)}{2\beta_{d}}})$, $\left(|\Delta_{s}|,|\Delta_{d}|\right)=(\sqrt{\frac{-(\alpha+\delta\alpha)}{2\beta_{s}}},0)$ and $\left(|\Delta_{s}|,|\Delta_{d}|\right)=\left(\sqrt{\frac{\alpha(\beta_{1}-\beta_{d})-\delta\alpha(\beta_{1}+\beta_{d})}{2(\beta_{s}\beta_{d}-\beta_{1}^{2})}},\sqrt{\frac{\alpha(\beta_{s}-\beta_{1})-\delta\alpha(\beta_{s}+\beta_{1})}{2(\beta_{1}^{2}-\beta_{s}\beta_{d})}}\right)$. The last pair of solutions which gives rise to a time reversal symmetry broken (TRSB) phase ($|\Delta_{s}|+i|\Delta_{d}|$) has the lowest free energy but is obviously a stable solution only when $(|\Delta_{s}|,|\Delta_{d}|)$ are both real. Hence, the conditions $\frac{\alpha(\beta_{1}-\beta_{d})-\delta\alpha(\beta_{1}+\beta_{d})}{2(\beta_{s}\beta_{d}-\beta_{1}^{2})}>0$ and $\frac{\alpha(\beta_{s}-\beta_{1})-\delta\alpha(\beta_{s}+\beta_{1})}{2(\beta_{1}^{2}-\beta_{s}\beta_{d})}>0$ define the TRSB phase. The $s-$wave and $d-$wave pairings are stable when $-(\alpha+\delta\alpha)$ and $-(\alpha-\delta\alpha)$ are positive respectively. These phases are shown in Fig. 3 (panels a-c) of the main text. The presence of the linear terms can be treated perturbatively and yields Fig 3 (d) as described in the main text when $\beta_{s}=\beta_{d}$.