Low energy inelastic response in the superconducting phases of $\mathbf{PrOs_{4}Sb_{12}}$

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.

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

	$\displaystyle\epsilon_{1}(\vec{k})$	$\displaystyle=$	$\displaystyle\epsilon_{1}+A_{g}(\vec{k})-\mu$
	$\displaystyle\epsilon_{2}(\vec{k})$	$\displaystyle=$	$\displaystyle\epsilon_{2}+A_{g}(\vec{k})-2\vec{E}_{g}(\vec{k})\cdot\vec{E}_{g}(\vec{k})-\mu$
	$\displaystyle\epsilon_{3}(\vec{k})$	$\displaystyle=$	$\displaystyle\epsilon_{3}-\vec{T}_{u}(\vec{k})\cdot\vec{T}_{u}(\vec{k})-\mu$		(1)

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

	$\displaystyle A_{g}(\vec{k})$	$\displaystyle=$	$\displaystyle 2(c_{x}+c_{y}+c_{z})$
	$\displaystyle\vec{E}_{g}(\vec{k})$	$\displaystyle=$	$\displaystyle 2\left(c_{x}-c_{y},1-c_{z}-\frac{1}{2}(2-c_{x}-c_{y})\right)$
	$\displaystyle\vec{T}_{u}(\vec{k})$	$\displaystyle=$	$\displaystyle 4\left(s_{\frac{x}{2}}(c_{y}-c_{z}),s_{\frac{y}{2}}(c_{z}-c_{x}),s_{\frac{z}{2}}(c_{x}-c_{y})\right).$		(2)

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.

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

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

$\displaystyle=$

$\displaystyle\beta\sum_{\begin{subarray}{c}\alpha\beta\\ \textbf{k}\sigma\end{subarray}}\bar{c}_{\alpha\sigma}(\textbf{k})\left(-ik_{n}\delta_{\alpha\beta}-\epsilon_{\alpha\beta}(\vec{k})\right)c_{\beta\sigma}(\textbf{k})-\beta g\sum_{\begin{subarray}{c}\textbf{k}\textbf{k}^{\prime}\alpha\end{subarray}}\bar{c}_{\alpha\uparrow}(\textbf{k})D(\vec{k})\bar{c}_{\alpha\downarrow}(-\textbf{k})c_{\alpha\downarrow}(-\textbf{k}^{\prime})D(\vec{k}^{\prime})c_{\alpha\uparrow}(\textbf{k}^{\prime}).$

(3)

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))

	$\displaystyle\mathscr{F}_{s}$	$\displaystyle=$	$\displaystyle\alpha_{s}|\Delta_{s}|^{2}+\alpha_{d}|\Delta_{d}|^{2}+\beta_{s}|\Delta_{s}|^{4}+\beta_{d}|\Delta_{d}|^{4}+4\beta_{1}|\Delta_{s}|^{2}|\Delta_{d}|^{2}$
			$\displaystyle+\alpha^{\prime}(\Delta_{s}^{*}\Delta_{d}+\Delta_{s}\Delta_{d}^{*})+\beta_{2}\left(\Delta_{s}^{2}\Delta_{d}^{*2}+\Delta_{s}^{*2}\Delta_{d}^{2}\right)+2\sum_{\nu=\pm}g_{\nu}(|\Delta_{s}|^{2}+\nu|\Delta_{d}|^{2})(\Delta_{s}\Delta_{d}^{*}+2\Delta_{s}^{*}\Delta_{d}),$

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

	$\displaystyle D_{s}$	$\displaystyle=$	$\displaystyle\sqrt{\frac{\alpha(\beta_{1}-\beta_{d})-\delta\alpha(\beta_{1}+\beta_{d})}{2(\beta_{s}\beta_{d}-\beta_{1}^{2})}}$		(4)
	$\displaystyle D_{d}$	$\displaystyle=$	$\displaystyle\sqrt{\frac{\alpha(\beta_{s}-\beta_{1})-\delta\alpha(\beta_{s}+\beta_{1})}{2(\beta_{1}^{2}-\beta_{s}\beta_{d})}}.$		(5)

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).

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

	$\displaystyle\frac{\partial^{2}\sigma}{\partial\omega\partial\Omega}$	$\displaystyle=$	$\displaystyle\frac{\omega_{S}}{\omega_{I}}r_{0}^{2}S_{\gamma\gamma}(\vec{q},\omega)$		(6)
	$\displaystyle S_{\gamma\gamma}(\vec{q},\omega)$	$\displaystyle=$	$\displaystyle-\frac{1}{\pi}[1+n_{B}(\omega)]Im\chi_{\gamma\gamma}(\vec{q},\omega),$		(7)

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

