Nematicity and magnetism in FeSe and other families of Fe-based superconductors

Nematicity and magnetism are two key features in Fe-based superconductors, and their interplay is one of the most important unsolved problems. In FeSe, the magnetic order is absent below the structural transition temperature $T_{str}=90$K, in stark contrast that the magnetism emerges slightly below $T_{str}$ in other families. To understand such amazing material dependence, we investigate the spin-fluctuation-mediated orbital order ($n_{xz}\neq n_{yz}$) by focusing on the orbital-spin interplay driven by the strong-coupling effect, called the vertex correction. This orbital-spin interplay is very strong in FeSe because of the small ratio between the Hund's and Coulomb interactions ($\bar{J}/\bar{U}$) and large $d_{xz},d_{yz}$-orbitals weight at the Fermi level. For this reason, in the FeSe model, the orbital order is established irrespective that the spin fluctuations are very weak, so the magnetism is absent below $T_{str}$. In contrast, in the LaFeAsO model, the magnetic order appears just below $T_{str}$ both experimentally and theoretically. Thus, the orbital-spin interplay due to the vertex correction is the key ingredient in understanding the rich phase diagram with nematicity and magnetism in Fe-based superconductors in a unified way.


I. INTRODUCTION
In Fe-based superconductors, the origin of the electronic nematic state and its relation to the magnetism have been a central unsolved problem. Recently, the nonmagnetic nematic state in FeSe has attracted increasing attention as a key to solve the origin of the nematicity. FeSe undergoes a structural and superconducting transitions at T str = 90K and T c = 9K, respectively, whereas the magnetic transition is absent down to 0 K [1]. The strength of the low-energy antiferro-magnetic (AFM) fluctuations is very weak above T str , while it starts to increase below T str [2][3][4][5][6][7]. In stark contrast, the magnetic transition occurs at T mag slightly below T str in other undoped Fe-based superconductors. Since the relation T str > T mag is unable to be explained by the random-phase-approximation (RPA), we should develop the microscopic theory beyond the mean-field-level approximations.
Up to now, two promising triggers for the structure transition have been discussed intensively: In the spinnematic scenario [8][9][10][11][12], the trigger is the spin-nematic order. This spin-fluctuation induced spin-quadrupole order could emerge above T mag in highly magnetically frustrated systems. In the orbital order scenario [13][14][15][16], the trigger is the ferro-orbital (FO) order n xz = n yz . Above T str , the strong orbital or spin-nematic fluctuations are observed by the measurements of shear modulus C 66 [2,17,18], Raman spectroscopy [19][20][21][22], and in-plane resistivity anisotropy [23,24]. The nematic orbital fluctuations originate from the strong orbital-spin mode-coupling due to the strong-coupling effect, which is described by the Aslamazov-Larkin vertex correction (AL-VC). The electronic nematic state studied in singleorbital models [25] is more easily realized in multiorbital systems thanks to the AL-VC mechanism [16]. Except for the presence or absence of magnetism below T str , FeSe and other Fe-based superconductors show common electronic properties. Below T str , in both FeSe and BaFe 2 As 2 , large orbital polarization ∆E ≡ E yz − E xz ∼ 50 meV [26][27][28][29][30][31][32][33] is observed. Such large ∆E originates from the electron-electron correlation since the lattice distortion (a − b)/(a + b) is just 0.2 ∼ 0.3%, as we discuss based on band calculation in Appendix A. Above T str , the electronic nematic susceptibility is enhanced in both BaFe 2 As 2 [17,19,23] and FeSe [2,18], following the similar Curie-Weiss behavior. These facts indicate that the common microscopic mechanism drives the nematic order and fluctuations in all Fe-based superconductors, in spite of the presence or absence of the magnetism.
The realistic multiorbital Hubbard models for Fe-based superconductors, which are indispensable for the present study, were derived by using the first-principles method in Ref. [34]. To understand the absence of the magnetism below T str in FeSe, one significant hint is the smallness of the ratio between the Hund's and Coulomb interactions, J/Ū , since the Hund's coupling enlarges (suppresses) the intra-site magnetic (orbital) polarization, which is verified by the functional renormalization-group (fRG) theory [35,36]. Another significant hint is the absence of the d xy -orbital hole-pocket in FeSe, which is favorable for the orbital-spin interplay on the (d xz , d yz )-orbitals due to the AL-VC mechanism.
The goal of this paper is to explain the amazing variety of the electronic nematic states in Fe-based superconductors, especially the non-magnetic nematic state in FeSe, on the same footing microscopically. For this purpose, we study the spin-fluctuation-mediated orbital order by applying the self-consistent vertex-correction (SC-VC) method [16] to the first-principles models. In FeSe, the orbital-spin interplay is significant because of the smallness ofJ/Ū and the absence of d xy -hole pocket. For this reason, the orbital order is realized even when the spin fluctuations are substantially weak. The rich variety of the phase diagrams in Fe-based superconductors, such as the presence or absence of the magnetic order in the nematic phase, are well understood by analyzing the vertex correction seriously. The SC-VC theory had been successfully applied to explain the phase diagram in LaFeAs(O,H) [37], nematic CDW in cuprates [38,39], and triplet superconductivity in Sr 2 RuO 4 [36].
We comment that the localized spin models have been successfully applied to the nematic order, stripe magnetic order, and so on [40]. On the other hand, weak-coupling theories have also been applied to Fe-based superconductors satisfactorily [41]. In the present study, we study the mechanisms of the nematicity and magnetism in various Fe-based superconductors in terms of the itinerant picture, by taking the strong-coupling effect due to the AL-VC into account. The significant role of the AL-VC on the orbital fluctuations has been confirmed by the fRG theory [35,36]. The AL-VC is important to reproduce the Kugel-Khomskii-type orbital-spin interaction [16].

