Spin-Driven Nematic Instability of the Multi-Orbital Hubbard Model: Application to Iron-Based Superconductors
M. H. Christensen, Jian Kang, B. M. Andersen, R. M. Fernandes
SSpin-Driven Nematic Instability of the Multi-Orbital Hubbard Model:Application to Iron-Based Superconductors
Morten H. Christensen,
1, 2
Jian Kang, Brian M. Andersen, and Rafael M. Fernandes School of Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA Niels Bohr Institute, University of Copenhagen, Juliane Maries Vej 30, DK-2100, Denmark
Nematic order resulting from the partial melting of density-waves has been proposed as themechanism to explain nematicity in iron-based superconductors. An outstanding question, however,is whether the microscopic electronic model for these systems – the multi-orbital Hubbard model –displays such an ordered state as its leading instability. In contrast to usual electronic instabilities,such as magnetic and charge order, this fluctuation-driven phenomenon cannot be captured by thestandard RPA method. Here, by including fluctuations beyond RPA in the multi-orbital Hubbardmodel, we derive its nematic susceptibility and contrast it with its ferro-orbital order susceptibility,showing that its leading instability is the spin-driven nematic phase. Our results also demonstratethe primary role played by the d xy orbital in driving the nematic transition, and reveal that high-energy magnetic fluctuations are essential to stabilize nematic order in the absence of magneticorder. The elucidation of electronic Ising-nematic order [1] –the state in which electronic degrees of freedom sponta-neously lower the point-group symmetry of the system –has become an important problem in unconventional su-perconductors [2, 3]. In both pnictides [4–9] and cuprates[10–12], the experimentally observed nematic order hasbeen proposed to arise from the partial melting of anunderlying spin density-wave (SDW) [13–16] or chargedensity-wave (CDW) [17–19] stripe-order. This mecha-nism is based on robust symmetry considerations. Con-sider for concreteness the stripe SDW case: the groundstate has an O (3) × Z degeneracy, with O (3) denot-ing the direction of the magnetic order parameter in spinspace, and Z denoting the selection of the SDW orderingvector Q X = ( π,
0) or Q Y = (0 , π ) (in the CDW case, thesystem has an O (2) × Z degeneracy). Fluctuations inlayered systems suppress the continuous ( O (3) or O (2))and the discrete ( Z ) symmetries differently, favoring anintermediate regime in which only the Z symmetry isbroken [13]. Because the Z symmetry distinguishes be-tween two ordering vectors related by a 90 ◦ rotation, itsbreaking implies a tetragonal-to-orthorhombic transition,and therefore nematic order.Although this mechanism for spin-driven (or charge-driven) nematic order has been established in simpli-fied low-energy models for pnictides [13–15, 20, 21] andcuprates [18, 19], it remains hotly debated whether morerealistic microscopic models display nematic order as theleading electronic instability. the cuprates, a sensiblemicroscopic model is the single-band Hubbard model,whose phase diagram has been reported to display ne-matic correlations in the strong-coupling regime [22, 23].For the pnictides, due to the 3 d configuration of Fe andto the small crystal field splittings, a five-orbital Hub-bard model, including Hund’s rule interactions, is a moreappropriate starting point [24, 25]. Furthermore, be-cause many pnictides display metallic behavior, a weak-coupling analysis of this intricate model can reveal impor-tant information about the underlying physics of thesematerials. Indeed, conventional RPA approaches have been employed to study the onset of SDW, CDW, andferromagnetism. However, in contrast to these usual elec-tronic instabilities, the standard RPA approach does notcapture the nematic instability even qualitatively, as weshow below, making it difficult to assess whether the re-alistic multi-orbital Hubbard model has a tendency to-wards nematic order.In this Letter, we extend the standard RPA ap-proach and derive the nematic susceptibility of an ar-bitrary multi-orbital Hubbard model. The fluctuationsincluded in this formalism arise solely from the non-interacting part of the Hamiltonian, such that interac-tions are treated at the same order as in the typical RPAmethod. We apply this formalism to the case of SDW-driven nematicity in iron pnictides, and establish that theleading instability of the five-orbital interacting model isa spin-driven nematic phase for a wide range of param-eters. In general, we find that nematic order exists ina narrow T range above the magnetic transition line, inagreement with experiments in the pnictides [26]. How-ever, magnetic fluctuations at higher energies can inducea sizable splitting between the two transitions, partic-ularly in the regime where the SDW transition is sup-pressed to zero. We propose that this effect may be rele-vant to understanding the unusual nematic phase of FeSe[27–31]. Previously, the investigation of the multi-orbitalHubbard model in Ref. [21] revealed the importance ofthe orbital content of the Fermi surface in the low-energyspin-nematic model of the pnictides. Here, we find fromthe orbitally-resolved nematic susceptibility that whereasthe d xz , d yz , and d xy orbitals contribute almost equally tothe SDW instability, the d xy orbital plays a stronger rolein driving the nematic instability. Finally, we comparethe nematic susceptibility with the RPA-derived ferro-orbital order susceptibility. work provides a promisingroute to search for nematicity in different compounds, asit is compatible with ab initio approaches and also withmethods that include the effects of moderate interactions,such as LDA+DMFT [32, 33].Our starting point is the multi-orbital Hubbard model a r X i v : . [ c ond - m a t . s up r- c on ] F e b with onsite interactions [25, 34]. The non-interactingpart is given by H = (cid:80) µ,ν ( (cid:15) µν ( k ) − ˜ (cid:15)δ µν ) c † k µσ c k νσ ,where c † k µσ creates an electron with momentum k andspin σ at orbital µ = 1 , ..., N orb and the hopping pa-rameters (cid:15) µν ( k ) are determined from tight-binding fitsto ab initio calculations (sums over spin and momen-tum indices are left implicit). The four onsite interac-tion terms correspond to the intra-orbital Hubbard term, H U = U (cid:80) µ n q µ ↑ n − q µ ↓ , the inter-orbital Hubbard term, H U (cid:48) = U (cid:48) (cid:80) µ<ν n q µσ n − q νσ (cid:48) , the Hund’s rule coupling, H J = J (cid:80) µ<ν c † k + q µσ c k νσ c † k (cid:48) − q νσ (cid:48) c k (cid:48) µσ (cid:48) , and the pair-hopping term H J (cid:48) = J (cid:48) (cid:80) µ<ν c † k + q µσ c † k (cid:48) − q µ ¯ σ c k (cid:48) ν ¯ σ c k νσ .These coefficients are related by U (cid:48) = U − J and J (cid:48) = J .Previous approaches considered the nematic susceptibil-ity of a spin-fermion model [35]; here, we will focus on theHubbard model within RPA. The mechanism in whichnematic order arises from the partial melting of an SDWor a CDW requires fluctuations at two momenta relatedby 90 ◦ , in general Q = (cid:0) πn , (cid:1) and Q = (cid:0) , πn (cid:1) , withinteger n . Although our formalism can be extended in astraightforward way to arbitrary n , hereafter we focus on n = 1. to make contact with the pnictides, we considerthe SDW channel. Performing a Hartree-Fock decouplingof H in both the q = 0 charge channel and the q = Q i SDW channel: H MF = (cid:88) k ( (cid:15) µν ( k ) − ˜ (cid:15) ν δ µν ) c † k µσ c k νσ − (cid:88) kq M i q µ · c † k − q + Q i µσ σ σσ (cid:48) c k µσ (cid:48) , (1)where ˜ (cid:15) ν incorporates the changes in the mean-field den-sities and M i q µ = (cid:80) k U ρµ (cid:104) c † k + q + Q i ρσ σ σσ (cid:48) c k ρσ (cid:48) (cid:105) are theSDW order parameters with i = X, Y . The interactionmatrix U ρµ is U aa = U and U ab (cid:54) = a = J . We consider onlyintra-orbital magnetism, since previous Hartree-Fock cal-culations revealed that in the ground state the intra-orbital SDW order parameters are the dominant ones[34]. In the standard RPA approach for the SDW in-stability, the electronic degrees of freedom are integratedout, yielding the quadratic magnetic free energy: F (2)mag [ M iµ ] = (cid:88) q,i = X,Y [ χ µνi ( q )] − M iq,µ · M i − q,ν , (2)with the magnetic propagator χ µνi ( q ) χ µνi ( q ) = (cid:2) ( U µν ) − + (cid:88) k G νµ ( k ) G µνi ( k + q ) (cid:3) − , (3)where G µνi ( k ) ≡ G µν ( k + Q i ) is the Green’s function inorbital basis, q = ( q , Ω n ), (cid:80) q = T /N q (cid:80) q (cid:80) Ω n , andΩ n = 2 nπT is the Matsubara frequency. The RPA mag-netic susceptibility (cid:10) M iq,µ · M i − q,ν (cid:11) is proportional to anddiverges at the same temperature as the magnetic prop-agator χ µνi ( q ). Note that the tetragonal symmetry of thesystem implies that a peak of χ µνi ( q ) at Q X = ( π,
0) will
FIG. 1. (Color online) Normal-state Fermi surface based onthe parameters of Ikeda et al . [44]. The colors indicate thedominant orbital contribution. be accompanied by an equal peak at Q Y = (0 , π ). There-fore, at this order in perturbation theory, the system doesnot distinguish the case in which either Q X or Q Y is se-lected (single- Q order) from the case in which both areselected (double- Q order), i.e. the standard RPA ap-proach is blind to nematicity. To remedy this problem,we go beyond the second-order expansion of the free en-ergy and calculate the quartic-order terms: F (4)mag [ M Xµ , M Yµ ] = 12 u ρνηµ (cid:0) M Xρ · M Xν + M Yρ · M Yν (cid:1) × (cid:0) M Xη · M Xµ + M Yη · M Yµ (cid:1) − g ρνηµ (cid:0) M Xρ · M Xν − M Yρ · M Yν (cid:1) × (cid:0) M Xη · M Xµ − M Yη · M Yµ (cid:1) + 2 w ρνηµ (cid:0) M Xρ · M Yν (cid:1) (cid:0) M Xη · M Yµ (cid:1) , (4)The quartic coefficients, whose expressions are shownexplicitly in Appendix A, depend only on the non-interacting Green’s functions. Although interactions canalso contribute to them, as shown in Refs. [36, 37],within the RPA approach these contributions are sub-leading and can be neglected. The most relevant coef-ficient for the nematic instability is g ρνηµ , whose termdistinguishes between single- Q and double- Q order.Specifically, a Hubbard-Stratonovich transformation ofthis term reveals the nematic order parameter (cid:104) φ µν (cid:105) ∝ (cid:10) M Xµ M Xν (cid:11) − (cid:10) M Yµ M Yν (cid:11) , a rank-2 tensor in orbital spacethat breaks the tetragonal symmetry of the system bymaking X (cid:54) = Y . The term with coefficient u ρνηµ is re-lated to Gaussian magnetic fluctuations in both SDWchannels, while the term with coefficient w ρνηµ mainlydistinguishes between the two types of double- Q order[36]. Eq. (4) is the multi-orbital generalization of themagnetic free energy previously obtained in effective low-energy models in the band basis, where the coefficient g becomes a scalar [13].It is now possible to compute the static nematic sus-ceptibility χ ρνηµ nem ∝ (cid:104) φ ρν φ ηµ (cid:105) in the paramagnetic phase FIG. 2. (Color online) The largest eigenvalues λ of (a) the bare nematic susceptibility χ ρνηµ nem , , the Q X/Y magneticpropagator χ mag , and (b) the full nematic susceptibility χ ρνηµ nem as a function of T for the case n = 6. The inset in (a) showsthe upturn of the magnetic susceptibility as it diverges. (see Appendix B for details of the derivation): χ ρνηµ nem = χ ηαµβ nem , (cid:16) δ ρβ δ να − g ρνγδ χ γαδβ nem , (cid:17) − , (5) χ ρνηµ nem , ≡ (cid:88) q,i = X,Y χ ρνi ( q ) χ ηµi ( − q ) . (6)The orbitally-resolved nematic susceptibility χ ρνηµ nem is arank-4 tensor that generalizes the scalar nematic suscep-tibility derived previously for effective low-energy models[38–40, 43]. The impact of the magnetic fluctuations en-coded in the coefficient g ρνγδ is clear: if this term wereabsent, then the (bare) nematic susceptibility would bemerely a higher-order convolution of the magnetic prop-agator, χ ρνηµ nem , , and therefore diverge at the same T asthe SDW susceptibility. To establish whether the ne-matic susceptibility diverges already in the paramagneticphase, one needs to compute its leading eigenvalue λ ( n ) from χ ρνηµ nem Φ ( n ) ρν = λ ( n ) Φ ( n ) ηµ , with n = 1 , ..., N . Thestructure of the corresponding eigen-matrix Φ ( n ) ηµ revealswhich orbitals promote the nematic instability, and whichorbitals favor a double- Q structure with no underlyingnematicity. We note that in principle the Gaussian fluc-tuations associated with the term with coefficient u ρνηµ can also renormalize the magnetic propagator χ ρνi . How-ever, because this effect merely renormalizes the SDWtransition temperature, we do not include it hereafter.Equation (5) is the RPA-generalized nematic suscep-tibility, which can be compared on equal-footing withother RPA instabilities of a weakly-interacting systemdescribed by a multi-orbital Hubbard model. We applythis formalism to a five-orbital model for the iron-basedsuperconductors and contrast the nematic susceptibilityto the ferro-orbital RPA susceptibility. The hopping pa-rameters are those from Ref. [44], whereas the interac-tions are set to U = 0 .
95 eV and J = U/ n = 6 ispresented in Fig. 1, consisting of three hole pockets atthe center and the corner of the Brillouin zone, and twoelectron pockets at the borders of the Brillouin zone. that FIG. 3. (Color online) Color plot of the normalized ele-ments of the eigen-matrix Φ ( n ) ηµ corresponding to the leadingeigenvalue of the bare (left) and of the full (right) nematic sus-ceptibilities. The dominant contributions arise from the d xz , d yz , and d xy orbital, with the d xy being the most importantfor nematicity. the d xy hole pocket at ( π, π ) is not present in all materi-als, as it depends on the Fe-As distance [41, 42].We evaluate Eqs. (3) and (5) numerically as functionsof T for various values of the occupation number n . Con-sider first n = 6: in Fig. 2(a), we plot the T dependenceof the largest eigenvalue of the static magnetic propaga-tor χ µνi (0) as well as the largest eigenvalue of the bare nematic susceptibility χ ρνηµ nem , . Despite having different T dependencies, both eigenvalues diverge at the same tem-perature T mag , confirming our assertion that the stan-dard RPA is blind to the nematic instability. In Fig. 2(b),we plot the largest eigenvalue of the full nematic suscep-tibility χ ρνηµ nem , as given by Eq. (5). Clearly, the eigen-value diverges at T > T mag : this is exactly the nematictransition temperature T nem .Interestingly, our results reveal a relatively small split-ting between T nem and T mag , with T nem ≈ . T mag ,which resembles the small T -range in which a nematic-paramagnetic phase is observed experimentally in theiron pnictides [26]. We caution, however, that this valueshould be understood as an upper boundary for the split-ting between the nematic and the actual magnetic tran-sition, since (cid:101) T mag calculated inside the nematic state isgenerally larger than T mag calculated in the tetragonalstate. Furthermore, the value for T mag obtained via RPAoverestimates the actual transition temperature due tothe absence of Gaussian fluctuations, as discussed above.While the largest eigenvalue λ ( n ) determines T nem , thestructure of the corresponding 5 × ( n ) ηµ re-veals the orbital-resolved nematic order parameter driv-ing the transition, since Φ ( n ) ηµ ∝ (cid:10) M Xη M Xµ (cid:11) − (cid:10) M Yη M Yµ (cid:11) .In Fig. 3 we plot the normalized elements of the leadingeigen-matrix Φ ( n ) ηµ for both the full and the bare nematicsusceptibility – which, as shown above, contains informa-tion only about the magnetic instability. In both cases,the dominant processes involve the d xz , d yz , and d xy or-bitals.There is however one important difference: the rela-tive weight of the d xy orbital is larger for χ ρνηµ nem than for FIG. 4. (Color online) (a)
Occupation number-temperature( n, T ) phase diagram for the bare magnetic and nematic phasetransitions, evidencing the narrow region displaying nematic-paramagnetic order. The solid T nem line takes into accountonly the contribution from low-energy (Ω n = 0) magnetic fluc-tuations, whereas the dashed line includes contributions fromhigher energies (Ω < Ω c = 1 eV). For n < .
75, an incom-mensurate magnetic order appears. (b)
Ferro-orbital ordersusceptibility χ oo a function of T for various values of the oc-cupation number n . In contrast to the nematic susceptibilityshown in Fig. 2, χ oo is nearly featureless and T -independentat low energies. χ ρνηµ nem , , i.e. while the three orbitals seem to contributeequally to drive the magnetic instability, the d xy orbitalplays a more important role in driving the nematic insta-bility. We interpret this in terms of the nesting propertiesof the orbital content of the Fermi surface in Fig. 1: whilethe d xy hole-pocket at ( π, π ) can form a single- Q SDWby combining with either the X or Y electron-pockets,since both have d xy spectral weight, the two d xz / d yz hole-pockets at (0 ,
0) can form a double- Q SDW by combiningwith both the X and Y pockets, since they have d yz and d xz spectral weight, respectively.Having analyzed the n = 6 case, we present in Fig. 4(a)the complete ( n, T ) phase diagram for the magnetic andnematic transitions. We restrict our analysis to n > . T mag is not peaked at n = 6. This is likely due to theabsence of disorder effects introduced by doping, whichare known to suppress T mag [45, 46]. Most importantly,across the entire phase diagram the nematic transitionline tracks closely the magnetic transition line, in agree-ment with the phase diagrams of the iron pnictides.An important issue in obtaining this phase diagram isthat, as shown in Eq. (6), the computation of the ne-matic susceptibility requires summing the magnetic fluc-tuations not only over the entire Brillouin zone, but alsoover energy (i.e. over Matsubara frequencies). Althoughthe propagator χ µνi ( q , Ω n ) is strongly peaked at Ω n = 0(see Appendix C), within RPA it saturates to a finitevalue for large energies [see Eq. (3)], requiring a fre-quency cutoff Ω c . Near a finite- T magnetic transition,due to the very sharp peak in χ µνi (cid:0) Q X/Y , Ω n (cid:1) , it is rea-sonable to take only the Ω n = 0 contribution – the low- energy magnetic fluctuations – resulting in the solid lineof Fig. 4. However, near the region where T mag → n (cid:54) = 0)is not justified. To address this problem, we introduce acutoff Ω c = 1 eV, at which the propagator reaches valuesclose to its saturation value, as shown in Appendix C.The corresponding nematic transition line is shown as adashed line in Fig. 4. Near the regime where the magnetictransition takes place at finite T , the only effect of thecutoff is to increase the nematic transition temperature,as expected. However, near the regime where T mag → T nem depends on the cutoff value, the main result isthat higher-energy magnetic fluctuations are essential topromote nematic order without magnetic order. In thisregard, it is interesting to note that, in FeSe, the onlyparent material in which nematic order is observed inthe absence of magnetic order, NMR measurements findno evidence for low-energy magnetic fluctuations [47, 48],whereas neutron scattering reports sizable fluctuations atmodest energy values [49, 50].A remaining question is whether or not the spin-drivennematic instability is the leading instability of the sys-tem. In particular, an ongoing debate [2, 37, 51–54] con-cerning iron-based materials is whether ferro-orbital or-der, signaled by an unequal occupation of the d xz and d yz orbitals, ∆ n ≡ n xz − n yz (cid:54) = 0, could drive the nematictransition, instead of the spin-driven mechanism exploredabove. To investigate this issue, we calculate the q = 0static component of the RPA orbital order susceptibility, χ oo ( q ) = (cid:104) ∆ n ( q )∆ n ( − q ) (cid:105) for the multi-orbital Hubbardmodel [31], of which a brief derivation is included in Ap-pendix D. As shown in Fig. 4(b), our results reveal anearly T -independent χ oo for the doping range and in-teractions investigated. This is not unexpected, since forreasonable values of U and J , there is no attraction inthe RPA charge channel. Therefore, within RPA, ferro-orbital order is unable to drive the nematic instability.Of course, once the coupling to magnetic fluctuations isincluded, which requires going beyond RPA, χ oo will di-verge at the same T as χ nem [21, 31, 55]. In this regard,by effectively decoupling these two channels, RPA pro-vides an interesting route to investigate which instabilityis the leading one – at least for weak or moderate inter-actions.In summary, we developed an appropriate extensionof the RPA approach to obtain the orbital-resolved spin-driven nematic susceptibility of an arbitrary multi-orbitalHubbard model. Application to the case of iron-basedsuperconductors reveals that the leading instability ofthe system is an interaction-driven nematic phase. The d xy orbital plays a leading role in promoting the ne-matic instability, and higher-energy magnetic fluctua-tions are essential to stabilize nematic order in the ab-sence of long-range magnetic order. Comparison withother RPA susceptibilities reveals that the nematic andmagnetic transitions follow each other closely, and thatthe ferro-orbital susceptibility does not diverge on itsown. More generally, our formalism can also be combinedwith first-principle approaches to search for other materi-als that may display electronic nematicity. Furthermore,because interactions appear only in the determination ofthe magnetic propagator, Eq. (3), this formalism canbe combined with other approaches that specifically in-clude moderate electronic interactions, such as DFT+Uor LDA+DMFT [32, 33]. ACKNOWLEDGMENTS
The authors are grateful to A. V. Chubukov, L. Fan-farillo, M. Sch¨utt, P. Orth, and B. Valenzuela for dis-cussions. MHC and BMA acknowledge financial supportfrom a Lunbeckfond fellowship (grant A9318). RMF andJK are supported by the U.S. Department of Energy,Office of Science, Basic Energy Sciences, under awardnumber de-sc0012336.
Appendix A: Fourth order coefficients
To derive the form of the free energy given in Eqs.(2) and (4) in the main text, we perform a Hubbard-Stratonovich (HS) decoupling thereby obtaining theelectron-mediated interactions between the magnetic or-der parameters. Formally the HS decoupling relies oninserting unity in the partition function, where unity, in the present case, is given by = (cid:90) D [ M Xµν , M Yµν ]exp (cid:20) − (cid:90) q (cid:0) M Xµν ( q ) (cid:0) U − (cid:1) µνρλ M Xρλ ( − q )+ M Yµν ( q ) (cid:0) U − (cid:1) µνρλ M Yρλ ( − q ) (cid:1)(cid:21) , (A1)and (cid:82) D [ M Xµν , M Yµν ] is chosen such that the path-integralevaluates to unity and q = ( q , Ω n ) (Ω n being a bosonicMatsubara frequency). The electrons are then integratedout resulting in an effective action for the magnetic order S eff [ M Xµν , M Yµν ] = (cid:88) i (cid:90) q M iµν ( q ) (cid:0) U − (cid:1) µνρλ M iρλ ( − q ) − Tr ln (cid:2) G µν ( k ) − − V µν ( q ) (cid:3) , (A2)where i = X, Y , µ and ν are orbital indices, k = ( k , ω n ), ω n = (2 n + 1) πT is the fermionic Matsubara frequency,and the trace is over all external indices (the spin indiceshave been suppressed, the Green’s function is diagonal inspin). G µν ( k ) is the matrix Green’s function, obtainedfrom the first term in Eq. (1) of the main text, and V originates from the coupling between the magnetic orderparameters and the electrons, the second term. In thebasis Ψ( k ) = ψ ( k ) ψ ( k + Q X ) ψ ( k + Q Y ) ψ ( k + Q X + Q Y ) (A3)these are given by the matrices G µν ( k ) = G µν ( k + q ) 0 0 00 G µν ( k + q + Q X ) 0 00 0 G µν ( k + q + Q Y ) 00 0 0 G µν ( k + q + Q X + Q Y ) (A4) V µν ( q ) = − M Xµν ( q ) · σ αβ − M Yµν ( q ) · σ αβ − M Xµν ( q ) · σ αβ − M Yµν ( q ) · σ αβ − M Yµν ( q ) · σ αβ − M Xµν ( q ) · σ αβ − M Yµν ( q ) · σ αβ − M Xµν ( q ) · σ αβ , (A5)where each element of the matrices should be understoodas an N orb × N orb matrix in orbital space, with the Greenfunction being G µν ( k ) = (cid:88) m (cid:104) µ | m (cid:105)(cid:104) m | ν (cid:105) iω n − ξ m ( k ) , (A6) where m refers to band basis and µ, ν refer to orbitalbasis. Expanding the trace-log to fourth order in themagnetic order parameters and applying the Pauli matrixidentity σ iαβ σ jβδ σ kδγ σ lγα = 2 (cid:0) δ ij δ kl − δ ik δ jl + δ il δ jk (cid:1) (A7)yields the magnetic free energy as written in Eqs. (2) and(4) of the main text, with the fourth order coefficients u ρνηµ = 116 (cid:88) k (cid:16) G µρ G ρνX G νη G ηµX − G µρ G ρηX G ην G νµX + G µρ G ρνX G νη G ηµY + G νρ G ρµX G µηX + Y G ηνX − G µρ G ρηX G ηνX + Y G νµY (cid:17) + ( X ↔ Y ) , (A8) g ρνηµ = − (cid:88) k (cid:16) G µρ G ρνX G νη G ηµX − G µρ G ρηX G ην G νµX − G µρ G ρνX G νη G ηµY − G νρ G ρµX G µηX + Y G ηνX + G µρ G ρηX G ηνX + Y G νµY (cid:17) + ( X ↔ Y ) , (A9) w ρνηµ = 116 (cid:88) k (cid:16) − G µρ G ρηX G ην G νµY + 2 G νρ G ρηX G ηµ G µνY − G ηρ G ρµX G µνX + Y G νηX + 2 G ηρ G ρνX G νµX + Y G µηX + G ρµ G µηY G ηνX + Y G νρX + G ρν G νηY G ηµX + Y G µρX + G µρ G ρνX G νηX + Y G ηµY + G νρ G ρµX G µηX + Y G ηνY (cid:17) , (A10)where repeated orbital indices are not summed. Hereall the Green functions are implicit functions of k and G µνj ( k ) = G µν ( k + Q j ) and (cid:80) k = T /N k (cid:80) k (cid:80) ω n . Appendix B: Nematic susceptibility
Preparing for an additional HS-decoupling we intro-duce two bosonic fields ψ ρν and φ ρν with the partitionfunction Z = (cid:90) D φ D ψ exp (cid:104)
12 ( u ρνηµ ) − ψ ρν ψ ηµ −
12 ( g ρνηµ ) − φ ρν φ ηµ (cid:105) , (B1)with integration measures chosen appropriately such that Z = 1. By performing the shifts ψ ρν → ψ ρν − u ρνηµ (cid:0) M Xη · M Xµ + M Yη · M Yµ (cid:1) , (B2) φ ρν → ψ ρν + g ρνηµ (cid:0) M Xη · M Xµ − M Yη · M Yµ (cid:1) , (B3)the terms quartic in M cancel accordingly. Following thestandard procedure we introduce a field ( h ρν ) conjugateto M Xρ · M Xν − M Yρ · M Yν and define (cid:101) φ ρν = φ ρν + h ρν .The resulting action is then S [ M iµ , ψ µν , φ µν ] = (cid:88) q,i = X,Y ( r µνi ( q ) + ψ µν ) M iµ · M iν −
12 ( u ρνηµ ) − ψ ρν ψ ηµ + 12 ( g ρνηµ ) − (cid:16) (cid:101) φ ρν − h ηµ (cid:17) (cid:16) (cid:101) φ ηµ − h ηµ (cid:17) − (cid:101) φ ρν (cid:0) M Xρ · M Xν − M Yρ · M Yν (cid:1) . (B4)Here r µνi ( q ) = ( U µν ) − + (cid:80) k G νµ ( k ) G µνi ( k + q ) and G µνi ( k ) ≡ G µν ( k + Q i ). It is now straightforward to com-pute the nematic susceptibility: χ ρνηµ nem = lim h → (cid:18) δ ln Z δh ρν δh ηµ (cid:19) = (cid:0) g ρνικ g ηµφλ (cid:1) − (cid:104) φ ικ φ φλ (cid:105) − ( g ρνηµ ) − , (B5) where we used the fact that (cid:104) φ ρν (cid:105) = 0 as we are abovethe nematic instability. To continue we note that δ Fδφ ρν δφ ηµ = (cid:104) φ ρν φ ηµ (cid:105) − , (B6)where the free energy is F = − T ln Z , (B7)obtained by integrating out the magnetic degrees of free-dom and taking the large N limit. We find the effectiveaction S eff [ ψ µν , φ µν ] = 12 ( g ρνηµ ) − φ ρν φ ηµ + 12 Tr ln (cid:2) χ − ικ,Y χ − κλ,X − φ ικ φ κλ + χ − ικ,Y φ κλ − φ ικ χ − κλ,X (cid:3) , (B8)where we have ignored the Gaussian fluctuations ψ ρν and( χ µνi ( q )) − = r µνi ( q ). Finally (cid:104) φ ρν φ ηµ (cid:105) − = ( g ρνηµ ) − − (cid:88) q,i = X,Y χ ρµ,i ( q ) χ νη,i ( − q ) (B9)and after some manipulations we arrive at the expressiongiven in the text for the nematic susceptibility. Appendix C: Frequency dependence of the magneticsusceptibility
In this section we illustrate the frequency depen-dence of the magnetic propagator at various tempera-tures for representative filling factors of the ( n, T ) phasediagram (Fig. 4(a) of the main text). Because themagnetic propagator peaks at ( π, / (0 , π ), we focus on Q X . For n = 5 .
90 as we approach the instability (at k B T = 45 meV), the frequency dependence of the prop-agator (cid:80) µν χ µνX ( Q X , Ω n ) has the form shown in Fig.5, where the bosonic Matsubara frequency is given byΩ n = 2 πnT . The gray area denotes the region included FIG. 5. Frequency dependence of the magnetic propaga-tor (cid:80) µν χ µνX ( Q X , Ω n ) for n = 5 .
90 at different temperatures.The parameters used are quoted in the main text. The mag-netic instability takes place at k B T = 45 meV. From (a) wesee that the contribution to the bare nematic susceptibilitycomes mostly from the zero frequency part of the magneticsusceptibility.FIG. 6. Frequency dependence of the magnetic propaga-tor (cid:80) µν χ µνX ( Q X , Ω n ) for n = 6 .
04 at different temperatures.The parameters used are quoted in the main text. As is ev-ident in (a), the peak broadens as zero temperature is ap-proached. However, even at higher temperatures, shown in(b) and (c), finite Matsubara frequencies provide considerablecontributions to the bare nematic susceptibility. in the cut-off Ω c = 1 eV, and the dotted line indicates (cid:80) µν χ µνX ( Q X , Ω n → ∞ ). The plots in Fig. 5 justifythe statement made in the main text that near a finite-temperature magnetic transition, one can safely neglectthe higher frequency contributions.To illustrate the importance of including high fre-quency contributions in the case where magnetic order isabsent, in Fig. 6 we also plot the frequency dependenceof the magnetic propagator for n = 6 .
04. It is clear thatthe peak is broadened, implying that it is no longer jus-tified to ignore the contributions originating from finite frequencies.
Appendix D: Derivation of the ferro-orbital ordersusceptibility
Ferro-orbital order is characterized by the breaking ofthe degeneracy between the d xz and d yz orbitals. In theitinerant framework this is seen by an inequivalent oc-cupation of the two orbitals, i.e. n xz (cid:54) = n yz . Defining∆ n ( q ) ≡ n xz ( q ) − n yz ( q ) as in the main text, the ferro-orbital susceptibility is given by (cid:104) ∆ n ( q )∆ n ( − q ) (cid:105) . Usingthe definition of ∆ n ( q ), we find that this is nothing but alinear combination of specific components of the chargesusceptibility, ( χ c ) µνρλ . In the standard RPA approach,the full expression is [31, 55] χ oo = ( χ cRPA ) xz,xzxz,xz + ( χ cRPA ) yz,yzyz,yz − ( χ cRPA ) xz,yzxz,yz − ( χ cRPA ) yz,xzyz,xz , (D1)where the RPA charge susceptibility is given by the usualexpression [25]( χ cRPA ) µνρλ = (cid:16) [1 + χ U c ] − (cid:17) µδργ ( χ ) δνγλ , (D2)where χ is the standard particle-hole bubble( χ ( q )) µνρλ = − (cid:88) k G µν ( k ) G ρλ ( k + q ) (D3)and U c is the interaction matrix in the charge channel.The latter differs from the interaction in the SDW chan-nel and is given by ( a (cid:54) = b )( U c ) aaaa = U , (D4)( U c ) aabb = 2 U (cid:48) − J = 2 U − J , (D5)( U c ) abab = 2 J − U (cid:48) = 4 J − U , (D6)( U c ) baab = J (cid:48) = J . (D7)We note that, due to the implicit summation over re-peated indices in Eq. (D2), all orbitals contribute to theRPA orbital order susceptibility. The static part of Eq.(D1) at q = 0 is the quantity plotted in Fig. 4(b) in themain text. [1] P. Chandra, P. Coleman, and A. I. Larkin, Phys. Rev.Lett. , 88 (1990).[2] R. M. Fernandes, A. V. Chubukov, and J. Schmalian,Nat. Phys. , 97 (2014).[3] S. A. Kivelson, I. P. Bindloss, E. Fradkin, V. Oganesyan,J. M. Tranquada, A. Kapitulnik, and C. Howald, Rev.Mod. Phys. , 1201 (2003).[4] J. Chu, J. G. Analytis, K. De Greve, P. L. McMahon, Z.Islam, Y. Yamamoto, and I. R. Fisher, Science , 824(2010).[5] T.-M. Chuang, M. P. Allan, J. Lee, Y. Xie, N. Ni, S. L.Bud’ko, G. S. Boebinger, P. C. Canfield, and J. C. Davis, Science , 181 (2010).[6] J. Chu, H. Kuo, J. G. Analytis, and I. R. Fisher, Science , 710 (2012).[7] C. Mirri, A. Dusza, S. Bastelberger, M. Chinotti, L. De-giorgi, J. H. Chu, H. H. Kuo, and I. R. Fisher, Phys. Rev.Lett. , 107001 (2015).[8] M. Yi, D. Lu, J.-H. Chu, J. G. Analytis, A. P. Sorini,A. F. Kemper, B. Moritz, S.-K. Mo, R. G. Moore, M.Hashimoto, W. S. Lee, Z. Hussain, T. P. Devereaux, I.R. Fisher, Z.-X. Shen, Proc. Nat. Acad. Sci. 2011 ,6878 (2011). [9] S. Kasahara, H. J. Shi, K. Hashimoto, S. Tonegawa, Y.Mizukami, T. Shibauchi, K. Sugimoto, T. Fukuda, T.Terashima, A. H. Nevidomskyy, and Y. Matsuda, Nature , 382 (2012).[10] R. Daou, J. Chang, D. LeBoeuf, O. Cyr-Choini`ere, F.Lalibert´e, N. Doiron-Leyraud, B. J. Ramshaw, R. Liang,D. A. Bonn, W. N. Hardy, and L. Taillefer, Nature, ,519 (2010).[11] V. Hinkov, D. Haug, B. Fauque, P. Bourges, Y. Sidis, A.Ivanov, C. Bernhard, C. T. Lin, and B. Keimer, Science , 597 (2008).[12] M. J. Lawler, K. Fujita, L. Jhinhwan, A. R. Schmidt, Y.Kohsaka, K. C. Koo, H. Eisaki, S. Uchida, J. C. Davis,J. P. Sethna, and K. Eun-Ah, Nature , 347 (2010).[13] R. M. Fernandes, A. V. Chubukov, J. Knolle, I. Eremin,and J. Schmalian, Phys. Rev. B , 024534 (2012).[14] C. Fang, H. Yao, W.-F. Tsai, J. P. Hu, and S. A. Kivelson,Phys. Rev. B , 224509 (2008).[15] C. Xu, M. M¨uller, and S. Sachdev, Phys. Rev. B ,020501(R) (2008).[16] E. Abrahams and Q. Si, J. Phys.: Condens. Matter ,223201 (2011).[17] S. A. Kivelson, E. Fradkin, and V. J. Emery, Nature ,550 (1998).[18] L. Nie, G. Tarjus, and S. A. Kivelson, PNAS , 7980(2014).[19] Y. Wang and A. Chubukov, Phys. Rev. B , 035149(2014).[20] I. Eremin and A. V. Chubukov, Phys. Rev. B , 024511(2011)[21] L. Fanfarillo, A. Cortijo, and B. Valenzuela, Phys. Rev.B , 214515 (2015).[22] S. Okamoto, D. S´en´echal, M. Chivelli, and A.-M. S.Tremblay, Phys. Rev. B , 180511(R) (2010).[23] S.-Q. Su, and T. A. Maier, Phys. Rev. B , 220506(R)(2011).[24] K. Kuroki, S. Onari, R. Arita, H. Usui, Y. Tanaka, H.Kontani, and H. Aoki, Phys. Rev. Lett. , 087004(2008).[25] A. F. Kemper, T. A. Maier, S. Graser, H-P. Cheng, P. J.Hirschfeld, and D. J. Scalapino, New J. Phys. , 062001 (2009); D. C. Johnston, Adv. Phys. , 803(2010); J. Paglione and R. L. Greene, Nature Phys. ,645 (2010); P. C. Canfield and S. L. Bud’ko, Annu. Rev.Cond. Mat. Phys. , 27 (2010); H. H. Wen and S. Li,Annu. Rev. Cond. Mat. Phys. , 121 (2011).[27] A. V. Chubukov, R. M. Fernandes, and J. Schmalian,Phys. Rev. B , 201105(R) (2015).[28] J. K. Glasbrenner, I. I. Mazin, H. O. Jeschke, P. J.Hirschfeld, R. M. Fernandes, and R. Valenti, NaturePhys. in press , doi:10.1038/nphys3434.[29] F. Wang, S. Kivelson, and D.-H. Lee, Nature Phys. inpress , doi:10.1038/nphys3456.[30] R. Yu and Q. Si, Phys. Rev. Lett. , 116401 (2015). [31] Y. Yamakawa, S. Onari, and H. Kontani,arXiv:1509.01161.[32] H. Park, K. Haule, and G. Kotliar, Phys. Rev. Lett. ,137007 (2011)[33] A. Georges, L. de’ Medici, and J. Mravlje, Annu. Rev.Cond. Mat. Phys. , 137 (2013).[34] M. N. Gastiasoro and B. M. Andersen, Phys. Rev. B ,150506(R) (2015).[35] S. Liang, A. Moreo, and E. Dagotto, Phys. Rev. Lett. , 047004 (2013).[36] X. Wang, J. Kang, and R. M. Fernandes, Phys. Rev. B , 024401 (2015).[37] S. Onari and H. Kontani, Phys. Rev. Lett. , 137001(2012).[38] R. M. Fernandes, L. H. VanBebber, S. Bhattacharya, P.Chandra, V. Keppens, D. Mandrus, M. A. McGuire, B.C. Sales, A. S. Sefat, and J. Schmalian, Phys. Rev. Lett. , 157003 (2010).[39] M. Khodas and A. Levchenko, Phys. Rev. B , 235119(2015).[40] Y. Gallais and I. Paul, arXiv:1508.01319.[41] V. Vildosola, L. Pourovskii, R. Arita, S. Biermann, andA. Georges, Phys. Rev. B , 064518 (2008).[42] M. J. Calder´on, B. Valenzuela, and E. Bascones, Phys.Rev. B , 094531 (2009).[43] H. Yamase and R. Zeyher, New J. Phys. , 073030(2015).[44] H. Ikeda, R. Arita, and J. Kuneˇs, Phys. Rev. B ,054502 (2010).[45] S. Liang, C. B. Bishop, A. Moreo, and E. Dagotto,arXiv:1506.00309.[46] R. M. Fernandes, M. G. Vavilov, and A. V. Chubukov,Phys. Rev. B , 140512(R) (2012).[47] A. E. B¨ohmer, T. Arai, F. Hardy, T. Hattori, T. Iye, T.Wolf, H. v. L¨ohneysen, K. Ishida, and C. Meingast Phys.Rev. Lett. , 027001 (2015).[48] S.-H. Baek, D. V. Efremov, J. M. Ok, J. S. Kim, J. vanden Brink, and B. B¨uchner, Nat. Mater. , 210 (2014).[49] M. C. Rahn, R. A. Ewings, S. J. Sedlmaier, S. J. Clarke,and A. T. Boothroyd, Phys. Rev. B , 180501(R)(2015).[50] Q. Wang, Y. Shen, B. Pan, Y. Hao, M. Ma, F. Zhou, P.Steffens, K. Schmalzl, T. R. Forrest, M. Abdel-Hafiez, D.A. Chareev, A. N. Vasiliev, P. Bourges, Y. Sidis, H. Cao,and J. Zhao, arXiv:1502.07544.[51] F. Kr¨uger, S. Kumar, J. Zaanan, and J. van den Brink,Phys. Rev. B , 054504 (2009).[52] C. C. Chen, B. Moritz, J. van den Brink, T. P. Devereaux,and R. R. P. Singh, Phys. Rev. B , 180418 (2009).[53] W. Lv, J. Wu, and P. Phillips, Phys. Rev. B , 224506(2009).[54] C. C. Lee, W. G. Yin, and W. Ku, Phys. Rev. Lett. ,267001 (2009).[55] Z. Wang and A. H. Nevidomskyy, J. Phys.: Condens.Matter27