Interacting spin-3/2 fermions in a Luttinger (semi)metal: competing phases and their selection in the global phase diagram

We compute the effects of electronic interactions on gapless spin-3/2 excitations that in a noninteracting system emerge at a bi-quadratic touching of Kramers degenerate valence and conduction bands, known as Luttinger semimetal. This model can describe the low-energy physics of HgTe, gray-Sn, 227 pyrochlore iridates and half-Heuslers. For the sake of concreteness we only consider the short-range components of the Coulomb interaction. By combining mean-field analysis with a renormalization group (RG) calculation, we construct multiple cuts of the global phase diagram of interacting spin-3/2 fermions at zero and finite temperature and chemical doping. Such phase diagrams display a rich confluence of competing orders, among which rotational symmetry breaking nematic insulators and time reversal symmetry breaking magnetic orders are the prominent excitonic phases. We also show that even repulsive interactions can be conducive for both s-wave and d-wave pairings. The reconstructed band structure inside the ordered phases allows us to organize them according to the energy (entropy) gain in the following (reverse) order: s-wave pairing, nematic phases, magnetic orders and d-wave pairings, at zero chemical doping. But, the paired states are always energetically superior over the excitonic ones for finite doping. The phase diagrams obtained from the RG analysis show that an ordered phase with higher energy (entropy) gain is realized at low (high) temperature. In addition, we establish a"selection rule"between the interaction channels and the resulting ordered phases, suggesting that repulsive interactions in the magnetic (nematic) channels are conducive for the nucleation of d-wave (s-wave) pairing. The proposed methodology can shed light on the global phase diagram of other strongly interacting multi-band systems, such as doped Dirac semimetal, topological insulators and the like.