II. MODEL HAMILTONIAN AND SC-VC THEORY
In the present study, we study the realistic d-p Hubbard models for M=LaFeAsO and FeSe by applying the SC-VC method [16]. In Eq. (1), is the 8-orbital d-p tight-binding (TB) model in k-space, which is obtained by using the WIEN2k and WAN-NIER90 softwares; see Appendix A for detailed explanation. σ is the spin index, and l, m are the orbital indices: Hereafter, we denote the five d-orbitals as d 3z  [42,43] studies. To reproduce experimental bandstructure of FeSe, we introduce the additional intra-orbital hopping parameters into H 0 FeSe , in order to shift the d xyorbital band [d xz/yz -orbital band] at (Γ, M, X) points by (0, −0.25, +0.24) [(−0.24, 0, +0.12)] in unit eV; see Appendix A. These energy shifts might be induced by the self-energy [44]. The constructed FSs in the FeSe model is shown in Fig. 1 (c). Since each Fermi pocket is very shallow, the superconductivity in FeSe could be close to a BCS-BEC crossover [45].
In Eq. (1), H U M is the first-principles screened Coulomb potential for d-orbitals given by the "constrained-RPA method" [34] given as where U m,l and J m,l are orbital-dependent Coulomb and Hund's interactions for d-electrons, respectively [34]. The averaged intra-orbital Coulomb interactionŪ ≡  [34]. Thus, the ratiō J/Ū = 0.0945 in FeSe is much smaller than the ratiō J/Ū = 0.134 in LaFeAsO. Such strong material dependence of (Ū ,J) is understood as follows: U l,m is strongly screened by the screening bands (excluding the 8 bands in H 0 M ) whereas the screening of J l,m is much weak, and the number of the screening bands is small in FeSe [46]. The factor r(< 1) in Eq. (1) is introduced to adjust the spin fluctuation strength. The ratio J l,m /U l,m is unchanged by introducing the factor r [47,48]. The 8 × 8 Green function in the orbital basis is given asĜ where k = (k, ǫ n = (2n + 1)πT ),ĥ 0 M (k) is the kinetic term in Eq. (2), andẑ −1 ≡ 1 − ∂Σ/∂ǫ| ǫ=0 represents the mass-enhancement due to the self-energy at the Fermi level. Here, we introduce the constant mass-enhancement factor for d-orbital 1/z l (≥ 1). Then, Eq. (4) gives the coherent part of the Green function, which mainly determines the low-energy electronic properties. In the present study, r and z l are the fitting parameters. In FeSe, the orbital order is obtained in the real first-principles Hamiltonian (r ≈ 1) by taking the experimental massenhancement factors z −1 l ≈ 4 into account, as shown later.
The d-orbital charge (spin) susceptibilities (per spin) is given in the following 5 2 × 5 2 matrix form: whereΦ c(s) (q) =χ 0 (q) +X c(s) (q) is the irreducible susceptibility for the charge (spin) channel. In the SC-VC theory, we employ the AL-VC asX c,s (q), and perform the self-consistent calculation with respect to the AL-VC and susceptibilities. Using the Green function in Eq. (4), the bare susceptibility is where q = (q, ω l = 2lπT ). Also, the AL-VC for the charge susceptibility is given as where p = (p, ω m ), andV s,c (q) ≡Γ s,c +Γ s,cχs,c (q)Γ s,c . The three-point vertexΛ(q; p), which gives the coupling between two-magnon and one-orbiton, is given as Λ l,l ′ ;a,b;e,f (q; p) and Λ ′ m,m ′ ;c,d;g,h (q; p) ≡ Λ c,h;m,g;d,m ′ (q; p) + Λ g,d;m,c;h,m ′ (q; −p − q). We stress that the strong temperature dependence of the three-point vertex is significant for realizing the orbital order. We include all U 2 -terms without the double counting in order to obtain quantitatively reliable results. Equation (7) means that the charge AL-VC becomes large in the presence of strong spin fluctuations. More detailed explanations are presented in the textbook [49] In the present study, we neglected the spin-channel VCs since it is expected to be unimportant as discussed in Ref. [50]. In Appendix B, we verify validity of this simplification in the present model by performing a time-consuming self-consistent calculation with respect to both charge-and spin-channel Maki-Thompson (MT) and AL-VCs.
As explained in Ref. [50], the development of χ c x 2 −y 2 (0) is mainly induced by the diagonal elements ofΦ c with respect to l = 2, 3. If we drop the off-diagonal elements ofΦ c approximately, χ c x 2 −y 2 (0) is given as where U ≡ U 2,2 = U 3,3 , J ≡ J 2,3 , Φ c ≡ χ 0 l,l;l,l (0) + X c l,l;l,l (0) (l = 2 or 3). Thus, the charge Stoner factor is , only the spin fluctuations develop since the relation α S > α C is satisfied for J > 0. However, the opposite relation α C > α S is realized if the relationΦ c ≫ χ 0 is satisfied due to the charge-channel AL-VC [37].