$$\chi_{\gamma\gamma}(\omega)=\int_{0}^{\beta}d\tau e^{-i\omega_{m}\tau}\langle T_{\tau}\tilde{\rho}_{\gamma}(\tau),\tilde{\rho}_{\gamma}(0)\rangle|_{i\omega_{m}\rightarrow\omega+i\delta}$$

(8)

where the effective (weighted) density is written as

$$\tilde{\rho_{\gamma}}=\sum_{\vec{k},\sigma}\sum_{n,m}\gamma_{n,m}(\vec{k})c_{n,\sigma}^{\dagger}(\vec{k})c_{m,\sigma}(\vec{k}),$$

(9)

$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

	$\displaystyle\gamma_{n}(\vec{k})$	$\displaystyle=$	$\displaystyle\hat{e}^{i}\hat{e}^{s}+\frac{1}{m}\sum_{j\neq n}\frac{\langle n,\vec{k}|\hat{e}^{s}p|j,\vec{k}\rangle\langle j,\vec{k}|\hat{e}^{i}p|n,\vec{k}\rangle}{\epsilon_{n}(\vec{k})-\epsilon_{j}(\vec{k})+\omega_{I}}$
			$\displaystyle+\frac{\langle n,\vec{k}|\hat{e}^{i}p|j,\vec{k}\rangle\langle j,\vec{k}|\hat{e}^{s}p|n,\vec{k}\rangle}{\epsilon_{n}(\vec{k})-\epsilon_{j}(\vec{k})-\omega_{S}}.$

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

$$\gamma_{n}(\vec{k})=\sum_{\mu}\gamma_{n}^{\mu}\mathscr{I}_{n}^{\mu}(\vec{k})$$

(10)

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)

$\displaystyle\chi_{\gamma\gamma}^{(n)}(i\omega)$

$\displaystyle=$

$\displaystyle\sum_{\vec{k}}|\gamma_{n}(\vec{k})|^{2}\lambda_{n}(\vec{k},i\omega)$

(11)

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

$$\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),$$

(12)

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.

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.

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:

	$\displaystyle S(\bar{c},c)$	$\displaystyle=$	$\displaystyle\beta\sum_{\begin{subarray}{c}\alpha\beta\\ \textbf{k}\sigma\end{subarray}}\bar{c}_{\alpha\sigma}(\textbf{k})\left(-ik_{n}\delta_{\alpha\beta}-\epsilon_{\alpha\beta}(\vec{k})\right)c_{\beta\sigma}(\textbf{k})$		(23)
			$\displaystyle-\beta g\sum_{\begin{subarray}{c}\textbf{k}\textbf{k}^{\prime}\\ \textbf{q}\alpha\end{subarray}}\bar{c}_{\alpha\uparrow}(\textbf{k}+\textbf{q})\bar{c}_{\alpha\downarrow}(\textbf{k}^{\prime}-\textbf{q})\tilde{V}(\vec{k},\vec{q})c_{\alpha\downarrow}(\textbf{k}^{\prime})c_{\alpha\uparrow}(\textbf{k}).$

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

	$\displaystyle S(\bar{\psi},\psi)$	$\displaystyle=$	$\displaystyle\beta\sum_{\begin{subarray}{c}\alpha\\ \textbf{k}\sigma\end{subarray}}\bar{\psi}_{\alpha\sigma}(\textbf{k})\left(-ik_{n}-E_{\alpha}(\textbf{k})\right)\psi_{\alpha\sigma}(\textbf{k})$		(24)
			$\displaystyle-\beta g\sum_{\begin{subarray}{c}\textbf{k}\textbf{k}^{\prime}\\ \textbf{q}\alpha\end{subarray}}\bar{\psi}_{\alpha\uparrow}(\textbf{k}+\textbf{q})\bar{\psi}_{\alpha\downarrow}(\textbf{k}^{\prime}-\textbf{q})V(\vec{k},\vec{q})\psi_{\alpha\downarrow}(\textbf{k}^{\prime})\psi_{\alpha\uparrow}(\textbf{k}),$

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

	$\displaystyle S(\bar{\psi},\psi)$	$\displaystyle=$	$\displaystyle\beta\sum_{\begin{subarray}{c}\alpha\\ \textbf{k}\sigma\end{subarray}}\bar{\psi}_{\alpha\sigma}(\textbf{k})\left(-ik_{n}-E_{\alpha}(\textbf{k})\right)\psi_{\alpha\sigma}(\textbf{k})$		(25)
			$\displaystyle-\beta g\sum_{\begin{subarray}{c}\textbf{k}\textbf{q}\alpha\end{subarray}}\bar{\psi}_{\alpha\uparrow}(\textbf{q})\bar{\psi}_{\alpha\downarrow}(-\textbf{q})D(\vec{q})D(\vec{k})\psi_{\alpha\downarrow}(-\textbf{k})\psi_{\alpha\uparrow}(\textbf{k}),$

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

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

