Magnetism, superconductivity, and spontaneous orbital order in iron-based superconductors: who comes first and why?

Magnetism and nematic order are the two non-superconducting orders observed in iron-based superconductors. To elucidate the interplay between them and ultimately unveil the pairing mechanism, several models have been investigated. In models with quenched orbital degrees of freedom, magnetic fluctuations promote stripe magnetism which induces orbital order. In models with quenched spin degrees of freedom, charge fluctuations promote spontaneous orbital order which induces stripe magnetism. Here we develop an unbiased approach, in which we treat magnetic and orbital fluctuations on equal footing. Key to our approach is the inclusion of the orbital character of the low-energy electronic states into renormalization group analysis. Our results show that in systems with large Fermi energies, such as BaFe2As2, LaFeAsO, and NaFeAs, orbital order is induced by stripe magnetism. However, in systems with small Fermi energies, such as FeSe, the system develops a spontaneous orbital order, while magnetic order does not develop. Our results provide a unifying description of different iron-based materials.

The proponents of either orbital or magnetic fluctuations put forward models in which the unwanted degree of freedom is quenched. In a class of models where spin degrees of freedom are quenched [6,8,24], density fluctuations with opposite signs on the Fe d xz and d yz orbitals are enhanced as the temperature is lowered. Consequently, below a temperature T s the occupations of the d xz and d yz orbitals become unequal, breaking the tetragonal symmetry of the system and triggering a structural transition. This C 4breaking orbital order can either have zero wave vector, in which case it does not break translational symmetry, or have a finite modulation. In the band basis, orbital order with zero momentum is a Pomeranchuk-type order in the d-wave charge channel (d-POM) [25,26]. Such an order has been extensively studied in recent years in the context of quantum criticality [8,27]. Orbital order with a finite momentum connecting hole and electron pockets also attracted interest [1,6,8,9] as fluctuations of such order parameter mediate superconductivity (SC) and favor a signpreserving s þþ SC order.
In models where orbital degrees of freedom are quenched, orbital order is a spin-off of stripe spindensity-wave (SDW) magnetism. Stripe SDW order breaks the tetragonal symmetry between the x and y directions in addition to breaking the spin-rotational symmetry [28]. It has been shown [10,18] that the breaking of the discrete tetragonal symmetry occurs prior to the breaking of the continuous spin-rotational symmetry, via the development of a composite Ising-nematic order. By symmetry arguments, this order induces orbital order [29]. Magnetic fluctuations that drive Ising-nematic order also favor a sign-changing s þ− SC [2,3,5,10,30].
Each set of models uses approximations that have been strongly questioned. Orbital models assume attractive local (Hubbard) intraorbital interaction [6,8], in variance with first-principles calculations [31], while magnetic models assume a priori that magnetic fluctuations are the strongest [4], and often neglect the orbital content of low-energy excitations [32,33]. In reality, however, magnetic and orbital degrees of freedom are coupled and affect each other [34,35].
In this work, we treat magnetism, superconductivity, and orbital order on equal footing. We use the analytical renormalization group (RG) technique, which is the most unbiased way to analyze how different interaction channels affect each other and what is the leading (and the subleading) instability in the system [25,32,33,[36][37][38][39]. We list potential instabilities in Fig. 1 and show how each reconstructs the fermionic states. We consider a model with repulsive intrapocket and interpocket interactions, like in earlier studies of the interplay between magnetism and superconductivity. However, in distinction to earlier works [18,22,32,33], we explicitly include into consideration the orbital composition of the low-energy electronic states. This allows us to consider fluctuations in the orbital channel on equal footing with fluctuations in the magnetic and superconducting channels.
The RG analysis allows one to understand how interactions in different channels evolve as fermions with higher energies are gradually integrated out, and how the system behaves at progressively lower energies, which for practical purposes means lower temperatures. In RG language, this evolution is described by the flow of the couplings in various channels. All previous RG-related studies, either analytical (parquet RG, denoted pRG) or numerical (functional RG, denoted fRG), found that at the beginning of RG flow the interaction in the SDW channel is the strongest, i.e., spin fluctuations increase first. Increasing spin fluctuations affect the interaction in the s þ− SC channel by providing an attractive pairing component. As a result, at low energies the s þ− channel is necessarily attractive, even if the bare interaction in this channel was repulsive. To see this, one does not even need to do full-scale RG analysis as attraction in the s þ− channel due to spin fluctuations appears already in RPA-type approaches, which use a subset of processes included in RG.
In this paper, we report new results that go beyond earlier studies. First, we argue that the "push" from magnetic fluctuations not only makes the SC interaction attractive, but it also gives rise to attraction in the orbital (d-POM) channel. The latter effect cannot be obtained within RPA and requires the full RG analysis. Second, if we compare running interactions in various channels, as was done in earlier RG works on FeSCs, we find that at low energies the largest attractive interactions are in the SDW and SC channels, while the interaction in the d-POM channel is a distant third [25]. We argue, however, that the result changes in two aspects once we extend the RG analysis and compute susceptibilities in SDW, SC, and d-POM channels and compare the exponents. Then we find that (i) the SDW susceptibility does not actually diverge down to the lowest energies as the growth of SDW fluctuations is halted by feedback effects from SC and d-POM fluctuations and (ii) the susceptibility in d-POM channel has a larger exponent than the one in s þ− SC channel. The outcome of this analysis is that, upon lowering the FIG. 1. Low-energy states and potential instabilities. The orbital content of the 2D Fermi surface of the Fe-based superconductors is plotted together with the changes in the fermionic excitations promoted by one of three electronic instabilities-s þ− superconductivity, stripe SDW magnetism, and nematicity (breaking of C 4 lattice rotational symmetry), which necessarily gives rise to orbital order. The low-energy excitations live near hole pockets centered at the Γ point (k x ¼ k y ¼ 0), and near electron pockets centered at ð0; πÞ ðπ; 0Þ in the 1-Fe Brillouin zone. Excitations near the hole pockets are made out of d xz and d yz orbitals, while the ones near the ð0; πÞ [ðπ; 0Þ] electron pockets are made out of d xz and d xy (d yz and d xy ) orbitals [3,4]. In some systems, there exists a third hole pocket (not shown) centered at ðπ; πÞ and made out of the d xy orbital. s þ− superconductivity gaps out low-energy excitations, and the superconducting order parameter changes sign between hole and electron pockets. Stripe SDW magnetism with momentum ð0; πÞ or ðπ; 0Þ (shown) mixes hole and electron states by band folding and splits hole and electron pockets into even smaller subpockets. Orbital order elongates the two hole pockets in opposite directions and makes one electron pocket larger and the other one smaller. temperature, the system first develops d-POM order (i.e., a spontaneous orbital order), then, at a lower T, it develops s þ− SC, while SDW order does not develop down to T ¼ 0. This sequence of transitions is consistent with what is observed in FeSe [24], and we argue that nematicity in FeSe is the consequence of a spontaneous orbital order. We apply our results to FeSe in Sec. V. We, however, do not argue that this sequence of transitions should hold in all FeSCs because the RG flow extends only down to energies comparable to the (largest of) Fermi energies E F . If E F is small enough, as was argued in Ref. [40] based on several measurements, the leading instability (in our case, d-POM) occurs at T ins ∼ E F and is captured within RG. If, however, E F is larger, the RG flow runs only down to OðE F Þ, as at smaller energies different channels no longer "talk" to each other. In this situation, the d-POM channel remains weak, and the system develops either SDW or s þ− SC order. This is the case for most ironbased systems (see, e.g., Ref. [41], which detected E F ∼ 200 meV on electron pockets in weakly doped BaFe 2 As 2 ). In this situation, the nematic order is most likely promoted by composite stripe SDW fluctuations, which force the C 4 symmetry to break at a higher temperature than O(3) spinrotation symmetry (an Ising-nematic order).
The structure of the paper is the following. In Sec. II, we introduce the model and discuss the approximations. In Sec. III, we consider the RG flow for the couplings, vertices, and susceptibilities for the model in which we approximate the orbital content of the electron pockets as either pure d xz or pure d yz . In Sec. IV, we do the same analysis for the model in which we approximate the orbital composition of the two electron pockets as pure d xy . We show that the outcome of RG flow is the same in both cases. In Sec. V, we compare our results with the experiments, and in Sec. VI, we compare our approach with other theoretical scenarios for FeSe. We present our conclusions in Sec. VII. To keep the presentation focused, some aspects of our analysis are moved into Appendixes. Some technical details of RG analysis are presented in the Supplemental Material (SM) [42].

A. Free-fermion Hamiltonian
We depart from the microscopic model containing several neighbors hopping among all five 3d Fe orbitals and convert it into the band basis in order to distinguish high-energy and low-energy states. The orbital composition of the excitations does not show up in the kinetic part of the Hamiltonian in the band basis, but it imposes angular dependencies on the four-fermion interaction terms. As we show, the terms with different angular dependencies flow differently under RG.
The fermionic structure in the band basis contains hole and electron pockets. Two hole pockets are centered at the Γ point (k x ¼ k y ¼ 0) and are constructed out of d xz and d yz orbitals [43] (Fig. 1). To construct the low-energy Hamiltonian for states near hole pockets, we follow Ref. [44] and introduce the low-energy spinor wave function ψ † σ ðqÞ ¼ ½f † 1;σ ðqÞ; f † 2;σ ðqÞ; d † 1;σ ðqÞ; d † 2;σ ðqÞ, where the subscripts μ, ν ¼ 1, 2 refer to the xz and yz orbital content, respectively. We use 1-Fe Brillouin zone notations and neglect processes with momentum transfer ðπ; πÞ, which may be present due to the difference between the hopping via pnictogen or chalcogen atoms above and below the iron layer [45].
The quadratic part of the Hamiltonian is expressed in terms of the components of the spinor ψ † σ ðqÞ as The effective Hamiltonian for hole pockets is specified by where θ k ¼ arctanðk y =k x Þ, and ϵ Γ , 1=m Γ , a, and c are the parameters of the model, which are determined by comparison with the angle-resolved photoemission (ARPES) experiments. We emphasize that these parameters already absorb regular renormalizations of the dispersions due to interactions, including the renormalizations of the bandwidths [46] and terms that shift separately chemical potentials of hole and electron bands [47]. Such terms are not double counted in our approach because they predominantly come from fermions with energies comparable to the bandwidth. The RG analysis, on the other hand, accounts for singular logarithmic renormalizations coming from fermions with progressively smaller energies, but neglects self-energy contributions from these fermions. The transformation to the band basis is a rotation: For a ¼ c, the rotation angle θ kσ is identical to the angle θ formed by the vector k with a given axis, and the two hole Fermi surfaces (FSs) are circular. For simplicity, below we focus on this case. The extension to a ≠ c is straightforward and does not change the results and the conclusions. The kinetic energy of excitations around hole pockets in the band basis is The two dispersions are not identical when m c ≠ m h , but are degenerate by symmetry at k ¼ 0 in the absence of spin-orbit coupling [44]. The degeneracy implies that both Γ-centered hole pockets must be present simultaneously already in the minimal model.
The two electron pockets are centered at Q 1 ¼ ð0; πÞ and Q 2 ¼ ðπ; 0Þ in the 1-Fe Brillouin zone (Fig. 1). The kinetic energy of the fermions near the electron pockets is x =ð2m y Þ þ k 2 y =ð2m x Þ − μ, with k measured with respect to Q i for f i;k ≡ f i;kþQ i . The two electron pockets are related by C 4 symmetry and transform into each other under a π=2 rotation. The band fermions f 1;kþQ 1 and f 2;kþQ 2 are linear combinations of d xz =d xy and d yz =d xy orbitals, respectively [3,4]. The relative amplitude of the spectral weights depends on system parameters.
We assume that there is a substantial energy range of metallic behavior and do not discuss orbitally selective tendency towards electron localization [46].

B. Interacting Hamiltonian
We depart from the local Hubbard-Hund interaction, in the notations of Ref. [4]: where the index i enumerates the iron sites located at R i and is the annihilation operator of an electron at the iron site located at R j with spin σ in the orbital state labeled by μ (μ ¼ 1 and μ ¼ 2 refer to xz and yz orbitals, respectively). Further, n iμσ ¼ ψ † iμσ ψ iμσ is the density operator, n iμ ¼ n iμ↑ þ n iμ↓ , and N is the number of iron atoms. Equation (4) can be rewritten in an SU(2) invariant form as Substituting Eq. (5) into Eq. (6), we obtain Here, the momenta of the fermion operators in each term, k 1 , k 2 , k 3 , k 4 , are omitted for clarity. P 0 stands for the summation over the spin indices (σ, σ 0 ) and over fermion momenta subject to momentum conservation k 1 − k 2 þ k 3 − k 4 ¼ 0, and also includes the normalization factor 1=N. The initial observation, which sets the stage for the RG analysis, is that Eq. (7) is not the most general one consistent with the tetragonal symmetry. The actual number of topologically distinct invariant combinations of 4fermion terms involving states near the Fermi pockets is much higher and equals 30 for a generic 4-band model in the absence of spin-orbit coupling [44]. The bare values of all 30 couplings are expressed in terms of U, U 0 , J, J 0 , but under RG the couplings flow to different values. To make the problem analytically tractable, in Sec. III we neglect interactions involving fermions from d xy orbitals on electron pockets. This effectively implies that we approximate ð0; πÞ [ðπ; 0Þ] pockets as purely d xz [d yz ], i.e., set f 1;kþQ 1 ¼ d xz;kþQ 1 and f 2;kþQ 2 ¼ d yz;kþQ 2 . This approximation reduces the number of couplings to a manageable 14. We emphasize that this does not imply that we change the tight-binding electronic structure, as the approximation involves only low-energy fermions.
In Sec. IV, we consider the opposite case, keeping only the d xy spectral weight on the two electron pockets. This last approximation also reduces the number of couplings to 14. We show that the outcome of the RG flow for the case of purely d xy electron pockets is the same as with pure d xz (d yz ) pockets. This gives us confidence that we capture the right physics even with the reduced number of couplings.

III. MODEL WITH d xz =d yz ELECTRON BANDS
As we explain above, there are 14 topologically different interaction terms. They are One can verify that each term in Eq. (8) obeys the tetragonal symmetry separately.
Out of the 14 interactions, four involve fermions near the two hole pockets (U 4 ,Ū 4 ,Ũ 4 ,Ũ 4 ), another four involve fermions near the two electron pockets (U 5 ,Ū 5 ,Ũ 5 ,Ũ 5 ), and six involve fermions near both hole and electron pockets ( . Comparing Eqs. (7) and (8), we find the bare values of the couplings These relations, however, hold only for bare couplings. As we see, they are not preserved under the RG flow. On the other hand, the RG flow does not generate new couplings in addition to the 14 in Eq. (8); i.e., the model with 14 coupling is fully renormalizable.
In the band basis, these 14 different interaction parameters are the prefactors for 14 combinations of the original 152 interaction terms in the band basis (96 involving c and h fermions, eight involving f 1 and f 2 fermions, and 48 cross terms), combined using the symmetry condition that under rotation by π=2, c k → −h k , h k → c k , and f 1 → f 2 . We do not present the full expression (it is too lengthy), but show a representative set from each combination: where Á Á Á stand for other terms in each of the 14 combinations in Eq. (9).
MAGNETISM, SUPERCONDUCTIVITY, AND SPONTANEOUS … PHYS. REV. X 6, 041045 (2016) In the mean-field approximation, the bare values of these 14 couplings are used to compute susceptibilities in SDW, SC, d-POM, and other channels. A simple analysis shows that in this case SDW wins over SC and orbital order. However, the mean-field approximation is strongly questionable because it effectively isolates each electronic channel, neglecting their interplay and mutual feedback. To overcome this limitation, below we implement the parquet RG approach and calculate how the couplings and the susceptibilities in different channels evolve as high-energy degrees of freedom are progressively integrated out. In this approach, each dimensionless coupling u i ¼ ðA i =4πÞU i , where A i are combinations of effective masses, acquires a dependence on the running energy or temperature scale E via L ¼ log W=E, where W is the upper energy cutoff for the low-energy theory, which in a generic case is of the order of the bandwidth.

A. RG equations
The derivation of the one-loop RG equations is tedious but straightforward. We present the details of the derivations and the full equations in the SM [42] and here list the 14 pRG equations in the approximation m c One can immediately verify that the running couplings flow to different values under the RG, and these new values cannot be reexpressed just in terms of running U, U 0 , J, J 0 . This implies that the model with only onsite interactions does not survive under renormalization and longer-range interactions emerge in the process of the RG flow. We discuss this in more detail later in this section.

B. RG flow
The analysis of Eq. (10) readily reveals that the last four RG equations decouple from the other ten, and thatũ 4;5 andũ 4;5 flow to zero under RG. The remaining ten RG equations are all coupled and have to be solved selfconsistently. We analyze the RG equations for nonzero bare couplings (i.e., nonzero U, U 0 , J, J 0 ) and find that the system flows towards a single stable fixed trajectory, along which u i andū i become equivalent. We show the RG flow in Fig. 2. The equivalence of u i andū i along the stable fixed trajectory implies that the termsū i , which were originally of order U 0 or even J, grow under RG and eventually become comparable to u i , which were originally of order U. In other words, the initial hierarchy of interactions disappears under the RG flow. When electron pockets are approximated as pure d xy , the behavior along the stable fixed trajectory is the same, but some of u i andū i turn out to be equal already at the bare level; see Sec. IV.
W is the bandwidth, and E is the running energy (temperature) at which one probes the system. The flow of other couplings is similar. The couplings u 1 − u 5 and u 1 −ū 5 all diverge as 1=ðL 0 − LÞ when L approaches L 0 , whose value depends on the initial conditions. The couplingsũ 4 ,ũ 4 ,ũ 5 , u 5 tend to zero at L → L 0 . (b) Flow of the ratios of the couplings. All ratios tend to fixed finite values as L approaches L 0 : The ratios u 2 =u 1 andū 2 =u 1 tend to zero as L approaches L 0 . The initial values used areū 1 ffiffiffiffiffiffiffiffiffiffiffi ffi m e m h p Þ ¼ 1.1 for definiteness. For the model with d xy electron pockets, the fixed trajectory is the same, but the system approaches it much faster (see Sec. IV).
Along the stable fixed trajectory, the ratios between various couplings become pure numbers: Solving Eq. (10) for u 1 and γ i , we obtain We emphasize that this behavior is completely universal and does not depend on initial parameters U, U 0 , J, J 0 , except for L 0 , which is the energy or temperature scale at which couplings diverge and the system becomes unstable towards either SC or SDW or d-POM order. We also emphasize that the RG equations are valid up to L F ¼ log W=E F , where, we remind the reader, E F is the largest of the Fermi energies associated with the different Fermi pockets. For E < E F , particle-particle and particlehole channels no longer talk to each other and the flow equation is different (see below). Note in passing that the different flow of the 14 couplings in Eq. (9) under RG implies that the system self-generates nonlocal interactions. The information extracted from the low-energy sector only is not sufficient to fully specify which nonlocal interactions are generated, but the model with 14 couplings can be constructed if one also adds to the local interactions U, U 0 , J, J 0 the interactions of the same Hubbard and Hund type, but involving fermions from different sites of each plaquette on a square lattice. Thus, five terms involving fermions from the same orbital d xz or d yz (U 1 , U 2 , U 3 , U 4 , and U 5 terms) appear with different couplings if we also introduce, in addition to the on-site U, the terms and analogous (symmetry-related terms) for the d yz orbital. In Eq. (12), a x and a y are the components of the lattice spacing a. The couplings U i (i ¼ 1-5) are now given by One can easily verify that the interactions within a given plaquette involving fermions from different orbitals split the U 0 , J, and J 0 terms into subsets, each consisting of three different interactions [there are five terms in each subset, like in Eq. (12), but there are only three nonequivalent combinations of different U 0 i , J i , and J 0 i ].

C. Competition between channels
We now use the results for the RG flow to find in which of the many electronic channels the system actually develops an instability upon lowering the temperature T. For this, we introduce order parameters in different channels Δ i [i ¼ SDW, SC, POM, CDW (charge-density wave)] and also introduce infinitesimally small vertices Γ 0;i for the coupling between fermions and these order parameters. Next, we identify the combinations of the couplings We list all potential order parameters, bilinear in fermions, and all U i in Appendix B. These bilinear combinations form irreducible representations of the D 4h group separately for electron and hole pockets. In particular, Pomeranchuk order parameters are even under inversion and belong to one of four one-dimensional irreducible representations A 1g , A 2g , B 1g , B 2g . They include intraorbital and interorbital s-wave (A 1g and A 2g ) and d-wave order (B 1g and B 2g ). An intraorbital A 1g s-wave Pomeranchuk order parameter is characterized by separate order parameters on the hole pockets Δ h POM;s ¼  We list the relevant combinations of U i for SC, SDW, and POM in Table I. We see that along the fixed trajectory the interaction in the d-POM channel is where, we remind the reader, U 4 and U 5 are densitydensity intrapocket interactions on hole and on electron pockets (sign convention is such that attraction corresponds to positive sign of U SDW ; U SC , or U POM ). The bare values U 4;0 ¼ U 5;0 ¼ U are positive; hence, if we use the bare values, we would find that the interaction in d-POM channel is repulsive. However, in the process of RG flow, the couplings U 4 ¼ U 5 change sign and become negative due to the "push" from the SDW channel onto the orbital channel [see Fig. 2(b)]. As a result, the interaction in the d-POM channel becomes attractive below a certain energy. We emphasize that this result could not be obtained in the RPA approximation. We also emphasize that the interaction in the d-POM channel changes sign when the interaction in the s þ− SC channel is already attractive. The interaction in the s þ− SC channel is U SC ¼ −U 4 þ U 3 , and the attraction develops when the running U 4 gets smaller than U 3 . The coupling U 3 contributes to both SC and SDW channels, and it increases under RG [see Fig. 2(b)]. As a consequence, the pairing interaction U SC changes sign when U 4 is still positive.
Earlier RG studies assumed that the channel with the largest U i along the fixed trajectory wins, setting the hierarchy of instabilities based on the values of U i . We argue that this is not always true, and to compare different channels one actually needs to obtain and solve another set of RG equations for Γ i , then compute the corresponding susceptibilities, find which ones diverge, and compare the exponents. The leading instability will be in the channel in which the exponent is the largest. This procedure has been applied to other problems [37,48,49], but it has not been applied yet to analyze the interplay between orbital order and other orders in FeSCs. The advantage of using analytical RG for this procedure is that (i) the number of interaction is set by symmetry rather than the number of patches on the Fermi surface and (ii) the RG equations for Γ i and for the susceptibilities can be obtained in a straightforward way.
The analysis of the susceptibilities is different for the SC and SDW channels and the POM channel. For the SC and SDW channels, Π i is logarithmic, and, to logarithmic accuracy, where Γ SDW ðL 0 Þ and Γ sþ− SC ðL 0 Þ are the fully renormalized SDW and SC vertices obtained from the solutions of the RG equations _ Γ i ∝ Γ i U i . Along the fixed trajectory, where γ 3;4 are given by Eq. (11) (see SM [42] for details). Solving these two equations and substituting the results into Eq. (14), we obtain with the exponents where C ¼ ðm e þ m h Þ=ð2 ffiffiffiffiffiffiffiffiffiffiffi ffi m e m h p Þ ≥ 1. We plot α i as a function of C in Fig. 3(a), and show the schematic behavior of the susceptibilities in the SDW and SC (s þ− ) channels along the fixed trajectory in Fig. 3(b). We see that for all values of C, 0 < α SC < 1, while α SDW < 0. This implies that only SC order develops. SDW order does not develop, despite that at the bare level the SDW channel is the only attractive channel.
The phenomenon in which the SDW interaction creates an attraction in the SC channel and eventually SC wins had already been found in earlier numerical and analytical RG studies in multiband FeSCs [33,38] and in doped graphene TABLE I. The interactions in different channels along the stable fixed trajectory. Positive sign of an interaction implies an attraction. All interactions scale as 1=ðL 0 − LÞ and diverge at RG scale L 0 . We use these interactions to compute vertices and susceptibilities in SDW channel, s þ− superconducting channel, and s-wave and d-wave Pomeranchuk channels (POM s and POM d). [50,51]. We find in the analysis of the susceptibilities that not only SC wins over SDW, but the feedback effect from the rising superconducting fluctuations halts the growth of the SDW susceptibility. We now turn to the POM channel. Here, the situation is different because the particle-hole polarization bubble at energy E is determined by fermions with energies of order E. As a result, the Pomeranchuk susceptibilities obey algebraic rather than differential equations (see SM [42] for details). We find that the susceptibilities in the s þ− and d þ− subchannels increase under RG and behave as For both susceptibilities, the exponent α POM ¼ 1 is larger than α SC < 1. Furthermore, for all values of C, L P s and L P d are smaller than L 0 . As a result, χ s;d POM diverge with a larger exponent and at a smaller L than the SC susceptibility [see Fig. 3(b)]; i.e., within one-loop RG, the first instability upon lowering the temperature actually occurs in the Pomeranchuk channel. Note that at L ¼ L P s;d , u 1 ∼ 1, and the corrections to one-loop RG may become relevant. Still, the comparison of the susceptibilities clearly favors the POM channel over the SC and SDW channels. Also, number wise, for L ¼ L P d and C ¼ 1, u 1 ¼ ð1=16Þ=ðL 0 − LÞ ¼ 1=6, which is still a small number.
Of the two Pomeranchuk susceptibilities, the larger one is in the s þ− channel. This order has been analyzed in Ref. [25], where it was noticed that the corresponding coupling U s−POM ¼ −U 4 þ 4U 1 (see Table I) contains U 1 , which also contributes to SDW channel and grows under RG; i.e., U s−POM is positive (attractive) at all energies. The s-POM order with Δ h POM;s ¼ −Δ e POM;s splits the chemical potentials on the hole and electron pockets, but conserves the total number of carriers. Because this does not correspond to a true symmetry breaking, the divergence of χ s POM must be softened by terms beyond RG, such as regular terms in the fermionic self-energy [47]. Yet, the relative chemical potential shift μ h − μ e must be enhanced near the temperature at which χ s POM diverges within the RG. Interestingly, the analysis of ARPES data for several FeSCs did find [52] evidence for temperature-dependent μ h − μ e , which behaves much like order parameter below 300 K.
The true Pomeranchuk instability is in the d þ− channel, signaled by the divergence of χ d POM . This order, with sgnΔ h POM;d ¼ −sgnΔ e POM;d implies that the symmetry between the on-site energies of d xz and d yz orbitals is spontaneously broken, hd † xz d xz i ≠ hd † yz d yz i, and hd † xz d xz − d † yz d yz i changes sign between hole and electron pockets. Note that in our solution, the functiongðkÞ is not exactly −1; hence, the order parameter has both d þþ and d þ− components.
We emphasize that the origin of the d þ− orbital order is not just an attraction in the POM channel, as proposed by other works. In our case, the bare interaction may well be repulsive (when U þ J > 2U 0 , see above), yet the POM ffiffiffiffiffiffiffiffiffiffiffi ffi m e m h p Þ. The largest exponent α POM ¼ 1 is in the Pomeranchuk channel. The exponent α SDW < 0, which implies that within RG χ SDW does not diverge. (b) The schematic behavior of susceptibilities in SDW, s þ− SC, and POM channels. We use the couplings along the fixed trajectories as inputs for the RG equations for the susceptibilities. The Pomeranchuk susceptibility actually diverges at L ¼ L P < L 0 . As a result, the leading instability upon lowering the temperature is towards d-wave orbital ordering. s þ− superconductivity develops at a smaller T, and SDW instability does not develop. This holds when L 0 is smaller than L F ¼ log W=E F , i.e., if the instability develops at an energy or temperature larger than E F . If L F < L 0 , the RG flow runs up to L ¼ L F , and at larger L, SDW and SC channels decouple and develop independent of each other, while the Pomeranchuk channel gets frozen. In this situation, the system first develops either SDW or SC order, depending on the interplay between L F and L 0 and the degree of nesting. channel becomes attractive in the process of the RG flow and eventually wins over SC and SDW. The attraction in the POM channel is driven by the coupling to magnetic fluctuations, and in this respect the RG scenario for orbital ordering falls into the orbit of "magnetic scenarios" for nematicity.
Therefore, the full one-loop RG analysis shows that the system first develops d þ− orbital order at T s and then becomes a superconductor at a lower T c . SDW order does not develop. This sequence of transitions is fully consistent with that in FeSe. In other FeSCs, however, the system does develop SDW order at T N at small dopings, and the nematic transition line follows T N very closely, suggesting that nematic order is a vestige of the SDW order.
To understand the difference between FeSe and other FeSCs, we note that in our analysis we assume that the RG flow reaches the fixed trajectory at L ¼ L P d ≈ L 0 , before the RG analysis breaks down at an energy comparable to the largest E F in the system, i.e., at L ¼ L F . This holds when L 0 < L F , i.e., when all Fermi energies are very small compared to the bandwidth. This is the case of FeSe [40,53]. If L F < L 0 , the RG flow runs up only to L ¼ L F , and at larger L, different instability channels decouple from each other. As a result, the divergence of the Pomeranchuk susceptibility is cut and this channel no longer competes with SC or SDW. Also, because the SC and the SDW channels do not mix below E F , each develops independently in a mean-field fashion with the couplings taken at L ¼ L F (Refs. [5,33]). If L F is small enough, these values are close to the bare ones and the system develops SDW order (and Ising-nematic order above it, if SDW order is a stripe). When doping gets larger (and nesting gets weaker), the SDW channel becomes less singular and SC order develops first. This last behavior is consistent with the one observed in, e.g., 122 systems, for which the largest E F ≤ 200 meV (Ref. [41]) well exceeds T N , T c ∼ 10 meV.
We also emphasize that the behavior that we find along the fixed trajectory of the RG flow is universal in the sense that it holds for arbitrary ratios of U, U 0 , J, and J 0 , as long as they are all positive. The values of Hubbard and Hund interactions determine only how quickly the system reaches the fixed trajectory and how quickly the susceptibility in the orbital channel gets larger than the ones in SDW and SC channels. To make this point explicit, in Appendix C, we show the behavior of the SDW, SC, and Pomeranchuk susceptibilities along the full RG flow (i.e., not only along a fixed trajectory) for two different values of the couplings: the one for which the bare interaction in the orbital channel is attractive, and the one for which it is repulsive. As expected, the interplay between different channels is the same in both cases, but the susceptibility in the Pomeranchuk channel diverges faster when the bare interaction in this channel is attractive.

IV. MODEL WITH d xy ELECTRON POCKETS
We now consider another approximation for the full 4band model, in which we approximate the two electron pockets as purely d xy . The hole pockets are made out of d xz and d yz orbitals, like before.
One can use the same reasoning as in the previous section and construct the most generic model describing the interaction between low-energy fermions. We argue that this model again contains 14 topologically different interaction terms and is described by the same Eqs. (8) and (7), only now f fermions correspond to d xy orbitals. The bare values of the couplings are, however, different from the case considered in the previous section, because now intraorbital Hubbard interaction acts within the subset of the two hole pockets and within the subset of the two electron pockets. The corresponding interaction terms are the U 4 , U 5 ,Ū 5 ,Ũ 5 ,Ũ 5 terms in Eq. (7). The bare values of all these couplings are the Hubbard U; i.e., Interorbital Hubbard terms include the density-density interactions U 1 andŪ 1 between hole and electron pockets and the interactionŨ 4 within the d xz and d yz components of hole pockets. Because the d xy orbital interacts equally with the d xz and d yz orbitals, the bare values of U 1 andŪ 1 are equal; i.e., The exchange Hund interaction J acts in the subspace of d xz and d yz orbitals (Ũ 4 term) and between d xy and d xz =d yz orbitals (U 2 andŪ 2 terms). Again, the d xy orbital interacts equally with the d xz and d yz orbitals; hence, the bare values of U 2 andŪ 2 are equal: Finally, the pair-hopping interaction J 0 also acts in the subspace of the d xz and d yz orbitals (Ū 4 term) and between the d xy and d xz =d yz orbitals (U 3 andŪ 3 terms). Like for other interactions, the bare values of U 3 andŪ 3 are equal: Because the structure of the low-energy electronic states is the same as in the model that we consider in the previous section, the 14 RG equations are the same as in Eq. (10).
The couplingsŨ 4 andŨ 4 still flow to zero if the bareŨ 4 exceeds the bareŨ 4 , which is the case when U 0 > J. The couplingsŨ 5 andŨ 5 remain equal under RG, and both tend to zero when U > 0. One can further make sure that the couplings u 1 andū 1 , u 2 andū 2 , u 3 andū 3 , and u 5 andū 5 , which are equal at the bare level, remain equal under pRG. This reduces the set of RG equations to The transformation from U i to dimensionless u i is the same as before, and we remind the reader that we immediately find that the equation for u 4− decouples from the rest and u 4− tends to zero under pRG. The fixed trajectory for the other five equations is the same as for the model with d xz =d yz electron pockets; namely, , and These are the same expressions as in Eq. (11). In distinction from the previous case, however, now u i ¼ū i , i ¼ 1, 2, 3, already at the bare level. As the result, the system approaches the fixed trajectory faster than in the model that we studied in the previous section, i.e., at a larger running energy. We show the RG flow of the ratios of the couplings in Fig. 4.
The analysis of the vertices and the susceptibilities proceeds in the same way as in Sec. III. Again, the SDW susceptibility is initially the largest one, but its growth is halted by superconducting fluctuations, and eventually χ SDW does not diverge. The superconducting s þ− susceptibility diverges, but with the exponent α SC < 1. The d-wave Pomeranchuk susceptibility is initially the smallest one, but it increases under RG with the exponent α POM ¼ 1. As a result, the instability towards Pomeranchuk order becomes the leading one upon lowering the temperature.
The d-wave Pomeranchuk order has two coupled components. One is orbital order on the hole pockets, n xz − n yz , another is the difference between the fermionic densities on n xy orbitals between X and Y pockets, n X xy − n Y xy . The last term is not associated with the orbital order, as only the d xy orbital is involved, but it breaks C 4 symmetry, just like the n xz − n yz term. In a generic case, when electron pockets are composed of both d xy and d xz (d yz ) orbitals, one can expect that states near the electron pockets develop both an orbital order, i.e., a nonzero hn xz − n yz i, and an order with nonzero hn X xy − n Y xy i. The simultaneous presence of both order parameters should be visible at the M point in the 2-Fe zone [i.e., ðπ; πÞ]. Without Pomeranchuk order, there are two doubly degenerate fermionic states at the M point, even in the presence of spin-orbit coupling-one degeneracy is in the d xz =d yz subset, another is in the d xy subset (Ref. [54]). Below the Pomeranchuk transition both have to split, one due to nonzero hn xz − n yz i, another due to nonzero hn X xy − n Y xy i.

V. COMPARISON WITH EXPERIMENTS
Our results agree with five sets of experiments on FeSe. First, a recent ARPES study [55] has found shifts of the chemical potentials of the hole and electron pockets in opposite directions at T ≥ T s . This is fully consistent with the enhancement of the s þ− Pomeranchuk susceptibility found here. Second, orbital order, which breaks the symmetry between the d xz and d yz orbitals, has been detected in ARPES [53], and x-ray or neutron diffraction studies [56] have shown that it does not break the translational symmetry of the system. Such an orbital order has d-wave character and has zero transferred momentum, in agreement with our results. Third, ARPES studies of detwinned single crystals of FeSe [57] found strong evidence for a sign inversion of the d xz =d yz on-site energy splitting between the hole and electron Fermi pockets. This is fully consistent with the sign-changing d þ− order parameter we find here. We note that there is some controversy at the moment about the magnitude of this splitting on the electron pockets relative to hole pockets [53,58,59]. To fully understand this issue, one has to extend the analysis to T < T s , which is beyond the scope of this work. Fourth, recent neutron or x-ray diffraction measurements revealed that the phase diagram of FeSe under moderate pressures is very similar to those of Fe pnictides, with the nematic phase existing only very close to the stripe SDW transition [60]. This crossover from FeSe to Fepnictide behavior with pressure is consistent with our results, because at least one Fermi energy is expected to increase with pressure, which would cause the RG flow to stop at a L ¼ L F , where only the SDW and SC susceptibilities are sizable. And fifth, very recent ARPES measurements of the spectra near the M point (of the 2-Fe Brillouin zone) in FeSe below the nematic transition at 85 K (Ref. [61]) have detected the splitting in both d xz =d yz and d xy subsets. This is fully in line with the results that we discuss at the end of Sec. IV.

VI. OTHER THEORETICAL SCENARIOS
There have been other proposals for nematicity in FeSe. One early proposal [22] was that the origin of nematicity in FeSe is essentially the same as in Fe pnictides (Isingnematic spin order). This paper associated the large difference between the nematic and the magnetic transition temperatures in FeSe at ambient pressure (T s ¼ 85 K and T SDW ¼ 0) with stronger fluctuations of a continuous order parameter, also because of smaller E F of FeSe. The authors of Ref. [62] proposed a different magnetic scenario, attributing the smallness of T SDW in FeSe to the softening of magnetic fluctuations in the direction transverse to the typical SDW ordering vectors ðπ; 0Þ or ð0; πÞ. The proposal for such soft magnetic fluctuations was motivated in Ref. [62] by first-principles calculations that showed that in FeSe several magnetic states with different ordering vectors have near-equal energy. The authors argued that this softening strongly reduces T SDW and also enhances T s .
Other magnetic scenarios for the nematic phase in FeSe were proposed in Refs. [63,64], in which the authors associated nematic order below T s with a more complex type of long-range magnetic order-antiferroquadrupolar order in Ref. [63] and ferroquadrupolar order in Ref. [64]. In yet another magnetic scenario, the authors of Ref. [65] used a local spin model and suggested that the nematic phase in FeSe is a nematic quantum paramagnetic phase of spin S ¼ 1 local moments with strongly frustrated exchange interactions-similar, but not equivalent, to the dimer phase of an S ¼ 1=2 antiferromagnet with frustrated interactions.
We do not believe that current experiments on FeSe have ruled out any of these magnetic scenarios. The recent pressure experiments of Ref. [60], however, offer important benchmarks to test these different scenarios. In our case, the quick suppression of T s observed under pressure, combined with the sudden appearance of stripe magnetism, is explained by an increase in the Fermi energy E F with pressure, signaling a crossover from FeSe to Fe-pnictide behavior. If quantum oscillation measurements under pressure can indeed confirm the predicted behavior of E F , that would lend strong support to our proposal for a spontaneous orbital order.
The possibility of spin-fluctuation-driven orbital order has also been recently discussed in Ref. [9]. These authors considered a large ratio of intraorbital Coulomb interaction U and Hund's interaction J, in which case the bare interaction in the orbital channel is attractive. They then analyzed how this attraction is enhanced by spin fluctuations by focusing on a particular Aslamazov-Larkin diagram. This diagram is one of many diagrams that give rise to the RG flow that we study in this paper. Our results are consistent with those of Ref. [9] (modulo the fact that we include many other terms besides the Aslamazov-Larkin diagram), but our analysis goes beyond this work in two important aspects: (i) we argue that at low energies the interaction in the orbital channel gets attractive even if the bare interaction was repulsive (as is the case if U=J is not too large), and (ii) we argue that the orbital susceptibility diverges not only faster than the susceptibility in the SDW channel (as the authors of Ref. [9] argued), but also faster than the susceptibility in the s þ− superconducting channel, which the authors of Ref. [9] did not consider. The latter is important for the elucidation of why orbital order develops prior to superconductivity. Nematicity due to orbital order was recently discussed in Ref. [66]. The authors of this work, however, assumed that an orbital order develops and studied how it is affected by impurities.

VII. SUMMARY
In this paper, we employ the analytical RG technique to analyze the interplay between SDW, SC, and orbital d-POM order in Fe-based superconducting materials. That magnetic fluctuations promote attraction in the s þ− SC channel is well known, and the RG analysis indeed confirms this. However, we find that the same magnetic fluctuations also promote attraction in the C 4 -symmetrybreaking d-POM orbital channel, even if the bare interaction in the d-POM channel is repulsive. We go beyond the analysis of the running couplings and compute running susceptibilities in SDW, SC, and orbital channels. We find that within the energy range where the RG approach is valid, the susceptibility towards a spontaneous orbital order increases in the process of RG flow and diverges with the largest exponent. The SC susceptibility also diverges, but with a smaller exponent, and SDW susceptibility does not diverge. As a consequence, if the system reaches the scale of the leading instability within the applicability range of parquet RG, the leading instability is towards orbital order, the subleading is towards SC, and SDW order does not develop. This hierarchy of instabilities agrees with the one observed in bulk FeSe at ambient pressure, and we conjecture that nematic order without SDW magnetism, observed in FeSe, is a spontaneous orbital order driven by magnetic fluctuations. We argue that the orbital scenario agrees with five sets of experiments on FeSe.
If the order does not develop down to the lower boundary of applicability of RG (which is set at an energy of order E F ), the leading instability (at a smaller energy or temperature) will be in the channel in which the susceptibility is the largest at this lower boundary. Because the susceptibility in the d-POM channel becomes the largest only near the instability, d-POM order does not develop in this situation. Instead, the competition becomes between SDW and SC, with SC winning at higher doping and SDW winning at smaller doping. In this situation, nematic order is most likely caused by composite stripe SDW fluctuations, which split the temperatures at which a discrete C 4 and a continuous O(3) symmetry get broken. This likely happens in LaFeAsO, BaFe 2 As 2 , and NaFeAs. From this perspective, our work provides a unified microscopic description of the behavior of different families of FeSCs.

APPENDIX A: DIMENSIONLESS INTERACTIONS
In this Appendix, we list the expressions of the dimensionless interactions in terms of the original interaction constants, introduced in Eq. (8). For m c ¼ m d ¼ m h and m x ¼ m y ¼ m e that we consider in the main text, we have The tetragonal symmetry allows us to decompose the running interactions in Eq. (8) into different channels. To achieve this goal we construct bilinear fermion operators that transform irreducibly under the symmetry group of the lattice. We consider separately the bilinear combinations in the particle-hole channel at zero momentum and at momenta Q 1;2 , and the bilinear combinations in the particle-particle channel at zero total momentum.
1. Bilinear fermion combinations in the charge and spin particle-hole channels at large momentum transfer The two possible order parameters that describe chargedensity wave order with momenta ðπ; 0Þ and ð0; πÞ are Another two possible order parameters with large momentum transfer describe anti-ferro-orbital order. The corresponding order parameters arē These order parameters differ from the ones in Eq. (B1) because they are off diagonal in the orbital index. The four possible SDW order parameters with the same momenta are The components of Eq. (8) that describe the interactions in CDW and SDW channels are and The bare values of the couplings in Eqs. (B5) and (B6) are 2. Bilinear fermion combinations in the particle-particle channel We focus on the singlet pairing with zero total momentum. These bilinear combinations form reducible representations of the D 4h group, separately for electrons f and holes d. Accordingly, we introduce the notations Spin-singlet bilinear combinations are even under inversion, and form one-dimensional irreducible presentations of the D 4h group: A 1g , A 2g , B 1g , and B 2g . We have The spin-singlet A 2g combination κ fðdÞ A 2 actually vanishes, as κ μμ 0 is invariant under simultaneous interchange of orbital indices and spin projections (the spin-triplet A 2g combination does not vanish).
The interaction component in the Cooper channel is obtained by setting k 1 ¼ −k 2 in Eq. (8). Expressing Eq. (8) in terms of the combinations Eq. (B9), we obtain The bare values of the couplings in Eqs. (B11)-(B13) are 3. Bilinear fermion combinations in particle-hole charge channel with zero momentum transfer The bilinear combinations of fermions with zero momentum transfer in the particle-hole charge channel are Again, these combinations form reducible representations of the D 4h group, separately for electrons f and holes d (Ref. [44] To obtain the interactions in the particle-hole charge channel at zero momentum transfer [the ones that renormalize bilinear combinations in Eq. (B16)], we set k 1 ¼ k 2 or k 1 ¼ k 4 in Eq. (B16). Expressing Eq. (8) in terms of the combinations Eq. (8), we obtain where where the contributions of direct and exchange processes in Eq. (8) are taken into account. We further have The bare values of the couplings in Eqs. (B18)-(B21) are a. Pomeranchuk susceptibilities We focus first on A 1g and B 1g channels, which we abbreviate as just as A 1 and B 1 . The equation for the vertex function reads where the polarization matrix has a diagonal element given by the static autocorrelation functions of the bilinear operators ρ f A 1g ;B 1g and ρ d A 1g ;B 1g , respectively, defined in Eq. (B16). We also have from Eqs. (B18) and (B19): The instability occurs when one of the eigenvalues of the matrix These expressions simplify on fixed trajectories defined by Eq. (11). And we consider now this case for the two channels separately. For the A 1g channel, we have from Eqs. (B27) and (11) For the B 1g channel, we similarly have from Eqs. (B27) and (11) The couplings λ þ− A 1 ¼ 2ð4u 1 C þ u 1 jγ 4 jÞ ¼ 2ð2u 1 C − u 4 Þ and λ þ− B1 ¼ 2u 1 jγ 4 j ¼ −2u 4 are dimensionless analogs of 4U 1 − U 4 and −U 4 in Table I. As we discuss in the main text, the equality Eq. (B29) holds strictly only on the fixed trajectory. In general, λ þ− B 1 > λ þþ B 1 , and the Pomeranchuk instability occurs in the d þ− symmetry channel. The scaling of susceptibilities is We now turn to the A 2 and B 2 channels. In this case, within the RPA the instability occurs when one of the two conditions is satisfied, 2ðu 5 AE 2ũ 5 ∓ũ 5 Þ ¼ 1; 2ðu 4 AE 2ũ 4 ∓ũ 4 Þ ¼ 1; ðB31Þ where the upper (lower) sign refers to A 2g and B 2g channels, respectively.
If we use the values of the couplings along the fixed trajectory, we find that the susceptibilities in these two channels occur at the same L as the one in the B 1 channel. However, once we include the corrections, we find that the instability in the B 1 channel occurs prior to the ones in the A 2 and B 2 channels.

APPENDIX C: RG FLOW OF SUSCEPTIBILITIES FOR DIFFERENT BARE VALUES OF THE COUPLINGS
In this Appendix, we present the results of the calculations of the flow of the susceptibilities in SDW, SC (s þ− ), and the B 1 Pomeranchuk channel for different values of the bare parameters. The goal here is to show that the system behavior at low energies is universal and is governed by the flow of the couplings towards the stable fixed trajectory, but the value of the RG scale L, at which the susceptibility in the Pomeranchuk channel becomes the largest, depends on the bare values of the couplings.
To demonstrate this, we obtain the flow of the couplings numerically, and then solve the RG equations for the vertices using the running couplings as inputs. We do not assume in this calculation that the system is already near the fixed trajectory.
We present the results in Figs. 5 and 6. The first is for the case when the bare interaction in the Pomeranchuk channel is repulsive, the second when it is attractive. We see that the flow of the susceptibilities is virtually the same in both cases, but for attractive bare interaction, the Pomeranchuk susceptibility becomes the largest at a smaller L. Note that for the parameters we have chosen, this happens before the RG flow reaches the fixed trajectory, as evidenced by the fact that the susceptibility in the SDW channel is still larger than the one in the s þ− superconducting channel.  6. The flow of the susceptibilities for attactive bare interaction in the orbital channel. We used J ¼ 0.02, U ¼ ð6.1ÞJ, U 0 ¼ U − 2J, J 0 ¼ J. For these parameters, the bare interaction in the B 1 Pomeranchuk channel is attractive [it is 2U 0 − ðU þ JÞ, Ref. [9]]. We see that the flow of the susceptibilities is very similar to the case when the bare B 1 interaction is repulsive, but the Pomeranchuk susceptibility diverges at a smaller L.