III. NUMERICAL RESULTS FOR LAFEASO AND FESE
First, we analyze the LaFeAsO model based on the SC-VC theory. For z = 1 for each l, the obtained χ s (q) and χ c x 2 −y 2 (q) are shown in Fig. 2 (a) and (b), respectively, for r = 0.41 (Ū = 1.74 eV) at T = 50 meV. Here, the number of k-meshes is 32 × 32, and the number of Matsubara frequencies is 256. Thus, both AFM and FO susceptibilities develop divergently, and the realized enhancement factors are S S ≈ 40 and S C ≈ 50. The rdependences of the enhancement factors at T = 50 meV are shown in the inset of Fig. 2 (c): Both S S and S C increase with r, and they are equivalent at r * = 0.41. The lower the temperature is, the smaller r * is, whereas the value of S S = S C at r * is approximately independent of T . Similar result is obtained in BaFe 2 As 2 model as shown in Appendix C.
The orbital-spin interplay due to the AL-VC is intuitively understood in terms of the strong-coupling picture U ≫ W band [37]: As shown in Fig. 2 (d), when the FO order n xz ≫ n yz is realized, the nearest-neighbor exchange interaction has large anisotropy J Then, the stripe AFM order with Q = (π, 0) appears if J (2) is not too small. Thus, the FO order/fluctuations and AFM order/fluctuations emerge cooperatively in the localized model, and such Kugel-Khomskii-type orbital-spin interplay is explained by the AL-VC in the weak-coupling picture.
Next, we analyze the FeSe model, in which the ratiō J/Ū is considerably small. In FeSe, the experimental mass-enhancement factor is ∼ 10 for d xy -orbital, and 3 ∼ 4 for other d-orbitals according to the ARPES study [27]. Therefore, we put z −1 l = z −1 for l = 4 and z −1 4 = 3z −1 in the present study. We find that the peak of χ s (q) moves from q = (π, π) to the experimental peak position q = (π, 0) [4][5][6][7] for z −1 4 ≥ 1.5z −1 , and the results of the SC-VC method are essentially unchanged for z −1 4 ≥ 1.5z −1 . Figures 3 (a) and (b) show the obtained χ s (q) and χ c x 2 −y 2 (q) for r = 0.25 (Ū = 1.76 eV) at T = 50 meV in the case of z = 1. We see that only the FO susceptibility develop divergently [S C ≈ 50], whereas the AFM susceptibility remains small [S S ≈ 8], consistently with experiments for FeSe. The r-dependences of the Stoner enhancement factors at T = 50 meV are shown in the inset of Fig. 3 (c): With increasing r, S C increases rapidly whereas S S remains small. Figure 3 (c) shows the temperature dependences of the enhancement factors at r = 0.25. We stress that S C approximately follow the Curie-Weiss behavior with the Weiss temperature θ C = 48 meV, which is consistent with the experimental Curie-Weiss behavior with positive θ C in FeSe [2]. Since the spin Weiss temperature takes a large negative value (θ S ∼ −20 meV), which is also consistent with experiments, one may consider that the orbital order in FeSe stems from causes other than spin fluctuations.