I. INTRODUCTION
The discovery of new phases of matter, and the study of the transitions between them, form the core of condensed matter and materials physics [1]. Much experimental ingenuity is devoted to realising and controlling different tuning parameters in the laboratorysuch as temperature, pressure, magnetic field, chemical composition-to drive a system from one phase to the other [1]. Most simply, the appearance of different phases can often be appreciated from the competition between energy and entropy. For example, with decreasing temperature transitions from water vapor to liquid to ice are intimately tied with the reduction of entropy or gain in energy. Similarly, in a metal, entropy of gapless excitations on the Fermi surface is exchanged for condensation energy as a superconducting gap opens [2].
With increasing complexity of quantum materials attained in recent decades, the richness of the global phase diagram of strongly correlated materials has amplified enormously. And various prototypical representatives, such as cuprates, pnictides, heavy fermion compounds display concurrence of competing orders, among which spin-and charge-density-wave, superconductivity, nematicity are the most prominent ones. Besides establishing the existence of -and, hopefully, eventually utilizing in technological applications -these phenomena, an obvious challenge is to discover any simplifying per-spective, or at least heuristic classification scheme, for predicting and classifying them: do there exist any organising principles among multiple competing orders that can shed light on the global phase diagram of strongly correlated materials? Restricting ourselves to a specific but remarkably rich metallic system, we here give a partially affirmative answer to this question. We study a collection of strongly interacting spin-3/2 fermions in three dimensions that in the normal phase display a biquadratic touching of Kramers degenerate valence and conduction bands at an isolated point in the Brillouin zone, see Fig. 1. This system is also known as Luttinger (semi)metal [3,4].
Such peculiar quasiparticle excitations can be found for example in HgTe [5], gray-Sn [6,7], 227 pyrochlore iridates (Ln 2 Ir 2 O 7 , where Ln is the lanthanide element) [8][9][10][11][12] and half-Heusler compounds (such as LnPtBi, LnPdBi) [13][14][15]. Bi-quadratic band touching has recently been observed in the normal state of Pr 2 Ir 2 O 7 [10] and Nd 2 Ir 2 O 7 [11] via angle resolved photo emission spectroscopy (ARPES). While most of the iridium-based oxides support all-in all-out arrangement of 4d Ir-electrons [16][17][18][19][20][21][22][23][24], a singular member of this family, namely Pr 2 Ir 2 O 7 , possibly resides at the brink of a metallic spin-ice or 3-in 1-out ordering and supports a large anomalous Hall conductivity [25][26][27][28]. While these materials harbor various competing magnetic orders (due to comparable Hubbard and spin-orbit inter- The red dot represents the Luttinger semimetal fixed point separating electron-and hole-doped Fermi liquid of spin-3/2 fermions. The shaded sector corresponds to the quantum critical regime associated with the z = 2 fixed point in d = 3 (the red dot), controlling scaling properties in this regime (see Sec. II A). Above a nonuniversal energy high energy cutoff EΛ (red dashed line), spin-3/2 fermions lose their jurisdiction. The crossover temperature T * (black dashed lines) above which critical Luttinger fermions are operative is given by T * ∼ |n| z/d , with z = 2 and d = 3, where n is the carrier density measured from the bi-quadratic band touching point. In this work we demonstrate the role of electron-electron interactions on the quantum critical regime (the shaded region), at finite temperature and chemical doping (see Sec. II B for summary of our results).
The effects of electronic interactions on spin-3/2 fermions are addressed within the framework of an extended Hubbardlike model, composed of only the short-range components of repulsive Coulomb interactions. 1 Due to the vanishing density of states (DoS) in a Luttinger system (namely, (E) ∼ √ E), sufficiently weak short-range interactions are irrelevant perturbations. Therefore, any ordering takes place at finite coupling. We here employ a renormalization group (RG) analysis to construct various cuts of the global phase diagram at zero as well as finite chemical doping [see Figs. 3,4 and 5], and combine it with mean-field analysis to gain insight into the organizing principle among distinct broken-symmetry phases (BSPs) 2 . A gist of our findings can be summarized as follows.
1. By computing the reconstructed band structure (within the mean-field approximation) we organize dominant BSPs according to the condensation energy and entropy gains. Results are summarized in Fig. 2. We note that while the stiffness (uniform or anisotropic) of the band gap measures the condensation energy gain, the amount of gapless quasiparticles (resulting in powerlaw scaling of DoS at low energies) measures the entropy inside the ordered phase. All cuts of the global phase diagram (obtained from a RG analysis) show that the low (high) temperature phase yields larger condensation energy (entropy) gain, see Fig. 3.
2. The quasiparticle spectra inside the s-wave paired state and two nematic phases (belonging to the T 2g and E g representations of the cubic or O h point group) are fully gapped. Hence, nucleation of these three phases leads to the maximal gain of condensation energy, and they appear as the dominant BSPs at zero temperature, as shown in Fig. 3.
3. At finite temperature condensation energy gain competes with the entropy, and phases with higher entropy are realized at higher temperatures. Onset of any nematicity results in an anisotropic gap, in contrast to the situation with an s-wave pairing. Thus former orderings are endowed with larger (smaller) entropy (condensation energy), and are found at higher temperatures in Figs. 3(a), 3(b). By contrast, the dominant magnetic orders (belonging to the A 2u and T 1u representations) produce gapless Weyl quasiparticles and re- 1 In this work we neglect the long-range tail of the Coulomb interaction. When the chemical potential is pinned at the band touching point, long-range Coulomb interaction can give rise to an infrared stable non-Fermi liquid fixed point [39,60,61], which, may however be unstable toward the formation of an excitonic phase [40,42,62]. By contrast, at finite chemical doping longrange interaction suffers conventional Thomas-Fermi screening and becomes short-ranged in nature. A systematic incorporation of the long-range interaction into our discussion is left for a future investigation. sult in (E) ∼ |E| 2 scaling of the DoS at low energies. Hence, these two magnetic orders carry larger entropy than the nematic phases or the s-wave pairing, and can only be found at finite temperature, see Figs. 3(c), 3(d). Luttinger semimetal (LSM)-A 2u magnetic order (results from the all-in all-out state in pyrochlore lattice [8,20]) transition at finite temperature is consistent with the experimental observation in Nd 2 Ir 2 O 7 [16], while the T 1u magnetic ordering (yielding a 3-in 1-out ordering in pyrochlore lattice and supporting anomalous Hall effect [28]) can be germane for Pr 2 Ir 2 O 7 [25][26][27] 3 . 4. Local four fermion interactions in the nematic channels are conducive for s-wave pairing (at zero and finite chemical doping), see Fig. 3(a), 3(b) and 4 , whereas short-range magnetic interactions give birth to topological d-wave pairing (only at finite chemical doping), see Fig. 5. We further elaborate such an emergent "selection rule" in Sec. II B. Confluence of magnetic order and d-wave pairing (resulting in gapless BdG quasiparticles, found in YPtBi [37]) is in (qualitative) agreement with the global phase diagram of LnPdBi [36].
The theoretical approach outlined in this work is quite general and can be extended to address the effects of electronic interactions in various strongly correlated multiband systems, among which two-dimensional Dirac and quadratic fermions (respectively relevant for monolayer and bilayer graphene) [63,64], three-dimensional doped topological, crystalline and Kondo insulators [65][66][67], Weyl materials [68], twisted bilayer graphene [69][70][71], are the most prominent and experimentally pertinent ones. In the future we will systematically study these systems.

A. Outline
The rest of the paper is organized as follows. In the next section we present an extended summary of our main results. The low-energy description of the Luttinger model and its symmetry properties are discussed in Sec. III. In Sec. IV we discuss the reconstructed band structure inside various excitonic and superconducting phases. In Sec. V we introduce the interacting model for spin-3/2 fermions and analyze the propensity toward the formation of various orderings within a mean-field approximation. Sec. VI is devoted to a renormalization group (RG) analysis of interacting Luttinger fermions at zero and finite temperatures and chemical doping. We summarize the main results and highlight some future outlooks in Sec. VII. Additional technical details are relegated to appendices.

II. EXTENDED SUMMARY
Our starting point is a collection of spin-3/2 fermions for which the normal state is described by a bi-quadractic touching of Kramers degenerate valence and conduction bands. The corresponding Hamiltonian operator is [3,4] whered j s are five d-wave harmonics in three dimensions, Γ j s are five mutually anticommuting four dimensional Hermitian matrices, and m/[cos α (sin α)] is the effective mass for gapless excitations in the T 2g (E g ) orbitals in a cubic environment. Momentum k and chemical potential µ are measured from the band touching point. Additional details of this model are discussed in Sec. III and Appendix C. The mass anisotropy parameter (α) lies within the range 0 ≤ α ≤ π 2 [28]. But, for the sake of concreteness we restrict our focus on the isotropic system with α = π 4 . For discussion on the evolution of phase diagrams with varying α, see Secs. V A and VI, and Figs. 6 and 16. The LSM is realized as an unstable fixed point at µ = 0, the red dot in Fig. 1. This fixed point is characterized by the dynamic scaling exponent z = 2, determining the relative scaling between energy and momentum according to E ∼ |k| z . The chemical potential is a relevant perturbation at this fixed point, with the scaling dimension [µ] = 2. Hence, the correlation length exponent at this fixed point is ν = 1/2. Therefore, the LSM can be envisioned as a quantum critical point (QCP) separating electron-doped (for µ > 0, the brown region) and holedoped (for µ < 0, the green region) Fermi liquid phases, as shown in Fig. 1. Our discussion is focused on the quantum critical regime (the shaded region).
The crossover temperature (T * ) separating the quantum critical regime accommodating gapless spin-3/2 excitations from the Fermi liquid phases can be estimated in the following way where ξ is a characteristic length scale and |n| ∼ |µ|/ξ d is the carrier density. Two critical exponents (z and ν) and the dimensionality of the system (d = 3) determine the scaling of various thermodynamic (such as specific heat, compressibility) and transport (such as dynamic conductivity) quantities in this regime.

A. Critical scaling in noninteracting system
The free-energy density (up to an unimportant temperature (T ) independent constant) inside the critical regime is given by (setting k B = 1) where Li represents the polylogarithimic function. The specific heat of this system is given by for µ/T 1 (ensuring that the system resides inside the critical regime), where a ≈ 1.7244 and b ≈ 0.6049. From the above expression of the free-energy density we can also extract the scaling of compressibility (κ), given by where c ≈ 0.00989. Therefore, the presence of finite chemical doping does not alter the leading power-law scaling of physical observables, such as C V ∼ T 3/2 , κ ∼ √ T , but only provides subleading corrections, which are suppressed by a parametrically small quantity µ/T . 4 Also note that is a universal ratio, capturing the signature of a z = 2 quantum critical point in d = 3. A detailed derivation of this analysis is shown in Appendix A. Qualitatively similar sub-leading corrections are also found in the scaling of the dynamic conductivity, which we discuss now. Gauge invariance mandates that the conductivity (σ) must scale as σ ∼ ξ 2−d . Hence for a collection of z = 2 quasiparticles (such as the Luttinger fermions), σ ∼ √ T or √ ω in three spatial dimensions. Indeed we find that the Drude (Dr) component of the dynamic conductivity in the Luttinger system is given by [see Appendix B] where F Dr (x) is a monotonically increasing universal function of its argument [see Eq. (B6) and Fig. 21] and ω is the frequency. On the other hand, the inter-band (IB) component of the optical conductivity reads as Hence, inter-band component of the optical conductivity vanishes as √ ω as ω → 0 and the LSM can be identified as a power-law insulator.
Therefore, even when the chemical doping is finite there exists a wide quantum critical regime, shown in Fig. 1, where the scaling of thermodynamic and transport quantities are essentially governed by z = 2 quasiparticles, and the chemical potential provides only sub-leading corrections. Next we highlight the imprint of finite temperature and chemical doping on the global phase diagram of interacting spin-3/2 fermions.

B. Electron-electron interactions in a Luttinger (semi)metal
In this work we compute the effects of electron-electron interactions on Luttinger fermions, occupying the critical regime of the noninteracting fixed point, see Fig. 1. In this regime any short-range or local four-fermion interaction (λ) is an irrelevant perturbation, since 4 The scaling of specific heat and compressibility is determined by the dimensionality (d) and dynamic scaling exponent (z) according to C V ∼ T 1+d/z and κ ∼ T d/z , respectively. due to the vanishing DoS, namely (E) ∼ √ E. We use a RG analysis, tailored to address the effects of electronic interactions on Luttinger fermions in d = 3, constituting a z = 2 band structure, to arrive at various cuts of the global phase diagram of this system. If, on the other hand, temperature and chemical doping are such that the system resides inside a Fermi liquid phase, the notion of z = 2 nodal quasiparticles becomes moot and our RG analysis loses its jurisdiction. 5 Furthermore, we aug- 5 Note that in the presence of a Fermi surface the DoS is constant and the interaction coupling λ is dimensionless. Consequently a Fermi liquid becomes unstable toward the formation of a super-ment the RG analysis with an organizing principle based on the competition between energy and entropy. To this end we rely on the computation of the reconstructed band structure inside the ordered phases within the mean-field approximation. Subsequently, we also promote a "selection rule" among neighboring phases in the global phase diagram, originating purely from their algebraic or symmetry properties. For the sake of simplicity we concenconductor (often non-s-wave) even in the presence of repulsive electron-electron interactions, following the spirit of the Kohn-Luttinger mechanism [72][73][74][75][76][77]  trate on the isotropic system (α = π 4 ) in the following three subsections. Nonetheless, our results hold (at least qualitatively) for any arbitrary value of α, as summarized in the last subsection.

IIB1. Organizing principle: Emergent Topology & Energy-Entropy
Let us first promote an organizing principle among BSPs according to their contribution to the energy and entropy gain and anticipate their presence in the global phase diagram. In what follows we highlight the reconstructed band structure inside the dominant ordered phases within a mean-field or Hartree-Fock approximation, which by construction undermines the ordered parameter fluctuations. The emergent band topology is computed by diagonalizing an effective single-particle Hamiltonian, composed of the noninteracting Luttinger Hamiltonian and corresponding order parameter, see Sec. IV C for details.
Perhaps it is natural to anticipate that at zero temperature strong electronic interactions favor the phases that produce the largest spectral gap, as the onset of these ordered states offers maximal gain of condensation energy. In a LSM there are three candidate BSPs that yield fully gapped quasiparticle spectra: (a) an s-wave superconductor, producing a uniform and isotropic gap, and (b) two nematic orders (belonging to the T 2g and E g representations), producing anisotropic gaps. As shown in Fig. 3, only these three phases can be found in an isotropic and interacting LSM at zero temperature.
The energy-entropy competition, leading to an organizing principle among competing phases at finite temperature, can be appreciated from the scaling of the DoS or the stiffness (uniform or anisotropic) of the spectral gap. As mentioned above, the s-wave pairing and nematic orders respectively produce uniform and anisotropic gaps, whereas two magnetic orders, belonging to the A 2u and T 1u representations, respectively produce eight [20] and two [28] isolated simple Weyl nodes, around which the DoS vanishes as (E) ∼ |E| 2 for sufficienlty low energies. On the other hand, each copy of the d-wave pairings accommodates two nodal loops for which the low-energy DoS scales as (E) ∼ |E| [see Sec. IV D] [49,51,53,54,59]. Since we are interested in energy or temperature scales much smaller than the ultraviolet cutoff or bandwidth (|E| 1), the structure of the spectral gap (isotropic or anisotropic) and powerlaw scaling of DoS carry sufficient information to organize the ordered phases according to their contribution to condensation energy and entropy, summarized in Fig. 2. In brief, existence of more gapless points (resulting in higher DoS near E = 0) yields larger entropy, while a more uniform gap leads to higher gain in condensation energy. From various cuts of the global phase diagram at finite temperature, see Fig. 3, we note the following common feature: Among competing orders, the one with maximal gain in condensation energy appears at low temperature, while the phase with higher entropy is realized at higher temperature, in accordance with the general principle of energy-entropy competition. Since the DoS in a LSM scales as (E) ∼ √ E (maximal entropy), it can always be found at sufficiently high (weak) enough temperature (interactions).
A similar conclusion can also be arrived at when the chemical potential is placed away from the band touching point. At finite chemical doping all particle-hole (two nematic and two magnetic) orders produce a Fermi surface (according to the Luttinger theorem [78]) and hence a finite DoS. By contrast, any superconducting order at finite doping maximally gaps the Fermi surface. Therefore, at finite chemical doping superconducting orders are energetically superior to the excitonic orders, and they can be realized at sufficiently low temperature even in the presence of repulsive electronic interactions. Concomitantly, the particle-hole orders are pushed to the higher temperature and interaction regime, see Figs. 4 and 5.
Even though we gain valuable insights into the organization of various BSPs in the global phase diagram of strongly interacting spin-3/2 fermions from the competition between energy and entropy inside the ordered phases (guided by emergent topology of reconstructed band structure), the phase diagrams shown in Figs. 3, 4 and 5 are obtained from an unbiased RG analysis, which systematically accounts for quantum fluctuations beyond the saddle point or mean-field approximation. Next we highlight the key aspects of the RG analysis.

IIB2. Methodology: Renormalization Group
The RG analysis we pursue in this work is controlled by a "small" parameter , measuring the deviation from the lower critical two spatial dimensions (d    The leading order RG analysis can be summarized in terms of the following set of coupled flow equations where is the logarithm of the RG scale, t = 2mT /Λ 2 and µ = 2mµ/Λ 2 are respectively the dimensionless temperature and chemical potential. For brevity we takeμ → µ in the above flow equations. Here {g i } represents the set of dimensionless four-fermion interactions with g i ∼ Λ λ i and H jk (α, t, µ) are functions of the mass anisotropy parameter α, t and µ. The RG flow equations for g i s are obtained by systematically accounting for quantum corrections to the quadratic order in the g i s. The relevant Feynman diagrams are shown in Fig. 14. A more detailed discussion of the RG analysis is presented in Sec. VI and Appendix E. Some salient features of the RG analysis are the followings.
2. Any weak local four-fermion interaction is an irrelevant perturbation and all orderings (realized when g i ( * ) → ∞) take place at finite coupling g i ∼ through quantum phase transitions (QPTs). Such QPTs are controlled by quantum critical points (QCPs) and all transitions are continuous in nature. The universality class of the transition is determined by two critical exponents, given by and for the physically relevant situation = 1.
Using the RG analysis we arrive at various cuts of the global phase diagram of interacting spin-3/2 fermions at (1) zero chemical doping [see . The universality class of the QPT leaves its signature on the scaling of the transition temperature (t c ). Note that t c ∼ |δ i | νz [81,82], where Hence, t c ∼ |δ| 2 for ν = 1 and z = 2, obtained from the leading order expansion, after setting = 1, irrespective of the choice of the coupling constant and the resulting BSP (see Fig. 17). We discuss this issue in detail in Sec. VI C.

IIB3. Competing Orders & Selection Rule
The correspondence between a given interaction coupling and the resulting phases can be appreciated by formulating the whole theory in the basis of an eight component Nambu-doubled spinor Ψ Nam , introduced in Sec. III C. In this basis the Luttinger Hamiltonian h L (k) → η 3ĥL (k). Pauli matrices {η µ } operate on the Nambu or particle-hole index. Any four-fermion interaction takes the form g I (Ψ † NamÎ Ψ Nam ) 2 and an order parameter (∆ O ) couples to a fermion bilinear according to ∆ O (Ψ † NamÔ Ψ † Nam ), whereÎ andÔ are eight dimensional Hermitian matrices. We argue that when g I is sufficiently strong, it can support only two types of ordered phases, for which 6 either (1)Ô ≡Î or (2) {Ô,Î} = 0.  The A2u and T1u magnetic phases are shown in purple and magenta. Due to a large separation of the interaction strength g 4 required for any ordering at µ = 0 and |µ| > 0, we display the µ = 0 cut of the phase diagram from panel (d) in Fig. 7. The region at weaker interaction and higher temperature is occupied by correlated Luttinger metal, without any long-range ordering. In panels (a) and (b) 100 g i → g i and in (c) 10 g i → g i . This outcome can be appreciated in the following way.
When an interaction coupling g I diverges toward +∞ under coarse graining (indicating onset of a BSP), it provides positive scaling dimension to an order parameter field ∆ O only when one of the above two conditions is satisfied. We substantiate this argument by considering the relevant Feynman diagrams [see Fig. 15] in Sec. VI E. Even though we arrive at this "selection rule" among competing orders from a leading order RG calculation, such a simple argument relies on internal symmetries among competing orders (breaking different symmetries) and is expected to hold at the non-perturbative level. We now support this claim by focusing on a specific example.
Let us choose a particular four-fermion interaction (in the T 2g nematic channel) From the phase diagrams shown in Figs. 3(a) for zero and finite temperature and 4(left) for finite chemical doping, we find that when this coupling constant is sufficiently strong, it supports two distinct BSPs.
Moreover we realize that the T 2g nematic order and the s-wave pairing together constitute an O(5) vector, of five mutually anticommuting matrices, reflecting the enlarged internal symmetry between these two orders. Following the same spirit we arrive at the following observations. A more detailed discussion supporting these scenarios is presented in Sec. VI E. Therefore, combining the energyentropy competition (obtained from the reconstructed band topology within a mean-field approximation) and an unbiased (controlled by a small parameter ) RG analysis with the selection rule we gain valuable insights into the nature of broken symmetry phases, competing orders and quantum critical phenomena in the global phase diagram of strongly interacting spin-3/2 fermions.

IIB4. Anisotropic Luttinger (semi)metal
So far, we centered our focus on the isotropic system [realized for α = π 4 in Eq. (1)]. Note that for α = π 4 the system enjoys an enlarged (but artificial) SO(3) rotational symmetry. However, in a cubic environment α = π 4 in general. Nonetheless, all the central results we quoted in the last three subsections hold (at least qualitatively) for any arbitrary value of α: 0 ≤ α ≤ π 2 . The discussion on the role of the mass anisotropy parameter α on the global phase diagram of interacting spin-3/2 fermions is rather technical, which we address in depth in Secs. V A, VI A, VI D. We here only quote some key results, which nicely corroborate with the rest of the discussion from this section.
1. We identify the mass anisotropy parameter as a valuable non-thermal tuning parameter, and for suitable choices of this parameter one can find (a) T 2g nematic and A 2u magnetic order respectively for strong enough g 1 and g 3 couplings (when α → π 2 ), see Figs. 6(a) and 6(c), (b) E g nematic and T 1u magnetic orders for strong enough g 2 and g 4 (when α → 0), as shown in Figs. 6(b) and 7 at zero temperature and chemical doping, respectively. These outcomes are in agreement with selection rule (1).  (2), as we argued previously for an isotropic Luttinger system.

At
We now proceed to a detailed discussion on each component of our work, starting from the noninteracting Luttinger model.

III. LUTTINGER MODEL
We begin the discussion with the Luttinger model describing a bi-quadratic touching of Kramers degenerate valence and conduction bands at an isolated point (here chosen to be the Γ = (0, 0, 0) point, for convenience) in the Brillouin zone. In this section, we first present the low-energy Hamiltonian and discuss its symmetry properties (Sec. III A). Subsequently, we introduce the corresponding imaginary time (τ ) or Euclidean action and the notion of the renormalization group (RG) scaling (Sec. III B). Finally, we define an eight-component Nambu-doubled basis that allows us to capture all, including both particle-hole or excitonic and particle-particle or superconducting, orders within a unified framework (Sec. III C).

A. Hamiltonian and Symmetries
The Hamiltonian operator describing a bi-quadratic touching of Kramers degenerate valence and conduction bands in three dimensions is given by [3,4] where Γ 0 is the four-dimensional identity matrix. Chemical potential µ and momenta k are measured from the band touching point. Here,d(k) is a five-dimensional unit vector that transforms under the l = 2 representation under the orbital SO(3) rotations. Henced(k) is constructed from the d-wave form factors or spherical harmonics Y m l=2 (θ, φ), as shown in Appendix C. The corresponding four-component spinor basis is given by where c k,j is the fermionic annihilation operator with momenta k and spin projection j = ±3/2 and ±1/2. The five mutually anti-commuting Γ matrices are defined as Two sets of two dimensional Pauli matrices {κ ν } and {σ ν } respectively operate on the sign (sgn[j]) and magnitude (|j| ∈ {1/2, 3/2}) of the spin projections, where ν = 0, 1, 2, 3. To close the Clifford algebra of all fourdimensional Hermitian matrices we also define ten commutators according to Γ jk = [Γ j , Γ k ] /(2i), with j > k and j, k = 1, · · · , 5. All sixteen four-dimensional matrices can be expressed in terms of the products of spin-3/2 matrices (J), as also shown in Appendix C. The energy spectra for Luttinger fermions are given by reflecting the quadratic band touching for µ = 0, which is protected by the cubic symmetry. For brevity we drop the explicit dependence of {d j } onk. Notice that the independence of E s (k) on s manifests the Kramers degeneracy of the valence and conduction bands, ensured by (1) the time-reversal (T ) and (2) the parity or inversion (P) symmetries. Specifically, under the reversal of time, k → −k and Ψ k → Γ 1 Γ 3 Ψ −k and hence T = Γ 1 Γ 3 K, where K is the complex conjugation, yielding T 2 = −1 (reflecting Kramers degeneracy of bands). Under the inversion P : k → −k and Ψ k → Ψ −k .
The "average" mass m and the mass anisotropy parameter α are respectively given by [28] Note that {d j } for j = 1, 2, 3 and j = 4, 5 belong to the T 2g (three component) and E g (two component) representations of the cubic or octahedral (O h ) point group, and m 1 and m 2 are effective masses in these two orbitals, respectively. The mass anisotropy parameter α allows us to smoothly interpolate between (1) the m 1 → ∞ limit when the dispersion of the T 2g orbital becomes flat, yielding α → 0 and (2) m 2 → ∞ when the E g orbital becomes non-dispersive, leading to α → π 2 . For α = π 4 or m 1 = m 2 , the Luttinger model enjoys an enlarged spherical symmetry. Any α = π 4 captures a quadrupolar distortion in the system (still preserving the cubic symmetry). In what follows, we treat α as a non-thermal tuning parameter to explore the territory of interacting Luttinger fermions.
The connection between the spin projections (j = ±3/2 and ±1/2) and the bands can be appreciated most economically by taking k = (0, 0, k). For such a specific choice of momentum axis, the Luttinger Hamiltonian takes a diagonal form, given bŷ From the above expression we can immediately infer that the pseudospin projections on the valence and conduction bands are respectively |j| = 3/2 and 1/2.

B. Lagrangian and Scaling
The imaginary time (τ ) Euclidean action corresponding to the non-interacting Luttinger model is given by The action remains invariant under the following rescaling of space-(imaginary)time coordinates and the fermionic field where z is the dynamic scaling exponent, measuring the relative scaling between energy and momentum according to E(k) ∼ |k| z . For Luttinger fermions z = 2. The parameter is the logarithm of the RG scale. In what follows in Secs. V and VI, we use the above scaling ansatz while addressing the effects of electronic interactions in this system. Under the above rescaling of parameters, the temperature (T ) and chemical potential (µ) scale as Therefore, the scaling dimension of these two quantities is [T ] = [µ] = z (same as that of energy). Throughout, we use the natural unit, in which = k B = 1.

C. Nambu doubling
To facilitate the forthcoming discussion we here introduce an eight-component Nambu-doubled spinor basis (suitable to capture both exitonic and superconducting orders within a unified framework) according to where Ψ k is a four-component spinor, see Eq. (14). In the lower block of Ψ Nam we absorb the unitary part of the time-reversal operator T , ensuring that the eightcomponent Nambu spinor (Ψ Nam ) transforms the same way as the original four component spinor Ψ k under the SU (2) pseudospin rotation. In this basis the eightdimensional Luttinger Hamiltonian takes a simple form and the time-reversal operator becomes T Nam = η 0 Γ 1 Γ 3 K. The newly introduced set of Pauli matrices {η ν } operates on the Nambu or particle-hole indices, with ν = 0, 1, 2, 3. Therefore, by construction while the excitonic orders assume block-diagonal form, all superconducting orders are block-off-diagonal in the Nambu subspace. Note thatĥ Nam L (k) commutes with the number operatorN = η 3 Γ 0 .

IV. BROKEN SYMMETRY PHASES
Next we discuss possible BSPs in this system. We introduce various possible excitonic and superconducting orders in the Nambu basis (Ψ Nam ) in two subsequent sections. Finally, we discuss the reconstructed band structure and emergent topology inside the ordered phases.
A. Particle-hole or excitonic orders The effective single-particle Hamiltonian in the presence of all possible momentum-independent or local or intra-unit cell excitonic orders is given by wherê The ordered phases can be classified according to their transformation under the cubic (O h ) point group symmetry. Regular fermionic density (∆ 0 ) does not break any symmetry (hence does not correspond to any ordering) and transforms under the trivial A 1g representation. A three-component nematic order-parameter, constituted by ∆ 1 = ∆ 1 1 , ∆ 2 1 , ∆ 3 1 , transforms under the T 2g representation. By contrast, a two-component nematic order transforming under the E g representation is captured by Both of them break only the cubic symmetry, but preserve time-reversal and inversion symmetries. The ordered phase represents either a time-reversal invariant insulator or a Dirac semimetal, about which more in a moment [see Sec. IV C]. Since five Γ matrices transform as components of a rank-2 tensor under SO(3) rotations, the two nematic phases represent quadrupolar orders, see Appendix C.
All ordered phases shown in the second line of Eq. (25) break time-reversal symmetry and represent different magnetic orders. For example, ∆ 3 corresponds to an octupolar order (since Γ 45 ∼ J x J y J z ), transforming under the singlet A 2u representation. In a pyrochlore lattice of 227 iridates such an ordered phase represents the "allin all-out" arrangement of electronic spin between two adjacent corner-shared tetrahedra [8,20]. By contrast, "two-in two-out" or "spin-ice" magnetic orderings on a pyrochlore lattice are represented by a three-component vector ∆ 4 = ∆ 1 4 , ∆ 2 4 , ∆ 3 4 (accounting for six possible two-in two-out arrangements in a single tetrahedron).
j such an ordered phase contains a linear superposition of dipolar and octupolar moments, and transforms under the T 1u representation [28]. Any other magnetic ordering can be represented by a six component vector ∆ jk 5 with j = 1, 2, 3 and k = 4, 5. No physical realization of such multi-component magnetic ordering in any material is currently known, and we do not delve into the discussion on such ordering for the rest of the paper. B. Particle-particle or superconducting orders The effective single particle Hamiltonian in the presence of all possible momentum-independent or local or intra-unit cell superconducting orders reads [49,[51][52][53] and φ is the global U (1) superconducting phase. Any pairing proportional to η 1 (η 2 ) preserves (breaks) timereversal symmetry (recall that the time-reversal operator in the Nambu basis is T Nam = η 0 Γ 1 Γ 3 K). Here, ∆ p A1g is the amplitude of the s-wave pairing, transforming under the A 1g representation. The s-wave pairing breaks only the global U (1) symmetry, but preserves the cubic symmetry. On the other hand, ∆ p,j T2g captures the amplitude of three d-wave pairings (for j = 1, 2, 3) transforming under the T 2g representation, and ∆ p,j Eg for j = 4, 5 represents the amplitude of two d-wave pairings belonging to the E g representation. Notice {Γ j , j = 1, · · · , 5} can be expressed in terms of the product of two spin-3/2 matrices, and all five d-wave pairings break the cubic symmetry, while introducing a lattice distortion or electronic nematicity in the system. Hence, they stand as representatives of quadrupolar nematic superconductors.

C. Reconstructed band structure and emergent topology
Next we consider the reconstructed band structure inside different BSPs which provides valuable information regarding the emergent topology inside ordered phases. The onset of any ordering discussed in the previous sections destabilizes the bi-quadratic touching and gives rise to either gapped or gapless quasiparticles (see below). Furthermore, this exercise will allow us to appreciate the energy-entropy competition among different orderings [see Sec. IV D], which ultimately plays a decisive role in the organization of various phases in the global phase diagram of interacting Luttinger fermions.
1. T 2g nematicity: The three component orderparameter for the T 2g nematic phase gives birth to gap- less quasiparticles for the following four configurations These four phase lockings are respectively shown as blue, red, green and black points in Fig. 8(a). The gapless phase corresponds to a topological Dirac semimetal (since nematicity preserves the Kramers degeneracy of valence and conduction bands), similar to the ones recently found in Cd 3 As 2 [83] and Na 3 Bi [84]. The DoS in a Dirac semimetal vanishes as (E) ∼ |E| 2 . The Dirac points are located along the body diagonals (the C 3v axes) of a cubic system and respectively placed at where k 0 = [2m 1 ∆ 1 /3] 1/2 , as shown in Fig. 8(b). For any other phase locking within the T 2g sector the system becomes an insulator. The spectral gap in the insulating phase is anisotropic and it is energetically superior over the gapless Dirac semimetal phase.
2. E g nematicity: The two component E g nematic order is most conveniently described in terms of the following parametrization where φ Eg is the internal angle in the order-parameter space. Only for the quasi-particle spectra are gapless, as shown in Fig. 9(a), and the ordered phase represents a topological Dirac semimetal. Specifically, for φ Eg = 0, 2π/3, 4π/3, the Dirac points are respectively located on k z , k x and k y axes (the C 4v axes), see Fig. 9(b), and the separation of two Dirac points is given by 2k 0 , where k 0 = 2m 2 ∆ 2 / √ 2 1/2 . For any other phase locking within the E g sector, the system becomes an insulator. 3. A 2u magnet: In the presence of an octupolar A 2u ordering, the two-fold degeneracy of the valence and conduction band gets lifted and a pair of Kramers nondegenerate bands touch each other at the following eight points in the Brillouin zone [see Fig. 10 where k 0 = 2m 1 ∆ 3 /3. They represent simple Weyl points, which act as source (4 of them) and sink (4 of them) of Abelian Berry curvature of unit strength. However, due to an octupolar arrangement of the Weyl nodes, the net Berry curvature through any high-symmetry plane is precisely zero and this phase does not support any anomalous Hall effect. The DoS at low energies then scales as (E) ∼ |E| 2 [8,20,28] at (0, 0, ±k 0 ) where k 0 = √ 2m 2 ∆ 4 and a nodal-loop is found in the k x − k y plane, as shown in Fig. 10(b). Similarly, for j = 1 and 2 the Weyl nodes are separated along the k x and k y axes, and the nodal-loops are respectively found in the k y − k z and k x − k z planes. Due to the presence of two Weyl nodes, each configuration of twoin two-out magnetic order supports a finite anomalous Hall effect in the plane perpendicular to the separation of the Weyl nodes. However, any triplet magnetic order, represented by ∆ 4 = |∆ 4 |(±1, ±1, ±1)/ √ 3, gets rid of the nodal loop and supports only two Weyl nodes along one of the body-diagonals (C 3v axes). Hence, triplet T 1u magnetic orders are energetically favored over their uniaxial counterparts [28]. 7

5.
A 1g or s-wave pairing: Notice that the matrix operator representing an s-wave pairing fully anti-commutes with the Luttinger Hamiltonian (for any value of α) and thus corresponds to a mass for Luttinger fermions. The quasiparticle spectra inside the paired state is fully gapped, but the phase is topologically trivial.
6. T 2g pairing: Three d-wave pairings, proportional to Γ 1 , Γ 2 and Γ 3 matrices, belong to the T 2g representation 7 The low energy DoS in the presence of a nodal loop and two point nodes (due to a uniaxial T 1u order) is dominated by the former and scales as (E) ∼ |E|, while in a triplet T 1u state the DoS scales as (E) ∼ |E| 2 (due to the point nodes). Hence, formation of the triplet ordering causes power-law suppression of the DoS and increases the condensation energy gain.
Pairing IREP. Equations for nodal loops Symmetry We display the symmetry of each d-wave pairing in the proximity to the Fermi surface (realized on the conduction or valence band) in the last column. Note that two nodal loops for dxy, dxz, dyz and d x 2 −y 2 parings can be rotated into each other, while those in the presence of d 3z 2 −r 2 pairing are disconnected from the remaining ones, see Fig. 11. For the sake of simplicity we here assume m1 = m2 = m, for which the nodal loops are circular in shape. For m1 = m2, the nodal loops become elliptic. Here, ∆ is the amplitude of d-wave pairings. and respectively possess the symmetry of d yz , d xz and d xy pairings. Each component supports two nodal loops in the ordered phase, as shown in the first three rows of Table I [49,53]. Two nodal loops for the Γ 3 or d xy pairing are shown in Fig. 11(a). The two nodal loops for Γ 1 or d yz and Γ 2 or d xz pairings can respectively be obtained by rotating the ones shown for d xy pairing by π 2 , with respect to the k y and k x axes.
7. E g pairing: E g pairings proportional to Γ 4 and Γ 5 matrices respectively possess the symmetry of d x 2 −y 2 and d 3z 2 −r 2 pairings and each of them supports two nodal loops, as shown in the last two rows of Table I [49,53]. Note that two nodal loops for the d x 2 −y 2 pairing can be realized by rotating the ones for the d xy pairing by π 4 about the k z axis. However, two nodal loops for the d 3z 2 −r 2 pairing, shown in Fig. 11(b), cannot be rotated into the ones for d x 2 −y 2 pairing. Therefore, despite the fact that the d 3z 2 −r 2 and d x 2 −y 2 pairings belong to the same E g representation, they are not energetically degenerate [53,85]. Since the radius of the nodal loops for the d 3z 2 −r 2 pairing is the smallest, this paired state is the energetically most favorable among five d-wave pairings.  Table I.

D. Energy and Entropy Inside Ordered Phases
From the computation of the reconstructed band structure we can gain insight into the condensation energy (∆ F ) and entropy (∆ S ) inside the ordered phases. While the stiffness of the spectral gap measures the gain of condensation energy, the scaling of the DoS at low-energies (due to gapless quasiparticles) measures the entropy. Recall that the s-wave pairing yields fully gapped spectra (isotropic), while the nematic orders produce either an anisotropic gap or gapless quasiparticles. Hence, the former ordering is associated with higher (lower) gain in condensation energy (entropy). On the other hand, the DoS vanishes as (E) ∼ |E| and |E| 2 respectively in the presence of a nodal-loop and Dirac or Weyl points. We found that the A 2u magnetic order gives birth to eight Weyl nodes, while only two Weyl nodes can be found inside the triplet T 1u magnetic order. By contrast, all five d-wave pairings are accompanied by two nodal loops (see Table I). Therefore, we can organize these ordered phases according to their contribution to (a) condensation energy and (b) entropy gain, as shown in Fig. 2. The LSM, on the other hand, accommodates the largest amount of gapless fermionic excitations near zero energy, where the DoS vanishes as (E) ∼ √ E. Hence, the LSM is endowed with largest entropy. 9 We now present a simple prescription to estimate (at Therefore, for the s-wave pairing, two nematic orders, two magnetic orders and five d-wave pairings N M anti = 5, 4, 2, 1, while N M comm = 0, 1, 3, 4, respectively. Then for any ordered phase This correspondence can be anchored from a simple example. Let us choose s-wave pairing, for which N M anti = 5, N M comm = 0 and chemical potential, for which N M anti = 0, N M comm = 5, as two perturbations in a LSM. While the s-wave pairing yields an isotropic gap, a finite chemical doping creates a Fermi surface (producing a constant DoS). Consequently, the s-wave pairing (chemical doping) is accompanied by larger gain of condensation energy (entropy). Therefore, the analysis of reconstructed band structure and emergent topology inside BSPs allows us to organize them according to the gain of condensation energy and entropy. The RG analysis at finite temperature (in the regime where the dimensionless temperature t = 2mT /Λ 2 1) captures such energy-entropy competition, which we discuss in Sec. VI C, see also Fig. 3.
The hierarchy of the energy and entropy gains inside the ordered phases changes when the chemical potential is placed away from the band touching point (i.e., orders (producing Weyl nodes), or the stiffness of the spectral gap, such as between T 2g and Eg nematic orders (producing anisotropic gaps). A more microscopic analysis is needed to resolve these situations. µ = 0). Since any pairing operator anti-commutes with the number operator (N = η 3 Γ 0 ), superconducting orders maximally gap (either fully by the s-wave pairing or partially by the individual d-wave pairings) the Fermi surface. By contrast, any excitonic order always gives birth to a Fermi surface, according to the Luttinger theorem [78]. Hence, at finite doping all superconductors are energetically superior over the particle-hole orders, while excitonic orders are accompanied by larger entropy (due to presence of a Fermi surface). The energy-entropy competition at finite-µ is also captured by the RG analysis, discussed in Sec. VI D, leading to the phase diagrams shown in Figs. 4, 5 and 6.

V. ELECTRON-ELECTRON INTERACTIONS
Next we proceed to demonstrate the onset of various BSPs, discussed in the previous section, triggered by repulsive (at the bare level) electron-electron interactions. As mentioned earlier we will focus only on the local or short-range part of Coulomb interaction and neglect its long-range tail. For the sake of concreteness we assume that the local interactions are density-density in nature. Any generic local density-density interaction (such as the ones appearing in an extended Hubbard model, for example) can be captured by six quartic terms and the corresponding interacting Hamiltonian reads In this notation λ j > 0 corresponds to repulsive interaction. However, all four-fermion interactions are not linearly independent due to the existence of Fierz identity among sixteen four-dimensional Hermitian matrices, closing a U (4) Clifford algebra [see Appendix D] [86][87][88]. It turns out that any generic local interaction can be expressed in terms of only three quartic terms, and we conveniently (without any loss of generality) choose them to be λ 0 , λ 1 and λ 2 . Following the Fierz relations we can express local quartic terms proportional to λ 3 , λ 4 , λ 5 as linear combinations of above three, see Eq. (D4). Whenever we generate four-fermion interactions proportional to λ 3,4,5 during the coarse-graining (discussed in Sec. VI), they can immediately be expressed in terms of λ 0,1,2 , and the interacting model defined in terms of λ 0,1,2 [see Eq. (36) below] always remains closed under the RG procedure to any order in the perturbation theory. The imaginary time Euclidean action for the interact-ing system is given by where Ψ ≡ Ψ(τ, r) and Ψ † ≡ Ψ † (τ, r). Under the rescaling of space-time(imaginary) coordinates and the fermionic fields, see Eq. (20), the local four-fermion interaction scales as Hence the scaling dimension of quartic couplings is [λ j ] = z −d. For z = 2 and d = 3, [λ j ] = −1, and any weak local interaction is an irrelevant perturbation and leaves the Luttinger fermions unaffected. Therefore, any ordering sets in at an intermediate strength of coupling through QPT. In Sec. VI we demonstrate appearances of various BSPs using a RG analysis, controlled via an -expansion, where = d − 2, about the lower-critical two spatial dimension of this theory. Within the framework of the -expansion, such QPTs take place at a critical interaction strength λ * j ∼ , and in three spatial dimensions (d = 3) = 1. Before proceeding to the RG analysis, we seek to gain some insight into the propensity toward various orderings by computing the corresponding mean-field susceptibility for a wide range of the mass anisotropy parameter (α). Readers interested in the RG analysis may skip the following discussion and directly go to Sec. VI.

A. Mean-field susceptibility
To gain insights into the propensity toward the formation of various orderings, we first compute the bare mean-field susceptibility (χ M ) of all possible symmetry allowed fermionic bilinears Ψ † Nam M Ψ Nam , where M is an eight dimensional Hermitian matrix (see Sec. IV). For simplicity we set µ = 0. For zero external momentum and frequency this quantity is given by (38) The relevant Feynman diagram is shown in Fig. 12 and the "−" sign arises from the fermion bubble. Here G k (iω n ) is the fermionic Green's function in the Nambu doubled basis. The factor of 1/2 takes care of the artificial Nambu doubling. Results are displayed in Fig. 13. Next we discuss the scaling of χ M in different channels for a few specific values of the mass anisotropy parameter.  For α = π 4 , the effective masses for the T 2g and E g orbitals are equal (i.e. m 1 = m 2 ) and the system enjoys an enlarged spherical symmetry. Since each one of the five Γ-matrices (representing two nematic orders) anticommutes with four matrices and commutes with one matrix appearing in the Luttinger Hamiltonian, two nematic orders belonging to the T 2g (red curve) and E g (orange curve) representations possess equal susceptibility. On the other hand, all ten commutators (representing various magnetic orders) anti-commute with three and commute with two matrices appearing in this model. Hence, magnetic orders in the A 2u (purple curve) and T 1u (magenta curve) channels also possess equal susceptibility. Two copies of the d-wave pairing, transforming under the T 2g (dark green curve) and E g (dark yellow curve) representations, have degenerate susceptibilities, as all five d-wave pairing matrices commute with four matrices and anti-commute with only one matrix appearing in the Luttinger model. As the s-wave pairing fully anti-commutes with the Luttinger Hamiltonian, it always possesses the largest susceptibility (blue curve) for any α.
Susceptibilities for different BSPs (characterized by the fermion bilinear Ψ † Nam M Ψ Nam ) ∼ N M anti , the number of matrices inĥ Nam L (k) anti-commuting with M . This simple correspondence is operative irrespective of the choice of α. Note that in an isotropic system N M anti = 5, 4, 3 and 1 for the s-wave pairing, nematicity, magnetic orders and d-wave pairings, respectively, hence Computation of the bare susceptibility suggests a strong propensity toward the formation of s-wave pairing and two nematic orders in the world of interacting spin-3/2 fermions with isotropic dispersion. The magnetic orders and d-wave pairings are expected to be suppressed near α = π 4 , at least at zero temperature [see Figs. 3 and 16]. We also note that the gain in free-energy [see Sec. IV D] and the mean-field susceptibility follow the same hierarchy ∆ F , χ M ∼ N M anti .

Anisotropic Luttinger Semimetal near
When the effective mass in the T 2g orbital becomes sufficiently large, the Luttinger model simplifies to This Hamiltonian possesses an emergent SU (2) ⊗ U (1) chiral symmetry, where {Γ 45 Γ j } with j = 1, 2, 3 are the three generators of an SU (2) rotation, whereas a U (1) rotation is generated byN = η 3 Γ 0 , the number operator. In this limit the Hamiltonian is similar to the one for spinless fermions in Bernal-stacked bilayer graphene, which altogether supports six masses, given by 10 (40) Note that three components of the T 2g nematicity break the continuous SU (2) chiral symmetry, while the s-wave pairing breaks the global U (1) symmetry. On the other hand, the A 2u magnet transforms as a scalar under the chiral rotation and breaks only time-reversal-symmetry. These three mass orders possess the largest and equal susceptibilities as α → π 2 , see Fig. 13. Therefore, repulsive interactions favor two excitonic masses for zero [see  Also note that each member of the following vector anti-commutes and commutes with one matrix appearing in lim α→ π 2ĥ Nam L (k). Hence, E g nematicity and d-wave pairing have identical susceptibilities as α → π 2 , but As a result such an anisotropic system can accommodate an E g d-wave pairing at finite chemical doping, see Fig. 6(c). However, T 1u magnet and T 2g d-wave superconductor have exactly zero susceptibility as they fully commute with lim α→ π 2ĥ Nam L (k). Hence, onset of these two orders is unlikely when α ≈ π 2 .
3. Anisotropic Luttinger Semimetal near α = 0 Finally, we compare the susceptibility for various orders when the mass of the E g orbital becomes sufficiently large. The Luttinger Hamiltonian then takes the form 11 which possesses an emergent U (1)⊗U (1) symmetry, generated by η 0 Γ 45 and η 3 Γ 0 . The mass orders in this limiting scenario constitute the following vector and consequently the E g nematicity and s-wave pairing acquire an identical and the largest susceptibility as α → 0, see Fig. 13. Hence, a competition between these two ordered phases can be anticipated near α = 0 [see Figs. 6(b) and 16(b)]. Any order parameter from the following vector anti-commutes with two matrices and commutes with one matrix appearing in Eq. (42) and they also possess degenerate susceptibilities, but As a result, a competition between T 1u magnet and T 2g nematicity can also be observed around α = 0, see Fig. 16(d). On the other hand, the T 2g d-wave pairing matrices anti-commute with one matrix and commute with two matrices appearing in Eq. (42) and its susceptibility is smaller than the orders appearing in M 0 and M 0 . Nonetheless, when assisted by finite chemical doping, the T 2g d-wave pairing can be realized even for repulsive magnetic interaction in the T 1u channel, as shown in Fig. 6(d). Finally, we note that A 2u magnet and E g d-wave pairing fully commute with lim α→0ĥ (k) and possess zero susceptibility. Hence, onset of these two orders around α = 0 is unlikely. Note that various matrices appearing in M π 2 , M π

VI. RENOMRALIZATION GROUP ANALYSIS
After gaining insight into the propensity toward various orderings in the Luttinger system, next we seek to investigate the onset of different BSPs and the competition among them within the framework of an unbiased RG analysis. This will allow us to go beyond the mean-field analysis, presented in the last section, and systematically incorporate fluctuations. In what follows we here restrict ourselves to the leading order in the expansion, where = d − 2, and account for corrections to the bare interaction vertices (λ j s) to quadratic order in the coupling constants. The relevant Feynman diagrams are shown in Fig 14. After performing summation over fermionic Matsubara frequencies ω n = (2n + 1)πT with −∞ ≤ n ≤ ∞, we integrate out a thin Wilsonian shell Λe − < |k| < Λ to arrive at the following RG flow equations for i = 0, 1, 2, where β X ≡ dX/d , in terms of dimensionless quantities defined as The prime symbol in the summation indicates that j > k.    (17)]. Note for α = π 4 the Luttinger model possesses an emergent spherical symmetry, while for α → 0 ( π 2 ) the band dispersion in the Eg(T2g) orbitals becomes almost flat. For α = π 4 the coupled RG flow equations support only one QCP, denoted by QCP 1 π 4 , while for α → π 2 and α → 0 we find two QCPs, respectively denoted by QCP j π 2 and QCP j 0 for j = 1, 2, see also Ref. [43]. For these two limiting cases the flow equations also support one BCP, respectively identified as BCP π 2 and BCP0. Existence of such BCP in the presence of two QCPs is necessary to ensure the continuity of the RG flow trajectories and separate the basins of attraction of two QCPs. All coupling constants at the fixed points are measured in units of , where = d − 2 measures the deviation from the lower critical two spatial dimensions, where all local quartic interactions are marginal. Note that QCP 1 π 4 , QCP 1 π 2 , QCP 1 0 represent the same QCP, only its location shifts as we tune α. By contrast, QCP 2 π 2 and QCP 2 0 are solely introduced by mass anisotropy and bear no analog to the fixed points found around α = π/4. Besides the above QCPs and the BCPs, there always exists a trivial Gaussian fixed point, representing the non-interacting Luttinger semimetal, at (g * 0 , g * 1 , g * 2 ) = (0, 0, 0), endowed with three stable directions. See Appendix F for details.
where t(0) and µ(0) are the bare values. We supply these solutions to Eq. (45) to find the phase diagram of interacting Luttinger fermions using the following prescription. 12 While the divergence of at least one of the quartic couplings (i.e., g i → +∞) under coarse graining indicates the onset of a BSP, to unambiguously determine the pattern of symmetry breaking we also account for the leading order RG flow of all source terms, associated with different BSPs, discussed in Sec. IV. The effective action in the presence of all symmetry allowed fermionic bilinears [see Eqs. (24)- (27)] is given by with Ψ † Nam ≡ Ψ † Nam (τ, x) and Ψ Nam ≡ Ψ Nam (τ, x). Relevant Feynman diagrams are shown in Fig. 15. The re- 12 We should note that Feynman diagram (b) in Fig. 14 provides interaction driven corrections (linear in g i ) to µ. However, to maintain the order by order correction to quartic interactions, we neglect such corrections in Eq. (45) within the framework of the leading order -expansion. sulting RG flow equations take the following schematic form See Appendix E for explicit form of these flow equations. The quantities appearing on the right hand side of each equation, represent the scaling dimension of the corresponding order-parameter. We simultaneously run the flow of the quartic couplings (g j s) and the source terms (∆ j s). When at least one of the quartic couplings diverges and flows toward → +∞ (thus indicating onset of a BSP), we identify the source term (say ∆ j ) that diverges toward → +∞ fastest (assuming a possible scenario when more than one source term diverge toward +∞). The BSP is then characterized by the order-parameter ∆ j = 0. We use this strategy to determine various cuts of the global phase diagram of spin-3/2 Luttinger fermions, displayed in Fig. 16 (for t = µ = 0), Fig. 3 (for µ = 0, but finite-t) and Figs. 4, 5 and 6 (for finite-t and finite-µ). Next we discuss these cases in three subsequent sections.

A. Quantum criticality in Luttinger semimetal
We first discuss the effects of electronic interactions on a LSM (i.e. when the chemical potential is pinned at the band touching point) at zero temperature. The RG flow equations for µ = 0 and t = 0 can be derived by taking the limit µ → 0 and then t → 0 in Eq. (45), suggesting that weak interactions are irrelevant perturbations and any ordering takes place at finite coupling g i ∼ through a QPT. Next we discuss the following three cases separately (i) isotropic Luttinger system (α = π 4 ), large mass for (ii) the T 2g orbital (α → π 2 ) and (iii) the E g orbital (α → 0). Such systematic analysis will allow us to anchor our anticipations regarding the nature of BSPs from the mean-field susceptibility, discussed in Sec. V A. Note at t = µ = 0, the system is devoid of any natural infrared cutoff as t * , µ * → ∞. Hence, we run the flows of quartic couplings up to an RG time * → ∞ to determine the stability of LSM and the flows of the source terms to pin the pattern of symmetry breaking.

1.
Isotropic Luttinger semimetal (α = π 4 ) For α = π 4 the coupled RG flow equations support only one QCP, reported in Table II and identified as QCP 1 π 4 , besides the trivial (and fully stable) Gaussian fixed point at g * 0 = g * 1 = g * 2 = 0 (representing the stable LSM). This QCP controls QPTs from LSM to various BSPs (depending on the interaction channel). To gain insight into the nature of the candidate competing BSPs, we compute the scaling dimensions for all fermion bilinears at this QCP. The results are summarized in Table III. Note that the s-wave pairing has the largest scaling dimension at this QCP, while two nematic orders possess degenerate but second largest (and positive) scaling dimensions. But, the rest of the fermion bilinears possess negative scaling dimensions. Notice that the scaling dimensions for different orders at this QCP follow the same hierarchy as the mean-field susceptibilities, discussed in Sec. V A 1.
The phase diagrams in an interacting LSM are displayed in Fig. 3 for various interaction channels. Around α = π 4 strong repulsive nematic interactions (g 1 and g 2 ) favor s-wave pairing even in the absence of a Fermi surface (since µ = 0), see Figs. 3(a) and 3(b). On the other hand, strong magnetic interactions in the A 2u (g 3 ) and T 1u (g 4 ) channels respectively support T 2g [see Fig. 3

(d)]
and E g [see Fig. 3(c)] nematicities. However, we could not find any magnetic ordering or d-wave pairing in the very close vicinity to α = π 4 , at least when t = 0, see Fig. 16. Therefore, the computation of mean-field susceptibilities and scaling dimensions of fermion bilinears, in corroboration with our unbiased RG calculation, show that the strongly interacting isotropic LSM becomes unstable toward the formation of two nematic orders and s-wave pairing at the lowest temperature. ] and one bi-critical point [denoted by BCP π 2 ], see Table II. The BCP possesses two unstable directions. Notice QCP 1 π 2 is the same as QCP 1 π 4 , only shifted toward weaker coupling, which can be verified from the fact that the signs of the coupling constants and scaling dimensions for all fermion bilinears are identical at these two QCPs [see Tables II and III].
On the other hand, QCP 2 π 2 is new and bears no resemblance to any fixed points we found for α = π 4 . This QCP is induced by the mass anisotropy of Luttinger fermions. At this QCP the A 2u magnetic (E g nematic) order possesses the largest (second largest) scaling dimension, see Table III. Therefore, when repulsive interaction in the A 2u channel dominates among various finite range components of the Coulomb interaction, the Luttinger semimetal can display a competition between these two orders as we approach α = π 2 from an isotropic system, see Fig. 16(c) [consult also Sec. VI E].
Among three possible local pairings the E g d-wave superconductor possesses the largest (and positive) scaling dimension at this QCP. Hence, the emergence of this paired state can be anticipated when the LSM is doped away from the charge-neutrality point around α = π 2 , see  Finally, we approach the opposite limit, when the mass in the E g orbital becomes sufficiently large, i.e. m 2 m 1 or equivalently α → 0. In this regime the RG flow equations support two QCPs [denoted by QCP 1 0 and QCP 2 0 ] and a BCP [denoted by BCP 0 ], see , only shifted toward weaker coupling. By contrast, QCP 2 0 is induced by the mass anisotropy. At this QCP, the T 1u magnetic (T 2g nematic) order possesses the largest (second largest) scaling dimension, see Table III. Therefore, when repulsive interaction in the T 1u channel dominates we expect a strong competition between these two orderings, as one tunes toward α → 0 starting from an isotropic system, see Fig. 16(d).
Also note that among three local pairings, the d-wave one transforming under the T 2g representation possesses the largest (and positive) scaling dimension. Hence, around α = 0 we expect onset of this paired state when the chemical potential is placed away from the band touching point, see Fig. 6(d) and discussion in Sec. VI D.
Finally, we comment on the role of the BCPs possessing two unstable directions, see Table II. Note that a BCP can only be found when there exists two QCPs in the three dimensional coupling constant space (g 0 , g 1 , g 2 ). The existence of a BCP separates the basin of attraction of two QCPs and ensures continuity of the RG flow trajectories in coupling constant space.

B. Universality class and critical exponents
All QCPs, listed in Table II,  responding stability matrix [see Appendix F], defined as To the leading order in the -expansion the positive eigenvalue at all QCPs is exactly , which in turn determines the correlation length exponent (ν) according to For the physically relevant situation = 1 and we obtain ν = 1. The fact that ν is the same at all QCPs is, however, only an artifact of the leading order -expansion. Generically ν is expected to be distinct at different QCPs, once we account for higher order corrections in .
Since to the leading order in the -expansion, local interactions do not yield any correction to fermion self-energy, the dynamic scaling exponent (z) at all interacting QCPs is Together the correlation length and dynamic scaling exponents determine the universality class of all continuous QPTs from a LSM to various BSPs. In the next section, we will discuss the imprint of these two exponents on the scaling of the transition temperature (t c ) associated with the finite temperature order-disorder transitions.

C. Interacting Luttinger Semimetal at Finite Temperature
Even though our RG analysis in the previous section was performed at zero temperature, one can still find the imprint of various interacting QCPs at finite temperatures. The RG flow equations at finite temperature can be derived from Eq. (45) by taking the µ → 0 limit in H l jk (α, t, µ). Recall that temperature introduces a natural infrared cutoff for the flow equations t * [see Eq. (9)]. Physically such an infrared cutoff corresponds to a scenario when the renormalized temperature t( ) becomes comparable to the ultraviolet energy E Λ = Λ 2 /(2m) [see Fig. 1], beyond which the notion of quadratically dispersing fermions becomes moot and the flow equations from Eq. (45) lose their jurisdiction.
To capture the effects of electronic interactions in a LSM at finite-t, we run three quartic coupling constants up to a scale ≤ t * . Now depending on the bare strength of interactions, two situations arise respectively representing a disordered LSM (without any long-range ordering) or onset of a BSP at finite temperature. Hence, for a given strength of interaction g > g * , where g * is the requisite critical strength of interaction for a BSP at t = 0, we always find a temperature t c above (below) which the BSP disappears (appears). We identify t c as the critical or transition temperature. None of the coupling constants diverge for t > t c . All ordered phases can display true long-range order at finite t in three dimensions and t c corresponds to a genuine transition temperature. General scaling theory suggests that the transition temperature scales as [81,82] for δ 1, where δ = (g − g * ) /g * is the reduced distance from a QCP, located at g = g * (say). Hence, for interacting LSM, t c ∼ δ 2 for = 1 (i.e., the prediction from the leading order -expansion). The scaling of critical or transition temperature for various choices of coupling constants and resulting BSPs, and different choices of the mass anisotropy parameter α, are shown in Fig. 17, indicating a fairly good agreement with the field theoretic prediction t c ∼ δ 2 around all QCPs, reported in Table II.

Phase diagrams at finite temperature
Besides the scaling of the transition temperature, we also investigate the phase diagram of an interacting LSM at finite temperature, allowing us to demonstrate the competition between condensation energy gain and entropy. For concreteness we focus on the isotropic system (α = π 4 ), where this competition is most pronounced. As argued in Sec. IV D, the onset of s-wave pairing leads to the maximal gain in condensation energy, while the two nematic orders produce higher entropy in comparison to the former. Two specific cuts of the global phase diagram, Figs. 3(a) and 3(b), show that while s-wave pairing is realized at low temperature, nematicities set in at higher temperature as we increase the strength of nematic interactions (g 1 and g 2 ) in the system.
By contrast, when we tune the magnetic interactions (namely g 3 , g 4 ), an isotropic LSM becomes unstable in favor of two nematic orders at t = 0. Such an outcome can be substantiated from the simple picture of condensation energy gain, as A 2u and T 1u magnetic orders accommodate Weyl nodes (yielding more entropy), while nematicities produce anisotropic spectral gap (leading to higher gain of condensation energy), see Sec. IV D. As we tune the strength of the A 2u (T 1u ) magnetic interaction, E g (T 2g ) nematic order sets in at lower and A 2u (T 1u ) magnet at higher temperature, see Figs. 3(c) and 3(d).
The LSM is endowed with the largest entropy in the global phase diagram of interacting spin-3/2 fermions, since (E) ∼ √ E, in comparison to any BSP. Consequently, the requisite strength of interactions for any ordering increases with increasing temperature, irrespective of the nature of the BSP, see Fig. 3. Therefore, our RG analysis for an interacting LSM at finite temperature corroborates the energy-entropy competition picture and substantiates the following outcome: ordered phases providing larger condensation energy gain are found at lower temperature, while at higher temperature, phases with larger entropy are favored. This observation is also consistent with the notion of reconstructed band structure and emergent topology, discussed in Sec. IV C. Therefore, the phase diagram of an interacting LSM at finite temperature is guided by topological structure (gapped or nodal) of competing BSPs.
We close this section by answering the following question: Why do we find two magnetic orders in an isotropic LSM at finite temperature, since this system supports only one QCP, see Table II, where all the magnetic orders bear negative scaling dimension (see Table III)? Note that any QCP can only be accessed at t = 0, whereas finite-t introduces an infrared cutoff ( t * ) for the RG flow of the quartic couplings, and thus prohibits a direct access to any QCP. Hence, at finite-t, when magnetic interactions are sufficiently strong, the system can bypass the basin of attraction of QCP 1 π 4 and nucleate magnetic phases at moderately high temperature.

D. RG analysis in Luttinger metal
Finally we proceed to the RG analysis when the chemical potential (µ) is placed away from the band touching point. The chemical potential introduces yet another infrared cutoff µ * [see Eq. (10)], suggesting that FIG. 17: Scaling of dimensionless transition temperature (tc) with the reduced distance (δj = (g j − g * j )/g * j ) from the quantum critical point, located at g j = g * j for j = 1, 2, 3, 4, see Sec. VI C. The mass anisotropy parameter α and nature of the ordered phase for δj > 0 in each panel are quoted explicitly. Note that tc ∼ δ 2 j for δ 1, in agreement with field theoretic prediction. The QCPs controlling the Luttinger semimetal to ordered states quantum phase transitions at t = 0 are quoted in each panel. Red dots represent numerically obtained transition temperature and the blue lines correspond to least-square fit of tc ∼ δ 2 j .
the RG flow equations of the three quartic coupling constants should be stopped when the renoramlized chemical potential µ( ) reaches the scale of the band-width E Λ = Λ 2 /(2m). At finite temperature and chemical doping, two infrared scales compete and the smaller one * (say), given by * = min. µ * , t * determines the ultimate infrared cutoff for the RG flow of g j . Now depending on the bare strength of the coupling constants one of the following two situations arises (a) g j ( * ) > 1 or (b) g j ( * ) < 1.
While (a) indicates onset of a BSP, (b) represents a stable Luttinger metal. The resulting phase diagrams for various choices of chemical potential, temperature, coupling constants and the mass anisotropy parameter (α) are shown in Figs. 4, 5 and 6. We note that in an isotropic system and for strong enough nematic interactions an s-wave pairing can be realized even for zero chemical doping [see Figs. 3(a) and 3(b)]. With increasing doping the s-wave pairing occupies larger portion of the phase diagram, while two nematic phases get pushed toward stronger coupling, see Fig. 4. However, the d-wave pairings do not set in for zero chemical doping. Nonetheless, when the magnetic interactions in the A 2u and T 1u channels are strong, the presence of a Fermi surface is conducive to the nucleation of d-wave pairings, belonging to the E g and T 2g representations, respectively, see Fig. 5.
Now we focus on the anisotropic system. We chose the mass anisotropy parameter α such that at zero chemical doping the repulsive electronic interactions accommodate either nematic or magnetic orders, see Fig. 16. Specifically for α = 1.5 (close to π 2 ) the system enters into the and QCP 2 0 [see Sec. VI A and Table II]. At QPT 1 0 and QPT 1 π 2 the s-wave pairing possesses the largest scaling dimension among all possible local pairings [see Table III]. Hence, in the presence of finite chemical doping, repulsive interactions in the nematic channels become conducive to s-wave pairing, as shown in Figs. 6(a) and 6(b). On the other hand, at QPT 2 π 2 (QCP 2 0 ), the d-wave pairing belonging to the E g (T 2g ) representation possesses the largest scaling dimension. Therefore, repulsive interactions in the A 2u (g 3 ) and T 1u (g 4 ) magnetic channels become conducive to the nucleation of E g [for α = 1.5] and T 2g [for α = 0.1] dwave pairings, respectively, as shown in Figs. 6(c) and 6(d). We conclude that nematic and magnetic interactions among spin-3/2 Luttinger fermions are respectively conducive to the s-wave and d-wave pairings. Otherwise, at finite chemical doping the excitonic orderings set in only for stronger couplings. Such a generic feature is also consistent with the energy-entropy competition picture as the superconducting phases maximally gap the Fermi surface (yielding optimal gain of condensation energy), while exitonic orders are accompanied by a Fermi surface with constant DoS (producing more entropy).

E. Competing Orders and Selection Rule
So far we have presented an extensive analysis of the role of electronic interactions among spin-3/2 fermions. We showed multiple cuts of the global phase diagram at zero and finite temperature and chemical doping (see Figs. 3,4 and 5). In addition, we also addressed the imprint of the quadrupolar deformation (or the mass anisotropy parameter α) on a such phase diagram (see Figs. 6, 7 and 16). Altogether we unearth a rich confluence of competing orders in this system. In this context an important question arises quite naturally: Is there a selection rule among short-range interactions in different channels (such as nematic and magnetic) and various ordered states (such as s-and d-wave pairings, magnetic phases and nematic orders) for interacting spin-3/2 fermions? In this section we attempt to provide an affirmative (at least partially) answer to this question.
To this end it is convenient to express the quartic terms appearing in H int [see Eq. (35)] in the Nambu doubled basis. For concreteness, we focus on the relevant interaction channels, namely λ 1,2,3,4 . In the Nambu basis these four quartic terms take the form . 19: Schematic representation of SU(2) symmetry between the A2u magnetic order and d-wave pairing belonging to the Eg representation. Notations are the same as in Fig. 18.
The factor of 1/2 takes care of the artificial Nambu doubling. The above question can then be reformulated in the following way. When at least one of the coupling constants, say λ j , diverges toward +∞, what is the nature of the resulting BSP; or which one of the source terms, appearing in Eq. (48), simultaneously diverges toward +∞? This is our quest in this section. When the coupling constant λ of the quartic interaction λ Ψ † Nam M Int Nam Ψ Nam 2 diverges toward +∞ it nucleates a BSP, characterized by the order-parameter Ψ † Nam M or Ψ Nam , only if one of the following two conditions is satisfied (1) M or = M Int Nam or (2) M or , M Int Nam = 0. (56) If, on the other hand, M Int Nam and M or have more than one component, then selection rule (2) requires a slight modification, see footnote 6.
The above selection rule can be justified from the Feynman diagrams shown in Fig. 15. A source term diverges toward +∞ (indicating onset of a BSP) only when the net contribution from diagrams (b) and (c) is positive. Feynman diagram (b) gives non-zero and positive contribution only when condition (1) from Eq. (56) is satisfied (due to the Tr arising from the fermion bubble). Even though diagram (c) then yields a negative contribution, all togther they still produce a positive definite quantity, since the contribution from (b) dominates over (c), as the former one involves Tr. Hence, when condition (1) is satisfied, interaction λ can enhance the propensity toward the formation of a BSP with M or = M Int Nam . On the other hand, when condition (1) is not satisfied, only Feynman diagram (c) contributes [due to Tr involved in (b)]. The contribution from this diagram is positive only when condition (2) is satisfied. By contrast, when M or = M Int Nam,j and M or , M Int Nam,j = 0 the net contribution from (b) and (c) is negative. Interaction coupling λ then does not support such an ordered phase. All phase diagrams presented in this work support one of these two selection rules (discussed below).
The above selection rule can be phrased in a slightly different fashion. Two ordered phases respectively breaking O(M 1 ) and O(M 2 ) symmetries, 13 can reside next to each other (when we tune the strength of a particular interaction channel) if they constitute an O(M ) symmetric composite order parameter, where We invite the readers to verify the equivalence between Eqs. (56) and (57). Next, we discuss some prototypical examples to support our claims.  Fig. 19]. When the chemical potential is finite and we tune the strength of g 3 , the system accommodates a paired (magnetic) state at low (high) temperature, see Figs. 5 (top) and 6(c).
6. Finally, note that multiple copies of an O(3) vector can be formed by combining the components of T 1u magnetic order and T 2g d-wave pairing (an O(2) order-parameter), see Fig. 20. These two phases reside next to each other at finite chemical doping when we tune the quartic coupling g 4 , see Figs. 5 (bottom) and 6(d).
From the matrix representations of all quartic interactions [see Eq. (55)] and order parameters [see Eqs. (25) and (27)] the readers can convince themselves that the above examples are in agreement with our proposed selection rule from Eqs. (56) and (57). It is admitted that we arrive at the conclusion from a leading order RG analysis. However, the alternative version of the selection rule [see Eq. (57)] solely relies on the internal symmetry among competing orders. We therefore believe that our proposed selection rule is ultimately non-perturbative in nature, which can be tested in numerical experiments, for example.

VII. SUMMARY AND DISCUSSIONS
To summarize we have presented a comprehensive analysis on the role of electron-electron interactions in a three-dimensional Luttinger system, describing a biquadratic touching of Kramers degenerate valence and conduction bands (in the absence of chemical doping) of effective spin-3/2 fermions at an isolated point in the Brillouin zone. This model can succinctly capture the low-energy physics of HgTe [5], gray-Sn [6,7], 227 pyrochlore iridates [8][9][10][11][12] and half-Heuslers [13][14][15]. For concreteness, we focused only on the short-range components of the Coulomb interaction (such as the ones appearing in an extended Hubbard model), and neglected its long-range tail. Due to the vanishing density of states, namely (E) ∼ √ E, sufficiently weak, but generic local four-fermion interactions are irrelevant perturbations in this system. Consequently, any ordering sets in at an intermediate coupling through quantum phase transitions. We here address the instability of interacting Luttinger fermions at finite coupling within the framework of a renormalization group analysis, controlled by a "small " parameter , where = d − 2. Notice that in two spatial dimensions a bi-quadratic band touching (similar to the situation in Bernal-stacked bilayer graphene) yields a constant density of states, and local four-fermion interactions are marginal in d = 2. Hence, our renormalization group analysis is performed about the lower-critical two spatial dimensions. In this framework all quantum phase transitions from a disordered Luttinger (semi)metal to any ordered phase take place at g * j ∼ , where g j is the dimensionless coupling constant, and for three dimensions = 1. We here restrict ourselves to repulsive (at the bare level) electron-electron interactions and present multiple cuts of the global phase diagram at zero and finite temperature (see Figs. 3 and 16) and finite chemical doping (see Figs. 4, 5 and 6).
Using the renormalization group analysis, we show that in an isotropic system an s-wave pairing and two nematic orders are the prominent candidates for a broken symmetry phase at zero temperature and chemical doping. While the s-wave pairing only breaks the global U (1) symmetry and produces a uniform mass gap, two nematic orders, transforming under the T 2g and E g representations, produce lattice distortion along the C 3v and C 4v axes, respectively. Such ordered phases can describe either time-reversal symmetry preserving insulators or topological Dirac semimetals. However, a collection of strongly correlated gapless spin-3/2 fermions do not show any noticable propensity toward the nucleation of any magnetic order or d-wave pairings in an isotropic system at least when t = 0. This is so, because the magnetic orders (d-wave pairings) produce gapless Weyl nodes (nodal loops) around which the density of states vanishes as (E) ∼ |E| 2 (|E|), while the former three orders support gapped spectra. Hence, the magnetic orders and d-wave pairings are energetically inferior to s-wave pairing and electronic nematicities. However, with increasing temperature, one finds a smooth crossover from nematic to magnetic phases, as shown in Fig. 3. Hence, energetically superior orders are found at low temperatures, whereas at higher temperature broken symmetry phases possess larger entropy. This energy-entropy competition is discussed in Sec. IV D and summarized in Fig. 2.
We identify that the mass anisotropy (α), measuring the quadrupolar distortion in the Luttinger system, can be a useful non-thermal tuning parameter to further explore the territory of strongly interacting spin-3/2 fermions [28,43]. In particular, we find that strong quadrupolar distortions can be conducive for various magnetic orderings even at zero temperature. Specifically, when the electronic dispersion along the C 3v axes becomes almost flat (realized as α → π 2 ) the singlet A 2u magnetic order stabilizes at t = 0, see Figs. 6(c) and 16(c). On the other hand, when Luttinger fermions are almost non-dispersive along the C 4v axes (realized when α → 0), the system becomes susceptible toward the formation of a triplet T 1u magnetic order, see Figs. 7 and 16(d). For these two limiting scenarios the above two magnetic orders become (almost) mass [see Secs. V A 2 and V A 3], and their nucleation becomes energetically beneficiary even at zero temperature.
Irrespective of these details, we realize that all quantum phase transitions from the Luttinger semimetal to symmetry-breaking phases are continuous and controlled by various critical critical points, see Sec. VI A and Table II. To the leading order in the -expansion, the universality class of these transitions is characterized by the (a) correlation length exponent ν −1 = and (b) dynamic scaling exponent z = 2 [see Sec. VI B]. The presence of such quantum critical points manifests itself through the scaling of the transition temperature t c ∼ |δ| νz , yielding t c ∼ |δ| 2 for d = 3 or = 1, see Fig. 17, where δ is the reduced distance from the critical point.
Finally, we introduce (chemical) doping as another non-thermal tuning parameter to map out the global phase diagram of an interacting Luttinger metal, see Sec. VI D. Since any paired state maximally gaps the Fermi surface, its appearance at the lowest temperature is quite natural, at least when |µ| > 0. By contrast, excitonic phases (insulators or semimetals) become metallike (possessing a finite density of states) at finite chemical doping, according to the Luttinger theorem [78]. Therefore, particle-hole orders are accompanied by higher entropy due to the presence of a Fermi surface (with a constant density of states). To demonstrate the energyentropy competition in a metallic system, we choose the mass anisotropy parameter α such that at µ = 0 the system only supports excitonic orders. Upon raising (lowering) the chemical potential to the conduction (valence) band, we observe that a superconducting order develops at low temperature and the excitonic order gets pushed toward higher temperature and stronger interactions, see Figs. 5 and 6. Therefore, the overall structure of the global phase diagram is compatible with the energy-entropy competition, dictated by the emergent band topology of competing broken symmetry phases.
Furthermore, we also identify a definite "selection rule" among competing phases [see Sec. VI E]. From multiple cuts of the global phase diagram of a collection of strongly interacting spin-3/2 fermions, we find that two phases can reside in close vicinity of each other if the order-parameters describing two distinct phases can be combined to form a composite order-parameter. In other words, two ordered phases, respectively described by O(M 1 ) and O(M 2 ) symmetric order-parameters [see footnote 13], can reside next to each other only if one can construct an O(M ) symmetric composite-vector, where M 1 , M 2 ≤ M ≤ M 1 + M 2 , from the elements of two individual order parameters. This is so because when an interaction favors O(M 1 ) symmetry breaking order, it also enhances the scaling dimension of an O(M 2 ) symmetry breaking order and vice-versa, when the above selection rule is satisfied. Therefore, by tuning a suitable parameter (such as temperature, mass anisotropy parameter, chemical potential) one can induce a transition between two competing phases. As an immediate outcome of this selection rule we realize that while repulsive interactions A few specific cuts of the global phase diagram corroborate with ones extracted experimentally in Ln 2 Ir 2 O 7 [16] and half-Heusler compounds [36]. For exmaple, a finite temperature phase transition from A 2u or all-in allout magnetic order to a Luttinger semimetal has been observed in majority of 227 pyrochlore iridates (except for Ln=Pr) [16], and Fig. 3(c) qualitatively captures this phenomena (for strong enough g 3 ). By contrast, Pr 2 Ir 2 O 7 supports a large anomalous Hall effect below 1.5K [25][26][27]. Note that a triplet T 1u or 3-in 1-out magnetic order supports anomalous Hall effect due to the presence of Weyl nodes in the ordered state (see Sec. IV C and Ref. [28]) and the phase diagram from Fig. 3(d) shows the appearance of T 1u magnetic order at finite temperature for strong enough interaction (g 4 ). It is worth recalling that ARPES measurements strongly suggests that isotropic Luttinger semimetal describes the normal states of both Nd 2 Ir 2 O 7 and Pr 2 Ir 2 O 7 [10,11]. On the other hand, half-Heusler compounds LnPdBi display a confluence of magnetic order and superconductivity [36] and the phase diagrams shown in Figs. 5, 6(c) and 6(d) capture such competition (at least qualitative). This connection can be further substantiated from a recent penetration depth (∆λ) measurement in YPtBi [37], suggesting ∆λ ∼ T /T c (roughly) at low enough temperature, where T c = 0.78K in the superconducting transition temperature. Such a T-linear dependence of the penetration depth can result from gapless BdG fermions yielding (E) ∼ |E|. Any d-wave pairing, producing two nodal loops (see Table I), is therefore a natural candidate for the paired state in half-Heuslers (see also Refs. [49,51,53,54]). It is admitted that more microscopic analysis is needed to gain further insights into the global phase diagram of strongly interacting spin-3/2 fermions in various materials, which we leave for future investigation.
The energy-entropy competition and the proposed selection rule among competing orders provide valuable insights into the overall structure of the global phase diagram of strongly interacting spin-3/2 fermions. Various cuts of the phase diagram, which we exposed by pursuing an unbiased renormalization group analysis, corroborate (at least qualitatively) the former two ap-proaches. These approaches are not limited to interacting spin-3/2 fermions living in three dimensions. The methodology is applicable for a large set of strongly interacting multi-band systems, among which two dimensional Dirac (semi)metals, doped Bernal-stacked bilayer graphene (supporting bi-quadratic band touchings in two dimensions) [63,64,80], twisted bilayer graphene near so called the magic angle [69][70][71], three-dimensional Dirac or Weyl materials [68], doped topological (Kondo-)insulators (described by massive Dirac fermions) [65][66][67] and nodal loop metals [89,90] are the prominent and experimentally pertinent ones. In the future we will systematically study these systems, which should allow us to gain further insights into the global phase diagram of correlated materials, appreciate the role of emergent topology inside various broken symmetry phases, and search for possible routes to realize unconventional high temperature superconductors. From the above expression of the free-energy density we arrive the expression for specific heat [see Eq. (4)] and compressibility [see Eq. (5)]. From these two expressions we finally arrive at the following universal ratio reported in Eq. (6). This number is a characteristic of a z = 2 scale invariant fixed point in d = 3.

Appendix B: Dynamic conductivity
This appendix is devoted to disclose some key steps of the computation of the dynamic conductivity in Luttinger system. We focus on an isotropic system (for the sake of simplicity) and explicitly compute the zcomponent of the conductivity (σ zz ). Due to the cubic symmetry σ zz = σ xx = σ yy ≡ σ (say). To this end we use the Kubo formula and compute the polarization bubble as a function of external (Matsubara) frequency Π zz (iω n ) = − e 2 β m k Tr ĵ z G 0 (ip m , k)ĵ z where e is the electronic charge, β is the inverse temperature and the current operator along the z direction is given bŷ In the spectral representation the noninteracting Greens function reads as where Components ofd are shown in Eq. (C1). After performing the analytic continuation iω → ω + iη, we find the dynamic conductivity using the Kubo formula σ zz (ω) = Π zz (iω → ω + iη) ω . (B5) From the above expression we obtain both Drude and inter-band components of the dynamic conductivity, respectively shown in Eqs. (7) and (8). The universal function F Dr (x) appearing in the expression of the Drude component in Eq. (7) is given by The scaling of this function is shown in Fig. 21.
In this Appendix we present the Fierz reduction of the number of linearly independent quartic terms for the interacting Luttinger system. To perform this exercise for generic local density-density interactions we introduce a six component vector constituted by the quartic terms appearing in H int , see Eq. (35). The Fierz identity allows us to rewrite each quartic term appearing in X as a linear combination of the remaining ones according to where M and N are four-dimensional Hermitian matrices and for local interactions x = y. The set of sixteen matrices {Γ a , a = 1, · · · , 16} closes the U (4) Clifford algebra of four-dimensional matrices, and we choose Γ † a = Γ a = (Γ a ) −1 . The above Fierz constraint then takes a compact form F X = 0, where F is the Fierz matrix given by (D3) The rank of F is 3. Hence, the number of lienarly independent quartic terms is 3 = 6 (dimensionality of F ) -3 (rank of F ). We chose four-fermion interactions proportional to λ 0 , λ 1 and λ 3 as three linearly independent quartic terms. The remaining three quartic terms can then be expressed in terms of linear combinations of the above three according to Therefore, during the RG analysis whenever we generate any one of the above three quartic terms we can immediately express them in terms of the ones proportional