(26)

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

$$\pi ge^{gab}=\int d\phi d\bar{\phi}\hskip 8.53581pte^{a\phi+b\bar{\phi}}e^{\frac{-|\phi|^{2}}{g}}$$

(27)

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

	$\displaystyle\mathscr{Z}$	$\displaystyle=$	$\displaystyle\int\mathscr{D}(\psi,\bar{\psi})\mathscr{D}(\phi,\bar{\phi})exp\left[-\frac{\beta|\phi|^{2}}{g}-\beta\sum_{\begin{subarray}{c}\textbf{k}\alpha\sigma\end{subarray}}\bar{\psi}_{\alpha\sigma}(\textbf{k})\left(-ik_{n}-E_{\alpha}(\vec{k})\right)\psi_{\alpha\sigma}(\textbf{k})\right]$		(28)
			$\displaystyle\times exp\left[\beta\sum_{\begin{subarray}{c}\textbf{k}\alpha\end{subarray}}\left(\phi\bar{\psi}_{\alpha\uparrow}(\textbf{k})\bar{\psi}_{\alpha\downarrow}(-\textbf{k})+\phi^{*}\psi_{\alpha\downarrow}(-\textbf{k})\psi_{\alpha\uparrow}(\textbf{k})\right)D(\vec{k})\right],$

which can be further simplified using matrix notation as

$$\mathscr{Z}=\int\mathscr{D}(\psi,\bar{\psi})\mathscr{D}(\phi,\bar{\phi})exp\left[-\beta\left(\frac{|\phi|^{2}}{g}+\sum_{\textbf{k}\alpha}\hat{\Psi}^{\dagger}_{\textbf{k}\alpha}\hat{\mathscr{G}}^{-1}_{\textbf{k}\alpha}\hat{\Psi}_{\textbf{k}\alpha}\right)\right].$$

(29)

Here we have defined

$$\hat{\mathscr{G}}^{-1}_{\textbf{k}\alpha}=\begin{pmatrix}-\left(ik_{n}+E_{\alpha}(\vec{k})\right)&-D(\vec{k})\phi\\ -D(\vec{k})\phi^{*}&-\left(ik_{n}-E_{\alpha}(\vec{k})\right)\end{pmatrix}$$

(30)

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

$$\mathscr{Z}=\int\mathscr{D}(\phi,\bar{\phi})exp\left[\frac{-\beta|\phi|^{2}}{g}+\sum_{\textbf{k}\alpha}Tr\left(Log\left[\hat{\mathscr{G}}_{\textbf{k}\alpha}^{-1}\right]\right)\right].$$

(31)

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,

$$\hat{\mathscr{G}}^{-1}_{0\textbf{k}\alpha}=\begin{pmatrix}-\left(ik_{n}+E_{\alpha}(\vec{k})\right)&0\\ 0&-\left(ik_{n}-E_{\alpha}(\vec{k})\right)\end{pmatrix}$$

(32)

and

$$\hat{\Sigma}=\begin{pmatrix}0&-D(\vec{k})\phi\\ -D(\vec{k})\phi^{*}&0.\end{pmatrix}$$

(33)

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

$$\mathscr{F}_{s}\simeq\frac{|\phi|^{2}}{g}+\frac{1}{\beta}\sum_{\textbf{k}\alpha}\left(\frac{1}{2}Tr\left[\left(\hat{\mathscr{G}}_{0\textbf{k}\alpha}\hat{\Sigma}\right)^{2}\right]+\frac{1}{4}Tr\left[\left(\hat{\mathscr{G}}_{0\textbf{k}\alpha}\hat{\Sigma}\right)^{4}\right]\right).$$

(34)

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

$$\mathscr{F}_{s}^{(2)}=T\sum_{\textbf{k}\alpha}\frac{|\phi_{s}+\phi_{d}|^{2}}{\left(E_{\alpha}(\vec{k})^{2}+k_{n}^{2}\right)},$$

(35)

and the fourth order contribution is given by

$$\mathscr{F}_{s}^{(4)}=\frac{T}{2}\sum_{\textbf{k}\alpha}\frac{|\phi_{s}+\phi_{d}|^{4}}{\left(E_{\alpha}(\vec{k})^{2}+k_{n}^{2}\right)^{2}}.$$

(36)