IV. ORIGIN OF THE RELATION SC ≫ SS IN FESE
In this section, we discuss why the relation S C ≫ S S (θ C > 0 and θ S < 0) is realized in FeSe. First, we focus on the ratio between the Hund's and Coulomb interac-tionsJ/Ū . It is intuitively obvious that the ratioJ/Ū is an important control parameter for the orbital nematicity: For largerJ/Ū , the local configuration of the twoelectrons in the (d xz , d yz )-orbitals is |d xz , ↑ ⊗ |d yz , ↑ , where the magnetic moment is s z = 1 whereas the orbital polarization is n xz − n yz = 0. Thus, the smallness ofJ/Ū in FeSe is favorable for the emergence of the orbital order without magnetization.
Microscopically, as we discuss in Sec. II, the charge Stoner factor for χ c , where X c is the charge AL-VC for orbital 2 or 3 at q = 0. SinceJ/Ū = 0.0945 in FeSe, the orbital order is realized by relatively small AL-VC; X c ∼ 0.9χ 0 (0). In LaFeAsO, in contrast, large AL-VC of order ∼ 2χ 0 (0) is required to realize the orbital order. The obtained AL-VCs in both systems are shown in Fig. 9 (c) in Appendix D.
We discuss why the AL-VC is important in the FeSe model with θ S < 0: As we explain in Appendix D. the T -dependence of the AL-VC is given as X c ∼ Λ 2 T S S , where Λ is the three-point vertex that represents the interference between two short-living magnons. We find the relation Λ 2 ∝ T −a with a ≈ 1 at low temperatures due to the good nesting between h-FSs and e-FSs [20]. Thanks to the strong enhancement of Λ at low temperatures, the orbital order (α C = 1) is realized even if θ S is negative. (Note that T S S decreases as T → 0 when θ S < 0.) Thus, serious diagrammatic analysis of the AL-VC is necessary to understand the rich normal-state phase diagrams in Fe-bases superconductors.
The enhancement of the nematic susceptibility due to the significant T -dependence of Λ 2 (∝ T −a ) had been discussed in Refs. [20,21,51,52]. However, the reported exponent a is not universal, since it depends on the bandstructure and temperature range. In Appendix E, we show the T -dependence of the three-point vertex for LaFeAsO and FeSe models for wide temperature range. It is found that a ≈ 1 for T = 20z ∼ 100z[meV], where z < 1 is the band-renormalization factor. Due to such large T -dependence of a, χ c x 2 −y 2 (0) obtained by the present study follows the Curie-Weiss law only approximately.
Finally, we stress the importance of the orbital dependence of the spin fluctuation strength. Since the is enlarged by the spin fluctuations on the (d xz , d yz )-orbitals. More detailed analysis is given in Appendix D.

V. EFFECT OF THE MASS-ENHANCEMENT FACTOR
Here, we study the effect of the mass-enhancement factor: We study the FeSe model in the case of z −1 = 4 for (d xz , d yz )-orbitals. The obtained S C,S as functions of r are shown in the inset of Fig. 3 (d) at T = 12.5 meV.
Here, S S remains small even for r ∼ 1 since the bare susceptibility is suppressed by z. In contrast, S C is enlarged to 50 for r ≈ 0.97, which is very close to the exact first-principles Hubbard model H FeSe (r = 1). The Tdependences of S C,S are shown in Fig. 3
To understand the similarity between the results in Fig. 3 (c) for z = 1 and the results in Fig. 3 (d) for z −1 = 4, we prove that both α C and α S are independent of z under the rescaling T → zT and (U, J) → (U, J)/z. Here, we assume that z −1 l = z −1 , and neglect the Tdependence of µ for simplicity. Under the scaling T → zT , the Green functionĜ(k, n) at Matsubara integer n given in Eq. (4) is independent of z. For this reason, the bare susceptibility χ 0 (q) = −T k,n G(k + q, n)G(k, n) is proportional to z. By following the same procedure, the three-point vertex Λ is scaled by z, and therefore the AL-VC X c (0) ∼ T U 4 q Λ(0; q) 2 χ s (q) 2 is proportional to z under the scaling T → zT and (U, J) → (U, J)/z. Thus, both spin and charge irreducible susceptibilities are proportional to z, and both α S and α C are unchanged under the rescaling T → zT and (U, J) → (U, J)/z. That is, the Weiss temperatures θ S(C) are scaled by z. The validity of these scaling relations are confirmed by the numerical study in Fig. 3.
It is possible to obtain z −1 by calculating the self-energyΣ(k) together withχ s,c (q) self-consistently. In this case, fine tuning of r will be unnecessary since the  relation α S,C < 1 is assured ifΣ(k) andχ s,c (q) are calculated self-consistently in two-dimensional systems (Mermin-Wagner theorem) [53]. This is our important future issue. Here, we study the electronic states in the FO order n xz = n yz established below the structure transition temperature T str , at which the shear modulus C 66 reaches zero. According to the linear-response theory, the electronic orbital susceptibility given by the SC-VC theory, and g is the phonon-mediated Jahn-Teller energy [54]. Therefore, C 66 ∝ (T − T str )/(T − θ C ), and T str = θ C + g is slightly higher than θ C due to the weak electron-phonon coupling (g ≈ 10 ∼ 50 K) [2,17,18]. Figure 4 (a) shows the T -dependence of S S given by the RPA for LaFeAsO and FeSe for z = 1. Here, we introduce the orbital polarization −∆E/2 (∆E/2) for the d xz(yz) -level. We put S S = 20 (5) for LaFeAsO (FeSe) at T str = 50 meV, and assume a mean-field-type T -dependence; ∆E = ∆E 0 tanh(1.74 T str /T − 1) with ∆E 0 = 80 meV. (For z −1 = 4, the renormalized orbital polarization z∆E 0 is just 20 meV.) In both LaFeAsO and FeSe, S S are enhanced by ∆E, since α S increases linearly with ∆E at q = (π, 0) as discussed in Ref. [54]. In LaFeAsO, the magnetic order temperature T mag increases from θ S to just below T str since S S is already large at T str . In contrast, in FeSe, the enhancement of χ s (π, 0) is much moderate [55].
We also perform the self-consistent analysis of the orbital-polarization (∆E xz (k), ∆E yz (k)) and anisotropic χ s (q), which is a natural extension of the SC-VC theory into the orbital ordered state [56]. The obtained S S and k-dependent orbital polarization are shown in Figs. 4 (a) (inset) and (b), respectively. The parameters are r = 0.256 and 1/z 4 = 1.6. The difference ∆n = n xz − n yz is 0.2%. The hole-pocket around Γ-point becomes ellipsoidal along the k y -axis due to the "sign-reversing orbital polarization", in which ∆E xz (0, k) − ∆E yz (k, 0) shows the sign reversal as shown in Fig. 4 (c). Due to this sign reversal, S S in the inset of Fig. 4 (a) tends to saturate below 40 meV [33]. Also, two Dirac-cone FSs appear around X-point when ∆E yz (π, 0) > 50 meV. These results are essentially consistent with the recent ARPES studies reported in Refs. [27][28][29][30][31][32][33]. The obtained orbitalpolarization (∆E xz (k), ∆E yz (k)) belongs to B 1g representation, and therefore it is consistent with the "d-wave orbital order" discovered in Ref. [31]. The d-wave orbital order is theoretically obtained by the mean-field approximation by introducing phenomenological long-range interaction [57], whose microscopic origin might be the AL-VC studied in this paper.
According to Ref. [4], χ s (q, ω) shows the broad maximum at q = (π, 0) at low-energies (ω 10 meV), and its strength is almost independent for T > T str . The magnitude of the low-energy spin susceptibility in FeSe  . The FO order is introduced below Tstr = 50 meV. Inset: T -dependences of SS for the FeSe model obtained by calculating the k-dependent orbital polarization and χ s (q) self-consistently. SS tends to saturate below 40 meV due to the sign-reversing orbital polarization. (b) Self-consistent solution of the orbital polarization (∆Exz(k), ∆Eyz(k)) in the orbital ordered state in the FeSe model at T = 50 meV. The shape of the C2-symmetric FSs in (b) is consistent with the experimental reports [27][28][29][30][31][32][33]. We also show (c) the ∆E xz(yz) (k) along the k y(x) -axis, and (d) the C2-symmetric χ s (q) in the orbital-ordered state.
is one order of magnitude smaller than that in BaFe 2 As 2 [58], whereas its magnitude would be comparable to that in LiFeAs [59]. This experimental report in FeSe will be consistent with the present theoretical result with the moderate S S ∼ 10 in Figs. 3 (c) and (d). Note that experimental dispersion relation in χ s (q, ω) for ω 100 meV is qualitatively understood based on the present FeSe model by considering the band-renormalization factor [7].
ofJ/Ū . Numerical study for NaFeAs and BaFe 2 As 2 are presented in Appendix C. In NaFeAs and FeSe, in which J/Ū is smaller, the obtained θ S /θ C decreases to 0.4 and −0.4, respectively. In Fig. 5 (a), experimental values of T mag /T str and θ NMR /T str are also shown, where θ NMR is the Weiss temperature of 1/T 1 T above T str . Since T str = θ C + g (g ≈ 10 ∼ 50 K) and θ NMR = θ S , the relation θ NMR /T str θ S /θ C is expected theoretically. In addition, the relation T mag /T str θ S /θ C is expected since T mag is substantially higher than θ S in the FO ordered state. These two theoretically predicted relations are verified in Fig. 5 (a). Thus, the ratio θ S /θ C is well scaled by the parameterJ/Ū , consistently with the discussion in Sec. IV. Figure 5 (b) shows the critical value of the spin Stoner factor for α C ≈ 1 in each model, α cr S . It is found that α cr S increases withJ/Ū qualitatively. In addition, we plot α cr S for the FeSe (LaFeAsO) TB model with different Coulomb interactions: H 0 FeSe(LaFeAsO) + rH U M . In both FeSe and LaFeAsO TB models, α cr S monotonically increases withJ/Ū , whereas α cr S is clearly small for FeSe TB model. There are two reasons why α cr S is smaller for the FeSe bandstructure. One reason is the absence of the d xy -orbital h-FS in FeSe: As we discussed in Sec. IV, the d xy -orbital spin fluctuations are unnecessary for the development of χ c x 2 −y 2 (0) due to the AL-VC. Another reason is the smallness of the FSs in FeSe: We found numerically that α cr S decreases when the size of the FSs is smaller, since the three-point vertex Λ m ≡ δχ 0 m (q)/δ∆E m , which is odd with respect to G, increases in magnitude when the particle-hole asymmetry is large: In fact, we analyzed the undoped LaFeAsO model with tiny FS pockets by introducing the positive/negative potentials around the electron/hole FSs, and verified that the orbital order is realized by small α S . Recently, the advantage of the small FSs for the nematicity had been stressed by the renormalization group study in Ref. [63].

C. Summary
The emergence of the electronic nematic order has attracted increasing attention as a fundamental phenomenon in strongly correlated metals. In this paper, we studied the origin of the nematicity in Fe-based superconductors, by paying the special attention to the nonmagnetic nematic order in FeSe. By applying the orbital+spin fluctuation theory to the first-principles dp Hubbard models, we succeeded in explaining the rich variety of the phase diagrams in Fe-based superconductors, such as the nonmagnetic/magnetic nematic order in FeSe/LaFeAsO. The key model parameter to realize rich phase diagram is J/U ; the ratio between the Hund's and Coulomb interactions. In addition, the ratio θ S /θ C tends to decrease as the size of the FSs shrinks, as discussed in Sec. VI B.
In both FeSe and LaFeAsO, strong orbital susceptibility χ c x 2 −y 2 (0) ∝ (T − θ C ) −1 with positive θ C is realized by the strong orbital-spin interplay due to the strong-coupling effect, called the Aslamazov-Larkin vertex correction in the field theory. In the FeSe model, ferro-orbital order is established even when the spin Weiss temperature θ S is negative as shown in Fig. 3, since the three-point vertex (=the coupling between twomagnon and one-orbiton) increases at low temperatures as Λ ∝ T −0.5 . In contrast, the spin-nematic susceptibility driven by the spin susceptibility should be T -independent if θ S < 0, as discussed in Ref. [7]. Therefore, we conclude that the nematicity in FeSe originates from the orbital order/fluctuations.
The nematic orbital fluctuations might play important roles in the pairing mechanism in Fe-based superconductors [64]. In FeSe, T c increases from 9 K to 40 K under pressure, accompanied by the enhancement of spin fluctuations [1]. At the same time, the system approaches to the orbital critical point since T str decreases to zero under pressure. These facts indicate the important role of the spin+orbital fluctuations in FeSe.

Acknowledgments
We are grateful to A. Chubukov, P.J. Hirschfeld, R. Fernandes, J. Schmalian, Y. Matsuda, T. Shibauchi and T. Shimojima for useful discussions. This study has been supported by Grants-in-Aid for Scientific Research from MEXT of Japan.
Appendix A: Eight-orbital models for FeSe and LaFeAsO Here, we introduce the eight-orbital d-p models for FeSe and LaFeAsO analyzed in the main text. We first derived the first principles tight-binding models using the WIEN2k and WANNIER90 codes. Crystal structure parameters of FeSe and LaFeAsO are given in Refs. [65] and [66], respectively. The obtained bandstructure and FSs in the LaFeAsO model are shown in Fig. 1 in the main text. In deriving the FeSe model, we introduce the k-dependent shifts for orbital l, δE l , in order to obtain the experimentally observed FSs. In FeSe, we introduce the intra-orbital hopping parameters into We also explain the orbital-dependent Coulomb interaction. The bare Coulomb interaction for the spin channel in the main text is Also, the bare Coulomb interaction for the charge channel is otherwise. (A2) Here, U l,l , U ′ l,l ′ and J l,l ′ are the first principles Coulomb interaction terms given in Ref. [34].
Finally, we perform the band calculations for the orthorhombic phase of FeSe and LaFeAsO, based on the experimental crystal structures. In both compounds, the obtained band splitting is too small to explain the large orbital polarization (∼ 60 meV) observed by ARPES studies. This result means that the orbital order originate from the electron-electron correlation, which is not included in the band calculation. Figure 6 (a) is the non-magnetic bandstructure in the orthorhombic LaFeAsO obtained by the WIEN2k software. The spin-orbit interaction is not taken into account. The crystal structure parameters in the orthorhombic phase is given in Ref. [66]. The orthorhombic structure deformation (a − b)/(a + b) is 0.3%. Due to the electron-phonon interaction, the four-fold symmetry of the bandstructure is slightly violated: The splitting between the d xz -and d yz -bands, ∆E band ≡ E yz − E xz , is 16 meV at X-point, and ∆E band = 2 meV at Γ-point. Figure 6 (b) is the bandstructure in the orthorhombic FeSe. In the orthorhombic phase, the nearest Fe-Fe length is a = 2.6716Å and b = 2.6610Å, so (a − b)/(a + b) is 0.2% [65]. Here, the k-dependent orbital shift to fit the ARPES bandstructure introduced above is not taken into account. In FeSe, ∆E band = 14 meV at X-point, and ∆E band = 3 meV at Γ-point. Thus, the sign reversing orbital splitting observed in Ref. [33] cannot be explained by the band calculation.
The splitting is reduced by the renormalization factor z due to the self-energy. Since z ∼ 1/3 in FeSe and LaFeAsO, the renormalized splitting at X-point is z∆E band ∼ 5meV, which is one order of magnitude smaller than the experimental orbital splitting. Therefore, it is confirmed that the origin of the electronic nematic state in Fe-based superconductors is the electronelectron correlation. In the original SC-VC theory, the spin and charge susceptibilities are calculated self-consistently, by including the MT-VC and AL-VC for the spin and charge susceptibilities [16,49]. The strong orbital fluctuations are induced by the charge-channel AL-VC in Fe-based SCs, Ru-oxides and cuprate SCs [16,37,50]. In the main text, we studied the eight-orbital d-p Hubbard models based on the SC-VC theory, by taking the charge-channel AL-VC into account self-consistently. The obtained χ s (q) is equivalent to the RPA since the spin-channel VCs are dropped. It is easy to verify that the charge-and spinchannel MT-VCs are negligible in the present model. However, the smallness of the spin-channel AL-VC is verified only in the two-orbital Hubbard model in Ref. [50].
Here, we study the FeSe model using the SC-VC method, by taking the MT-VC and AL-VC for both spinand charge-channels in order to confirm the validity of the numerical study in the main text. The charge (spin) susceptibilities arê whereΦ c(s) (q) =χ 0 (q)+X MT,c(s) (q)+X AL,c(s) (q). The spin-channel AL-VC is given as where Λ ′′ m,m ′ ;c,d;g,h (q; p) ≡ Λ c,h;m,g;d,m ′ (q; p) − Λ g,d;m,c;h,m ′ (q; −p − q).
Also, the expressions of the charge-and spin-channel MT-VCs are given in Ref. [49]. The double-counting second-order terms with respect to H U inX MT,s(c) +X AL,s(c) should be subtracted [49] to obtain reliable results.  In the main text, we introduced the first principles models for LaFeAsO and FeSe, and analyzed these models by using the SC-VC method. Here, we also introduce the effective models for BaFe 2 As 2 and NaFeAs, and analyze them using the SC-VC method.
In both BaFe 2 As 2 and NaFeAs, the FSs have relatively large three-dimensional characters. In addition, the unfolding of the bandstructure in BaFe 2 As 2 cannot be exactly performed because of its body-centered tetragonal crystal structure. Here, we introduce an simple effective BaFe 2 As 2 TB model H 0 BaFe2As2 by magnifying the size of the d xy -orbital hole-FS around k = (π, π) in the LaFeAsO unfolded model, in order to reproduce the ARPES bandstructure in Ba122 compounds. Here, we shifted the d xy -orbital band at M point by +0.20 eV. As for NaFeAs, we just use H 0 LaFeAsO as an effective NaFeAs TB model, e.g., H 0 NaFeAs = H 0 LaFeAsO , considering that the FSs in NaFeAs in the k z = 0 plane are similar to the FSs in LaFeAsO. We use H U NaFeAs in place of H U LiFeAs given in Ref. [34]. The bandstructures and the FSs of the effective TB models for BaFe 2 As 2 and NaFeAs are shown in Figs. 8 (a)-(c). Here, we perform the SC-VC analysis for the models H M = H 0 M + rH U M (M=BaFe 2 As 2 , NaFeAs), where r(< 1) is the reduction parameter. We choose the parameter r to satisfy the charge Stoner factor is α C = 0.98; The obtained T -dependences of the spin and charge Stoner enhancement factors, S S ≡ (1 − α S ) −1 and S C ≡ (1 − α C ) −1 respectively, are shown in Fig. 8 (d) and (e). As for BaFe 2 As 2 , both spin and orbital fluctuations strongly develop at T ∼ 50 meV in the case of r = 0.36. This result is consistent with experimental re- lation T mag ≈ T str in BaFe 2 As 2 . As for NaFeAs, only orbital fluctuations strongly develop whereas spin fluctuations remain moderate at T ∼ 50 meV in the case of r = 0.287. This result is consistent with experimental results in NaFeAs [62], in which T mag (= 40K) is more than ten Kelvin smaller than T str (= 53K). Thus, normalstate phase diagrams in BaFe 2 As 2 and NaFeAs are well explained by analyzing their effective Hamiltonians using the SC-VC method. In the main text, we studied the first-principles d-p Hubbard models for LaFeAsO and FeSe by applying the SC-VC theory. In both models, strong spin-fluctuationdriven orbital fluctuations are induced by AL-VC. In FeSe, we found that very small spin susceptibility χ s max is sufficient to realize the orbital order, consistently with experimental results.
Here, we discuss why strong orbital fluctuations are induced by tiny spin fluctuations in FeSe. In Figs. 9 (a) and (b), we show the spin and orbital susceptibilities, χ s max ≡ χ s (Q) and χ c x 2 −y 2 (0) ≡ χ c 2,2;2,2 (q) + χ c 3,3;3,3 (q) − 2χ c 2,2;3,3 (0), in the FeSe model and LaFeAsO model obtained by the SC-VC theory. Here, 32 × 32 k-meshes and 256 Matsubara frequencies are used. In both models, the charge Stoner factor is α C = 0.98 at T = 50 meV, and the obtained orbital susceptibilities show similar Tdependence. We setŪ = 1.76 (r = 0.25) in FeSe, and U = 1.74 (r = 0.41) in LaFeAsO, as we did in the main text. As for the spin susceptibility, in LaFeAsO, strong spin fluctuations develop at T = 50 meV (α S = 0.98), consistently with previous theoretical studies [16,37]. In FeSe, in contrast, χ s max is almost constant till T = 50 meV (α S = 0.87), consistently with experimental reports in FeSe. Now, we discuss why the spin fluctuation strength required to realize α C ≈ 1 is so different from LaFeAsO to FeSe. One reason is the difference in the ratioJ/Ū : Figure 9 (c) shows the T -dependence of the AL-VC on d xzorbital, X AL,c 2 (0) ≡ X AL,c 2,2;2,2 (0), obtained in the LaFeAsO and FeSe models. In both models, α C = 0.98 is satisfied at T = 50 meV. At T = 50 meV, the AL-VC for FeSe is about one-half of that in LaFeAsO. Thus, small AL-VC is enough to induce large orbital fluctuations in FeSe, since the charge Stoner factor is α C ≈ (1 − 5J/Ū )Ū Φ c 2 (0). In Fig. 9 (d), we show that X AL,c: non−zero 2 (0) ≡ X AL,c 2 (0) − X AL,c: zero 2 (0) is very small for both FeSe and LaFeAsO. Here, "zero" represents the zero-Matsubara term (classical contribution) in Eq. (7) in Sec. II. Thus, non-zero Matsubara terms in the AL-VC are negligible in the present calculation (by chance). Note that the U 2 -term in AL-VC gives negative contribution.
Another reason for the relation χ s max (FeSe) ≪ χ s max (LaFeAsO) at α C ≈ 1 is the difference in the orbital dependence of the spin fluctuation strength: The AL-VC for the xz-orbital is approximately given as where we dropped the inter-orbital terms ofχ s andΛ, and leave only the zero-Matsubara term in the Matsubara summation in Eq. (7)    the same model parameters used in Fig. 9. As derived from Fig. 9 (a) and Fig. 10 (a), the ratio χ s 2 (Q)/χ s (Q) is just 0.22 in LaFeAsO, whereas the ratio increases to 0.53 in FeSe, since the relation χ s 4 (Q) ≪ χ s 2 (Q) (χ s 4 (Q) ∼ χ s 2 (Q)) is satisfied in FeSe (LaFeAsO) because of the absence (presence) of h-FS3. This orbital dependence of the spin fluctuations in FeSe is favorable for realizing the FO fluctuations.
To understand the model-dependence of the AL-VC in more detail, we calculate C 2 ≡ q χ s 2 (q) 2 and show the result in Fig. 10 (b): The ratio C LaFeAsO 2 /C FeSe 2 is just 1.35 since the width of the peak of χ s 2 (q) 2 around q = Q is much wider in FeSe. We also examine the square of the three-point vertex for d xz -orbital Λ 2 ≡ Λ 2,2;2,2;2,2 (q, k) at q = 0 and k = Q in Fig. 10 (c). In both models, the relation |Λ 2 | 2 ∝ T a with a ≈ 1 is satisfied for wide temperature range: Such strong T -dependence of the charge-spin coupling Λ 2 is essential for realizing the orbital fluctuations, so it should be taken into account in the numerical calculation. As results, we obtain a crude approximation for the AL-VC,X AL,c 2 ≡ 3U 4 |Λ 2 (0; (0, π))| 2 T C 2 , and show the result in Fig. 10 (d). This crude approximation qualitatively reproduces the exact numerical results for both FeSe and LaFeAsO given in Fig. 9 (c).
In summary, in both LaFeAsO and FeSe, strong orbital fluctuations are induced by AL-VC for the d xz(yz)orbital, X AL,c 2(3) (0). In FeSe, very small spin susceptibility χ s max is sufficient to realize the spin-fluctuation-driven orbital order, because of both the smallness ofJ/Ū and the largeness of C 2 . Strong T -dependence of Λ 2 is essential for realizing the orbital fluctuations due to AL-VC.
Appendix E: Strong T -dependence of the three-point vertex In this paper, we found that the strong orbital fluctuations in Fe-based superconductors originate from the AL-VC for the orbital susceptibility. The moderate increment of the AL-VC at low temperatures shown in Fig. 9 (c) gives rive to the Curie-Weiss behavior of χ c x 2 −y 2 (0). For the increment of the AL-VC, the strong T -dependence of the three-point vertex, shown in Fig. 10 (c), plays the significant role. Its strong T -dependence in Fe-based superconductors had been pointed out in Refs. [20,21,51,52].
Here, we calculate the three-point vertex for LaFeAsO and FeSe models for wide temperature range with high numerical accuracy, using 512×512 k-meshes and ∼ 2048 Matsubara frequencies. Figure 11 shows the square of the three-point vertex for d xz -orbital Λ 2 (0; Q) ≡ Λ 2,2;2,2;2,2 (0; Q) for T ≥ 10 meV. In both LaFeAsO and FeSe models, the coefficient a of |Λ 2 (0; Q)| 2 ∝ T a depends on the temperature range. In both models, a ≈ 1 for T = 20 ∼ 100meV, so the numerical result in Fig.  10 (c) is confirmed by this accurate calculation. When the band renormalization due to z < 1 is considered, the relation a ≈ 1 is realized for T = 20z ∼ 100z [meV].
As shown in Ref. 5 (b), the value of α cr S remains small (∼ 0.9) in the FeSe TB model withĤ U BaFe2As2 (J/Ū = 0.12) or withĤ U LaFeAsO (J /Ū = 0.134). In each case, the obtained T -dependences of S S and S C are qualitatively similar to those shown in Fig. 3 (c). Therefore, the main results of the present study are unchanged even ifJ/Ū in FeSe is slightly larger than 0.1.