Mechanism for nematic superconductivity in FeSe
Jian-Huang She Michael J. Lawler Eun-Ah Kim







^1Department of Physics, Cornell University, Ithaca, NY 14853, USA 

^2Department of physics, Binghamton University, Vestal, NY 13850, USA 

^3Kavli Institute for Theoretical Physics, Kohn Hall
University Of California Santa Barbara CA 93106-4030, USA





Despite its seemingly simple composition and structure, the pairing mechanism of FeSe remains an open problem due to several striking phenomena. Among them are nematic order without  magnetic order, nodeless gap and unusual inelastic neutron spectra with a broad continuum, and gap anisotropy consistent with orbital selection of unknown origin. Here we propose a 
microscopic description of a nematic quantum spin liquid that reproduces key features of neutron spectra.
 We then study how the spin fluctuations of the local moments lead to pairing within a spin-fermion model. 
We find the resulting superconducting order parameter to be nodeless s± d-wave within each domain. Further
we show that orbital dependent Hund's coupling can readily capture observed gap anisotropy. Our prediction 
calls for inelastic neutron scattering in a detwinned sample. 






 Mechanism for nematic superconductivity in FeSe
    Jian-Huang She^1, Michael J. Lawler^2, 1,3, and Eun-Ah Kim^1,3
    December 30, 2023 [file: ]
==================================================================



The pairing mechanism and gap symmetry of 
bulk<cit.> and single layer<cit.> FeSe is an open issue that 
inhibits an overarching understanding of iron-based superconductors. 
 Although a spin-fluctuation mediated pairing scenario is a broadly accepted mechanism in iron-based superconductors<cit.>, much debate continues to focus around two distinct perspectives: weak coupling and strong coupling.
Weak coupling approaches 
 are sensitive to the band structure and generally predict dominantly (π,0), (0,π) spin density wave fluctuations that couple hole pockets to electron pockets in all Fe-pnictides as well as in bulk FeSe <cit.>. Strong coupling approaches take strong electron-electron correlations to generate quasi-localized moments that would interact with itinerant carriers. 
 
 FeSe presents new challenges to both perspectives, including explaining its nematic order <cit.>(see Fig. <ref>(a)), its absence of magnetism, its gapped but active spin fluctuations at (π,π) in addition to (π,0) <cit.> and its nodeless superconducting gap. There have been much efforts to address these issues.  RPA based weak-coupling approaches focused on  implications of assumed nematic order <cit.>. Renormalization group approaches found the effective interactions promoting spin density wave to be also promoting orbital order <cit.>. Approaches focusing on sizable local moments <cit.> led to proposals of quadrupolar order accompanying nematic order <cit.> and the proposal of 
a quasi-one dimensional quantum paramagnet state<cit.>  of AKLT (Affleck-Kennedy-Lieb-Tasaki)<cit.> type. Nevertheless, strikingly unique inelastic neutron spectra (INS) of FeSe evade the approaches so far one way or another. 
 
 





The absence of the stripe order in FeSe has been attributed to the notion of frustration<cit.>.
Indeed FeSe is close to a classic situation for frustrated magnets in the much studied J_1-J_2 model<cit.>(see Fig. <ref>(b)). 
Interestingly, in systems that form stripe upon cooling, viewing the nematic state as thermally melted version of stripe was a very productive point of view <cit.>. Here we note that frustration from the competition between J_1 and J_2 has been long known to drive quantum melted versions of Neel and stripe orders giving rise to C_4 symmetric and C_2 symmetric (nematic) quantum spin liquids (QSL) respectively<cit.>. 
Moreover DMRG studies on J_1-J_2 model noted an intermediate paramagnetic phase between stripe order and Neel order state <cit.>. A recent DMRG study of J_1-J_2-K_1-K_2 spin model found a nematic quantum paramagnetic state between the Neel and stripe ordered states<cit.>. 
In this letter we propose a microscopic description of the frustration driven nematic quantum spin liquid (QSL)  state that amounts to quantum melted stripe and
captures the observed INS. We then investigate the implication of dramatically anisotropic spin-fluctuation spectra of the proposed state on the nature of superconductivity.  







In FeSe, there is evidence that local moments <cit.> coexist with itinerant carriers of all three t_2g orbitals<cit.>.
In order to capture the dual character <cit.> 
we turn to a spin-fermion model <cit.>: H= H_c+ H_S+ H_ int, where H_c and H_S describe the itinerant carriers and local moments respectively that are coupled through  H_ int. 
 For the spin model 

    H_S= ∑_ij J_ij S_i· S_j,

with exchange interactions J_ij on a square lattice (Fig.1b), the two dominant interactions 
are the nearest-neighbor J_1 and the next-nearest-neighbor J_2 exchange interactions as in other Fe-based superconductors <cit.>. But due to the near itinerancy of the core electrons, longer range terms  
are also expected <cit.>.


The J_1-J_2 model has been extensively studied both classically and quantum mechanically
(see Refs.[Chandra90, Misguich04, Jiang09, Jiang12]). Within classical models the role of frustration is clear from the fact that the model can be recast as
H_S=J_2∑ ( S_1+ S_2+ S_3+ S_4)^2 up to a constant at J_2=J_1/2 point, where S_1-4 are the four spins on each plaquette ⟨ 1234⟩ and the summation is over all plaquettes.
Classical ground state with vanishing total spin on each plaquette property leads to a zero mode at each wave vector on the Brillouin zone boundary<cit.> and so the model is highly frustrated.
With quantum effects of small spin S the frustration effects are not limited to the fine tuned point of J_2=J_1/2. Unfortunately, a controlled theoretical study for quantum spins for such frustrated spin systems is challenging. Hence we will restrict ourselves to mean field theories and choose an ansatz that (1) agrees with the observed inelastic neutron spectrum <cit.>, and (2) the ordering tendencies obey the classical condition of S_1+ S_2+ S_3+ S_4=0 on a plaquette.

 A prominent feature of the INS data <cit.> is its broad and gapped continuum of spectral weight (Fig.<ref>a) without any one-magnon branch. Intriguingly such a continuum is expected in a QSL with deconfined spinons in two-dimension in an insulating magnetic system <cit.>. Indeed it is a common feature of slave-particle mean field theories. So we will choose Schwinger boson mean field theory (SBMFT) <cit.> as our mean field theory approach. Additional features of Fig. <ref> (a-c) we aim to capture include:


  * The simultaneous presence of both (π, π) spin fluctuations and (π, 0), (0, π) spin fluctuations. 

  * The quasi-one-dimensional dispersion ω∼sin k_y  <cit.> found in the shape of the upper and lower bounds.

  * The observed cross-shaped spectrum around (π, π).









To find these features in a SBMFT, we turn to the known <cit.> SBMFT phase diagram of the J_1-J_2 model (Fig. <ref>). 
Note that the Neel and stripe long range order for small J_2/J_1 and large J_2/J_1 are expected<cit.> to melt into C_4 symmetric and C_2 symmetric QSL's respectively (see Fig. <ref>).
Hence the shaded region near the phase boundary between C_4 symmetric QSL, C_2 symmetric QSL and the stripe ordered phase will capture all of the above features. Specifically, states in this region will support a dynamic spin structure factor with 1d-like dispersion and cross-shaped spectrum 
assuming twin domains of the stripe state are averaged over in the INS data. To account for the itinerancy of the electrons, we extend  
an ansatz within the shaded region of Fig. <ref>  
with additional neighbor couplings.

To construct the ansatz, we now turn briefly to the specifics of SBMFT. In Schwinger boson representation, each spin S_ r  is represented by two bosonic operators b_ rσ, σ=↑,↓ with the constraint ∑_σ b^†_ rσb_ rσ=2S. The spin operator is then S_ r=1/2∑_σσ'b^†_ rσσ_σσ'b_ rσ', with σ the Pauli matrices. We can then expand H_ r, r'≡ J_ r, r' S_ r· S_ r' in terms of the spin singlet operator A^†_ r, r'=b^†_ r↑b^†_ r'↓-b^†_ r↓b^†_ r'↑ to obtain H_ r, r' = -J_ r, r'1/2A^†_ r, r'A_ r, r'+S^2. Finally, 
we mean-field decompose H_ r, r' and introduce mean fields ⟨ A_ r, r'⟩ using A^†_ r, r'A_ r, r' = ⟨ A^†_ r, r'⟩ A_ r, r' + A^†_ r, r'⟨ A_ r, r'⟩ - ⟨ A^†_ r, r'⟩⟨ A_ r, r'⟩. We will further assume the bosons do not
 condense for we are interested in the spin liquid phase.  


Defining A_μ̂≡⟨ A_ r, r+μ̂⟩, we keep  A_x̂≠ 0 and the diagonals A_x̂±ŷ≠ 0 and A_x̂± 2ŷ≠ 0
for states in the shaded region of Fig. <ref>. 
The fourth neighbor term can be understood as a result of the competition between Neel and stripe states: it is a bond that is favored by both the (π,π) Néel state and the (π,0)/(0,π) stripe state.
The result is a state with the same projective symmetry group as the Read and Sachdev state used in the phase diagram of Fig. <ref>. It is a “zero flux state" in that the smallest loop has zero “flux" obeying the so-called flux expulsion principle <cit.> and hence energetically competitive. 
Most importantly it is a state in which translational symmetry is restored by quantum melting stripe into C_2 symmetric nematic QSL state. 


We can then calculate 
the dynamic spin structure factor S_ qω≡ Im⟨ S^z( q, ω) S^z(- q, ω)⟩ associated with our ansatz.
 At T=0, it is of the form <cit.>

    S_ q, ω∼∑_ k{cosh[ 2(θ_ k+θ_ k+q̃)]-1}δ( ω_ k+ω_ k+q̃-|ω|),

where θ_ k is the angle in the Bogoliubov transformation of SBMFT (see SM1 for explicit expression), and q̃= q-(π, 0) arises because of a standard unitary transformation we carried out on the B sublattice for simplicity. The results summing over two domains are plotted in Fig. <ref>(d-f). They capture the basic features of the neutron spectra: (1) The spectrum is gapped (Fig. <ref>d), as a result of the absence of long range magnetic ordering. 
 (2) Both (π, π) and (π, 0)/(0, π) spin fluctuations are present (Fig. <ref>d, e).   
  (3) The spectrum displays the novel feature of continuum with the bounds exhibiting quasi-one-dimensional dispersion (Fig. <ref>d).
  
 A sharp prediction of our model is the dramatic suppression of 
 spectral weight around (0, q_y) in a detwinned sample ((q_x, 0) for the other domain). This means at low energies there are weights at say (π, π) and (π, 0), but not at (0, π). 
 By contrast, in an orbital order driven picture for nematic ordering, there is only a weak anisotropy in the spin-structure factor with the spectral weight at (π, π), (0, π) and (π, 0) of roughly the same magnitude even in a single nematic domain<cit.>.
 Such a distinction has profound implications for pairing.  When the degree of anisotropy in the momentum distribution of the spin spectra is mild, pairing interactions with different q-wavevectors compete, leading to nodes <cit.>. On the other hand, the strong anisotropy in the spectral weight distribution in our SBMFT ansatz quenches such competition removing any need for a superconducting gap node. 


 

 
We now turn to the itinerant degrees of freedom to study nematicity and superconductivity. Their kinetic energy is given by a tight-binding model:
 
    H_c=∑_ k,αβ,νϵ_αβ^μν( k)c^†_αμ( k)c_βν( k),

where c^†_αμ( k) creates an itinerant electron with momentum k, spin μ and orbital index α. 
The Fermi surface of FeSe consists of two electron pockets around the M points and one hole pocket around the Γ point <cit.>. Following <cit.>, we take a simple symmetry based approach of expanding the dispersion ϵ_αβ^μν( k) around the Fermi surface. It is known experimentally that the spectral weight of the low energy states are predominantly from d_yz and d_zx around the Γ point, from d_yz and d_xy around (π, 0), from d_zx and d_xy around (0, π). We consider the corresponding intra- and inter-orbital hopping terms. Furthermore we include on-site nematicity and spin-orbit coupling to produce the band splitting that gives rise to a single hole pocket around Γ. The resulting simplified Fermi surface is shown in Fig.<ref>a, see SM2 for explicit parameters. [Our rather simple band structure allows largely analytic calculation at the expense of missing out quantiative details such as large mismatch in the pocket sizes as found in quantum oscillation (see SM4).]


The itinerant electrons couple to the local moments via the ferromagnetic Hund's coupling <cit.>:

    H_ int=-∑_i,α,μνJ_α S_i· c^†_iαμσ_μνc_iαν,

where σ represents the vector of Pauli matrices, and J_α>0 denote the Hund's couplings. 
Since the Hund's couplings depend on the overlap of the itinerant electron wave function with the local moment, they are generally different for different orbitals. 
[Coupling to itinerant electrons generates a self-energy for the local moment propagator, giving rise to Landau damping (see SM3). However the strength of the Hund's couplings can be estimated to be much weaker than the spin exchange interaction (see SM3). Hence coupling to itinerant electrons will not significantly modify the local moment spin susceptibility.]




 Note that the proposed nematic QSL state induces nematicity in the charge sector. 
  For instance non-zero 
  ⟨ A_ r, r±x̂⟩ in the nematic QSL state generates 
  an interaction among conduction electrons along the x-direction,  
  which drives bond-centered nematic order with 
  φ_c≡⟨ c^†_ r+x̂,αc_ r,α -c^†_ r+ŷ,αc_ r,α⟩≠ 0 below the temperature at which the nematic QSL develops. The observed nematic transition at T_s∼ 90K <cit.> is consistent with this picture. Furthermore, φ_c linearly couples to 
   φ_o≡n_zx-n_yz/n_zx+n_yz, where n_zx, yz denote occupation of zx and yz orbitals, and φ_s ≡ M_x^2-M_y^2, where M represents the magnetic moment. These different measures of nematicity are consistent with orbital imbalance 
   observed in ARPES <cit.> (φ_o≠ 0) and the  observed NMR resonance line splitting <cit.> (φ_s≠ 0). 




       






  
Furthermore, the nematic spin fluctuations in the proposed QSL state mediate pairing 
and the resulting gap structure can be determined via standard a mean field procedure (see SM4). An immediate observation is that  non-universal aspects of the gap structure such as relative gap strength of each pocket and the T_c are sensitive to strength of the Hund's couplings J's (see Fig. <ref>b,c). Nevertheless the gap functions resulting from our model share the following generic features: (1) The gap is generically nodeless as a result of severe anisotropy of the spin fluctuations in the nematic QSL state. In particular, the near absence of spin fluctuations around say (0, π) for one nematic domain renders the determination of gap sign in different pockets unfrustrated. In contrast, in the itinerant model where (π, π), (π, 0) and (0, π) spin fluctuations are close in magnitude they compete for deciding the sign structure of the gap causing nodal gap structures. (2) The gap is deeply anisotropic due to the 
variation of orbital content around each Fermi pocket. The resulting nodeless but very anisotropic gap structure explains the seemingly contradictory experimental results of STM <cit.>, penetration depth <cit.>, thermal conductivity measurements <cit.>, observing low energy excitations <cit.> despite the evidence of a full gap <cit.>. (3) 
The gap changes sign from pocket to pocket. This is consistent (see SM4) with the observation of sharp spin resonance in the superconducting state <cit.>. 
More specifically, our gap function is a combination of d-wave as induced by (π, π) spin fluctuations and s_± as induced by (π, 0) spin fluctuations.
 We consider a single nematic domain, where the pairing interaction concentrates around (π, q_y). Two examples of the gap function (in arbitrary units) are shown in Fig.<ref>b,c, where J_xy=J_yz=J_zx and J_xy=0.4 J_yz=0.4 J_zx respectively. 



Now we turn to the question of orbital dependence of the Hund's coupling. Fig. <ref>b,c shows that the orbital dependent Hund's coupling can alter the relative magnitude and anisotropy of gap functions at different Fermi pockets (while the gap is predominantly d-wave in Fig. <ref>b, d- and s-wave are at par in Fig. <ref>c). 
Since the Hund's coupling requires overlap of the wave-function between the conduction electrons and local moments, significantly lower spectral weight of d_xy orbitals <cit.> implies J_xy≪ J_zx, J_yz.
 

Indeed, the gap function with such orbital dependent Hund's coupling shows remarkable resemblance to the gap structure observed by recent STM measurements<cit.> (see Fig. <ref>c,d).  In Ref.  Sprau16 the observed pocket specific gap anisotropy was interpreted as resulting from orbital-selective pairing of unknown microscopic origin.
In our model, such orbital selective pairing arise from orbital dependence in the Hund's coupling J_xy<J_yz=J_zx, reflecting much smaller weight of d_xy orbitals in the conduction electrons <cit.>. This orbital dependent Hund's coupling amplify the role of (π,0) spin fluctuation in pairing despite larger spectral weight at (π,π), which is consistent with the observation of sharp spin resonance at (π,0) <cit.> (see SM4 for further discussion).  


 
In conclusion, we propose a nematic QSL state description of FeSe that explains 
the basic phenomenology of FeSe: (1) Spin dynamics observed in Ref.Zhao16 assuming it is averaged over domains, (2) nematic transition without mangetic ordering, (3) highly anisotropic fully gapped superconducting gap. 
The central assumption that neutron scattering is averaging over domains could be tested in a detwinned neutron experiment. Orbital dependent Hund's coupling mechanisms for orbital selective pairing in bulk FeSe further offers new insight regarding higher T_c observed in 
mono-layer FeSe and K-doped FeSe. As we show in SM4, larger J_xy that enables conduction electrons to utilize (π,π) spin fluctuation with larger intensity and higher characteristic frequency leads to higher transition temperature (as high as  47K). Combined with the observation that spectral weight of the d_xy orbitals in the conduction electrons is much higher in the higher T_c settings of  mono-layer FeSe and K-doped FeSe <cit.>, it is conceivable these systems make better use of already more prominent (π,π) fluctuation to achieve higher T_c. We note here that the nematic QSL state we propose is distinct from the proposal of Ref. FWang15 in that it contains no one-magnon branch of excitations (see SM5) although both proposals start from strong coupling perspective and spin ground states lacking any form of magnetic order.
Finally, although we used SBMFT as a calculational crutch to capture the spin wave continuum, the ultimate fate of spinons in this spin system coupled to itinerant electrons needs further study. Interestingly, such a state with spinons coexisting with conduction electrons would resemble the FL* state first proposed in Refs. Senthil03, Senthil04 that has recently been revisited using DMRG <cit.>.




Acknowledgements We thank Andrey Chubukov, J.C. Davis, Rafael Fernandez, Yong Baek Kim, Steve Kivelson, Igor Mazin, Andriy Nevidomskyy, Subir Sachdev, Doug Scalapino, Qimiao Si, Fa Wang   for discussions. E-AK  and J-HS were  supported by the U.S. Department of Energy, Office of Basic Energy Sciences,  Division of Materials Science and Engineering  under  Award  DE-SC0010313.   E-AK  also  acknowledges Simons Fellow in Theoretical Physics Award#392182. E-AK and MJL acknowledge hospitality of the KITP supported by Grant No. NSF PHY11-25915.


apsrev











§ SM1: SCHWINGER BOSON MEAN FIELD THEORY


We show here that our ansatz state is a self-consistent solution of the J_1-J_2-J_3-J_4 spin model (see Fig.1a of main text for definition of J's). On a bipartite lattice, it is convenient to perform a unitary transformation by defining a_im≡ b_im on A sublattice, and a_j↑≡ b_j↓, a_j↓≡ -b_j↑ on B sublattice. The valence bond operator is then brought to the simpler form A^†_ij= ∑_σ a^†_iσa^†_jσ. Modular a constant, the spin Hamiltonian 
H_S=∑_ijJ_ij S_i· S_j
  can be written in terms of the valence bond operators as
   
    H_S=-1/2∑_ijJ_ijA^†_ijA_ij.

    We then apply mean field theory to the bosonic Hamiltonian <cit.>. Defining Q_ij=J_ij⟨ A_ij⟩≡ Q_δ, the quadratic part of the mean field Hamiltonian reads:  
    H_S^( MF)=λ∑_iσ a^†_iσa_iσ+1/2∑_iδσQ_δ( a^†_iσa^†_i+δ, σ+ a_iσa_i+δ, σ).

  For a given mean field ansatz, the mean field Hamiltonian can be diagonalized by the Bogoliubov transformation

    α_ kσ=coshθ_ ka_ kσ-sinhθ_ ka^†_- kσ,

   with tanh( 2θ_ k)=-Qγ_ k/λ. Here Qγ_ k denotes the Fourier transform of Q_δ: Qγ_ k=∑_δQ_δe^-i k·δ. 
 The resulting Hamiltonian reads 
 
    H_S^( MF)=∑_ kσω_ k(α^†_ kσα_ kσ+1/2),

with the dispersion ω_ k=√(λ^2-(Qγ_ k)^2).  Integrating out the bosonic fields, one obtains the free energy 
  
    F=∑_δ|Q_δ|^2/2J_δ -1/2(2S+1)λ+1/β∫d^2 k/(2π)^2ln[2sinh(1/2βω_ k)],

from which follow the self-consistency equations. In SBMFT, long-range order occurs through Bose-Einstein condensation (BEC) of the Schwinger bosons, and condensation gives rise to gapless spectrum due to the resulting Goldstone mode. A QSL state corresponds to a solution of the self-consistency equations with gapped spectrum, where there is no condensation of Schwinger bosons.

We start with decoupled 1d chains where only Q_x≠ 0, and 
 
    Qγ_ k=2Q_xcos k_x.

 Its free energy is
  
    F^(1)=Q_x^2/J_1-1/2(2S+1)λ+1/β∫d^2 k/(2π)^2ln[2sinh(1/2βω_ k)],

   and the self-consistency equations are 
 
    S+1/2   =   ∫d^2 k/(2π)^2λ/2ω_ k, 
     
    Q_x/J_1   =   ∫d^2 k/(2π)^2Qγ_ k/2ω_ kcos k_x.
  
This ansatz state is basically determined by the dimensionless parameters: Q_x/λ. We have plotted the spin structure factor taking Q_x/λ=0.498 (see Fig. <ref>a,b,c). Here S=0.677.
  
 We consider then our ansatz state, namely the quantum melted stripe state, where Q_x≠ 0, Q_x+y≠ 0, Q_x+2y≠ 0, and
     
    Qγ_ k=2Q_xcos k_x+4Q_x+ycos k_x cos k_y  + 4Q_x+2ycos k_x cos(2k_y).

Its free energy is
  
    F^(2)=Q_x^2/J_1+ 2Q_x+y^2/J_2+2Q_x+2y^2/J_4-1/2(2S+1)λ+1/β∫d^2 k/(2π)^2ln[2sinh(1/2βω_ k)],

   and the self-consistency equations are 
 
    S+1/2   =   ∫d^2 k/(2π)^2λ/2ω_ k, 
     
    Q_x/J_1   =   ∫d^2 k/(2π)^2Qγ_ k/2ω_ kcos k_x, 
    
    Q_x+y/J_2   =   ∫d^2 k/(2π)^2Qγ_ k/2ω_ kcos k_xcos k_y ,
    
    Q_x+2y/J_4   =   ∫d^2 k/(2π)^2Qγ_ k/2ω_ kcos k_xcos(2k_y).
  
This ansatz state is basically determined by the three dimensionless parameters: Q_x/λ, Q_x+y/Q_x and Q_x+2y/Q_x. We have plotted the spin structure factor taking Q_x/λ=0.398, Q_x+y/Q_x=0.025 and Q_x+2y/Q_x=0.1 (see Fig. <ref>d,e,f). Here S=0.153, J_2/J_1=0.904, J_4/J_1=0.975.


In addition, we find that the quantum melted stripe state has lower energy than the decoupled 1d chain state. We set J_1=1. For S=0.25,  J_2=0.88, J_4=0.944, we obtain F^(1)=-0.175,  F^(2)=-0.179. For S=0.26, J_2=0.864, J_4=0.925, we obtain F^(1)=-0.188, F^(2)=-0.192. For S=0.3, J_2=0.869, J_4=0.935, we obtain F^(1)=-0.224, F^(2)=-0.227. 

 






§ SM2: ITINERANT PART: THREE ORBITAL MODEL



For the itinerant part of the system, what matters for pairing are the low energy electronic states around the Fermi pockets. We take a phenomenological approach to expand the dispersion around the Fermi pockets. Since it is known experimentally that the spectral weight of the low energy states arises mainly from d_yz, d_zx, d_xy orbitals, we consider a band structure involving these three orbitals. Such orbital-projected band models have been studied in <cit.>. Consider first the Fermi pocket near the Γ point, a single hole pocket has been detected, with the d_yz and d_zx orbitals dominating the spectral weight. We introduce a spinor ψ^T_γ,  k=(c_yz,  k↑, -c_zx,  k↑, c_yz,  k↓, -c_zx,  k↓), and the kinetic energy part of the Hamiltonian is of the form 

    H_0,Γ = ∑_ kψ^†_γ,  k h_Γ( k)ψ_Γ,  k.

The Hamiltonian includes the on-site energy, intra-orbital hopping, inter-orbital hopping. To get the right orbital splitting, we include also the difference in the on-site energy for the two orbitals reflecting nematicity, and the spin-orbit coupling <cit.>. The result reads

    h_Γ( k)=(ε_Γ+k^2/2m_Γ) τ^0⊗σ^0 + [δε_Γ +b k^2cos (2θ_k)] τ^3⊗σ^0 + ck^2sin (2θ_k) τ^1⊗σ^0 + λτ^2⊗σ^3,

where τ and σ are Pauli matrices in orbital and spin space respectively, and k=(k_x, k_y)=k(cosθ_k, sinθ_k).


For the electron pocket near (π, 0), the d_yz and d_xy orbitals dominate the spectral weight.
We introduce a spinor ψ^T_X,  k=(c_yz,  k↑, c_xy,  k↑, c_yz,  k↓, c_xy,  k↓), and the kinetic energy part of the Hamiltonian is of the form 

    H_0,X = ∑_ kψ^†_X,  k h_X( k)ψ_X,  k,

with 
 
    h_X( k)=[ε_1 +k^2/2m_1-a_1 k^2cos(2θ_k)] τ^0+τ^3/2⊗σ^0 +[ε_3 +k^2/2m_3-a_3 k^2cos(2θ_k)] τ^0-τ^3/2⊗σ^0+2vksinθ_k τ^2⊗σ^0.

Here k is measured from (π, 0).
 

For the electron pocket near (0, π), the d_zx and d_xy orbitals dominate the spectral weight. We introduce a spinor ψ^T_Y,  k=(c_zx,  k↑, c_xy,  k↑, c_zx,  k↓, c_xy,  k↓), and the kinetic energy part of the Hamiltonian is of the form 

    H_0,Y = ∑_ kψ^†_Y,  k h_Y( k)ψ_Y,  k,

with 
 
    h_Y( k)=[ε_1 +k^2/2m_1+a_1 k^2cos(2θ_k)] τ^0+τ^3/2⊗σ^0 +[ε_3 +k^2/2m_3+a_3 k^2cos(2θ_k)] τ^0-τ^3/2⊗σ^0+2vkcosθ_k τ^2⊗σ^0.

Here k is measured from (0, π).

With a proper choice of the parameters, we can obtain a single hole pocket around Γ, a single electron pocket around (π, 0), and a single electron pocket around (0, π) as shown in Fig.4a of main text. The corresponding parameters are: ε_Γ =14, δε_Γ = 11, 1/2m_Γ=-350, b=-70, c=120, λ = 9, ε_1=-20, ε_3=-60 1/2m_1=75, 1/2m_3=160, a_1=100, a_3=-120, v=-60. Note that the band structure employed in our paper is simplified in order to use a closed form for hamiltonian that allows us to carry out the study of superconductivity semi-analytically. Although our band structure misses  quantitative details such as severe mismatch in the pocket sizes between two electron pockets (see SM 4 for further discussion), such details will not impact qualitative conclusions of the paper. 



§ SM3: HOW ITINERANT ELECTRONS AFFECT LOCAL MOMENTS: LANDAU DAMPING


Coupling to itinerant electrons generates a self-energy for the local moment propagator, giving rise to Landau damping. Since the Fermi pockets are small in size, and located near Γ- and M-points, the induced self-energy will be predominantly near q=(π, π), (π, 0) and (0, π) (see Fig.<ref>). We expect the neutron spectrum near these points to be smeared. For q∼ (π, 0), the self-energy is predominantly from d_yz orbitals, 
 
    D^(π, 0)( q, Ω) ∼∑_ k,ωJ_yz^2 G^(yz)( k, ω) G^(yz)( k+ q, ω+ Ω),

with fermion Green's function G( k, ω). Here k is at the Γ pocket, and k+ q at the M_x pocket. For q∼ (π, π), the self-energy is predominantly from d_xy orbitals, 
 
    D^(π, π)( q, Ω) ∼∑_ k,ωJ_xy^2 G^(xy)( k, ω) G^(xy)( k+ q, ω+ Ω),

where k and k+ q are at the two M-pockets. With the suppression of J_xy, we expect the Landau damping effect to be weaker near q=(π, π).





We can estimate the strength of the Hund's coupling from the resulting superconducting transition temperature: T_c∼ E_F e^-1/λ, with the dimensionless coupling λ= N_0 V∼ J_H^2/(E_FJ_ ex). The energy scales involved are: (1) Fermi energy E_F∼ 10 meV <cit.>, (2) the exchange interaction J_ ex∼ 100 meV <cit.>, and (3) T_c∼ 1 meV. From these, we obtain Hund's coupling J_H∼ 20 meV, which is much smaller than the exchange interaction J_ ex. Hence we expect coupling to itinerant electrons will not significantly modify the local moment spin susceptibility.



§ SM4: HOW LOCAL MOMENTS AFFECT ITINERANT ELECTRONS: NEMATICITY AND PAIRING


The dynamic spin fluctuations in the QSL affect the itinerant electrons. Since the spins have a gapped spectrum, we can integrate them out to obtain an effective interaction for the itinerant electrons. The induced action reads

    S_ int =-1/2∫_0^β dτ∑_α, α' J_αJ_α'χ_ij(τ) s_iα(τ)· s_jα'(0),

with the itinerant electron spin density s_iα≡∑_μνc^†_iαμσ_μνc_iαν, and the local moment spin correlation function χ_ij(τ)≡⟨ T_τS_i^a(τ)S_j^a(0) ⟩. The induced interaction is highly anisotropic, and the dominant interaction term is the nearest-neighbor interaction (say along the x-direction): J_H^2χ c^†_ rασ^a_αβc_ rβc^†_ r+x̂,α'σ^a_α'β'c_ r+x̂,β'. This interaction results in a phase transition to a nematic state with order parameter ⟨ c^†_ r+x̂,αc_ r,α⟩≠ 0, or more generally, φ_c≡⟨ c^†_ r+x̂,αc_ r,α -c^†_ r+ŷ,αc_ r,α⟩≠ 0.

Furthermore, the induced interaction leads to pairing among the itinerant electrons.
Since the spin fluctuations are antiferromagnetic, one expects pairing in the spin singlet channel. We then mean field decompose the induced interaction into spin singlet pairing channel with the corresponding pair operator h^†_αα'( k)=1/√(2)(c^†_ kα↑c^†_- kα'↓ -c^†_ kα↓c^†_- kα'↑). Due to the special form of spin susceptibility and band structure in FeSe, the pairing problem is largely simplified. The spin fluctuations enter the pairing problem through the spin susceptibility χ( q)≡χ( q, Ω_n=0), which can be obtained from the dynamic spin structure factor via χ( q)=-∫ dω S_ q, ω/ω. The special form of χ( q) in FeSe (see Fig. <ref>) results in only inter-band pairing correlation among the three Fermi pockets. Furthermore, since Hund's coupling is diagonal in orbital space, there are only pairing correlations between the same orbitals: in orbital basis, the pairing interaction is of the form H_ pair∼ J_α^2 χ( k- k')h^†_αα( k)h_αα( k'). 



Pairing occurs near the Fermi surface, which is naturally expressed in the band basis. We then transform from the orbital basis to the band basis: c^†_ kαμ=∑_a η^*_α a μ( k)d^†_ kaμ with band index a. Note that since the spin-orbit coupling here is in the σ^3 channel, different spins do not mix. The pair operator in the band basis is h^†_a( k)=1/√(2)(d^†_ ka↑d^†_- ka↓ -d^†_ ka↓d^†_- ka↑). Omitting the frequency dependence, the pairing Hamiltonian is of the form 

    H_ pair =∑_ k k'abΓ_ab( k,  k')h^†_a( k)h_b( k'),

 with the projected pairing interaction Γ_ab( k,  k')=1/2∑_αJ_α^2χ( k- k') M^*_α a( k) M_α b( k'). The orbital content is encoded in the form factor M_α a( k)= η_α a ↑( k) η_α a ↓(- k).
  The gap function is then defined as Δ_a( k) =∑_ k'bΓ_ab( k,  k')⟨ h_b( k')⟩. The gap symmetry function g_i( k)∝Δ_a( k) on the Fermi surface is determined by the eigen equation

    -∑_j∮_ FS_jd k'_∥/2π v_F( k')Γ_ij( k,  k')g_j( k') =λ g_i( k),

where k_∥ denotes momentum along the Fermi surface FS_j, and v_F( k)=|∇ E_a( k)| represents the Fermi velocity. We can then solve the above eigen equation to find the leading eigenvalue and the corresponding eigenvector, which determines the resulting gap structure within a single nematic domain. The inputs are (1) itinerant electron band structure encoded in ϵ_αβ^μν( k) (2) local moment spin susceptibility χ( q), and (3) Hund's couplings J_α. 


We show here more concretely how (π, π) spin fluctuation mediated pairing enhances the superconducting T_c.
We first estimate the value of λ in bulk FeSe from the observed T_c. Since Fermi energy is small compared to spin fluctuation scale (so called antiadiabatic limit), Fermi energy acts as cutoff in the T_c equation: T_c∼ E_F e^-1/λ <cit.>. With T_c∼ 8 K, E_F∼ 10 meV <cit.>, we obtain λ∼ 0.37. In bulk FeSe, pairing occurs predominantly among d_yz orbitals as mediated by (π, 0) spin fluctuations, while (π, π) spin fluctuation mediated pairing among d_xy orbitals is largely suppressed. This corresponds to taking (J_xy, J_zx, J_yz) ∼ (0, 1, 1). (Note that due to the near absence of (0, π) spin fluctuations, pairing among d_zx orbitals is suppressed for any coupling. So we just set J_zx =1.) When pairing is predominantly among d_xy orbitals, we have (J_xy, J_zx, J_yz) ∼ (1, 1, 0). We have obtained the resulting eigenvalue  λ'=2.98λ, which gives T_c∼ 47 K. Hence (π, π) spin fluctuation mediated pairing is indeed able to account for the much higher T_c in heavily doped FeSe (T_c∼ 48 K <cit.>), and a large part of the T_c increase in monolayer FeSe (T_c∼ 50-64 K <cit.>).  



 §.§ Neutron resonance


A characteristic feature of our gap function is that it has different signs at different Fermi pockets. Such a sign-changing gap function can give rise to resonances in the neutron spectrum. However the intensity of the resonances depends on the details of the band structure and the superconducting gap function. The resonance comes from the itinerant electron spin susceptibility, which contains a term of the form<cit.>

    χ”( q, ω)∼∑_ k( 1-Δ_ kΔ_ k+ q/E_ kE_ k+ q)
    [ 1-f(E_ k)-f(E_ k+ q) ]δ(ω-E_ k-E_ k+ q ).

The gap function changes sign between the Γ pocket and (π, 0) pocket, and between the (π, 0) pocket and (0, π) pocket. Since the sizes of the Fermi pockets are small, one expects the  resonances to be localized in momentum space around q=(π, 0) and q=(π, π). 


Due to lack of full knowledge of the gap function in the whole Brillouin zone, we approximate the neutron intensity at q=(π, 0) and q=(π, π) by summing over the corresponding Fermi surfaces. Furthermore we penalize the resulting term by an exponential factor depending on the difference between the momentum transfer and q, i.e. we replace ∑_ k by ∮_ k∮_ k'e^-| k- k'- q|^2/q_0^2, where the integrals are restricted to the Fermi surfaces, and q_0 is a parameter that basically measures the range of blurring in momentum space. Using the simplified band structure employed in the paper and specified in SM3 (Fig.4(a) of main text), and the resulting gap function with J_xy/J_yz=0.4 (see Fig.4(c) of main text), we obtain the neutron intensity as shown in Fig.<ref>(a) (with q_0=0.2). One can see that the resonance at q=(π, π) is suppressed due to worse nesting of the corresponding Fermi surface and gap functions. 

Actually one expects further suppression of (π, π) resonance with more realistic band structures.  In particular, it has been found in quantum oscillation measurements <cit.> that the Fermi pocket at (0, π) (with k_F≃ 0.13Å^-1) is much larger than the Fermi pocket at (π, 0) (with k_F≃ 0.043Å^-1). Indeed a simple check by enlarging the Fermi pocket at (0, π) by a factor of two while keeping the rest of the band structure and gap functions fixed suppreses the resonance at q=(π, π) to be almost vanishing (see Fig.<ref>(b)).    










§ SM5: AKLT CHAINS OR FRACTIONALIZED SPIN LIQUID IN FESE?

The theory of FeSe presented in the main manuscript consisted of coupling a nematic quantum spin liquid of local moments described by Schwinger boson mean field theory to itinerant electrons. In the absence of this coupling, the local moments would behave like they do in insulators where it is known <cit.> that this Schwinger boson mean field theory is unstable by confinement and forms chains of AKLT states. This confined state is presumably related to the AKLT chain states studied in Ref. FWang15 as a potential theory of the magnetism in FeSe. 

Here we attempt to find an experimental signature that could distinguish the deconfined state studied in the main manuscript and an AKLT chain state in a future experiment. Our approach will be to make use of known results from one dimensional physics. If we decouple chains in the AKLT chain state, we can use numerical results on the well studied one dimensional model to determine signatures in neutron scattering for such a state. Then after determining experimentally relevant features of this state we will step back and assess the implication of these results for the broader question of the distinction between the AKLT chain state and nematic quantum spin liquid states.








In one dimension, the S=1 spin chain with bilinear and biquadratic interactions has a variety of phases including the AKLT state and a variety of phase transitions. The Hamiltonian is

    H_S=J∑_i[cosθ( S_i· S_i+1) +sinθ( S_i· S_i+1)^2  ],

with -π≤θ≤π. As shown in Fig.<ref>, its phase diagram contains (see <cit.> and references therein): (1) a ferromagnetic phase for -π<θ<-3π/4 and π/2<θ≤π, (2) a dimerized phase for -3π/4<θ<-π/4, (3) a gapped and topologically ordered Haldane phase for -π/4<θ<π/4 (the AKLT state corresponds to θ=arctan1/3≃ 0.1024π) and (4) a gapless phase with antiferroquadrupolar (AFQ) correlations for π/4<θ<π/2 <cit.>. The dimerized phase and the Haldane phase are separated by a critical point, the Takhtajan-Babujian (TB) point <cit.>. At the TB point, the system possesses gapless spinon excitations. The spinon continuum is manifest in the dynamical spin structure factor (Fig.<ref>) as obtained using algebraic Bethe ansatz-based method <cit.>. As one moves away from the critical point into the Haldane phase, the spinons get confined, and the elementary excitations are magnons. However the magnons are strongly interacting and not always well defined: in addition to the one magnon branch, the two-magnon processes have important contributions to the dynamic spin structure factor (see Fig.<ref>) <cit.>. Returning back to the critical point, the two-magnon excitations merge with the one magnon excitations and only a continuum of spinon excitations remain<cit.>.

Presumably, the qualitative features of the one dimensional model would carry over to a coupled two dimensional chain model. This implies that an AKLT chain state would similarly consist of both a continuum of excitations and a one magnon branch similar to those found in Fig. <ref>. But the nematic quantum spin liquid in the main manuscript has no such one magnon branch: the spinons at the mean field level are deconfined. Distinguishing between the AKLT state and a nematic quantum spin liquid state is therefore a matter of finding evidence for the one magnon branch of excitations in neutron scattering. If such a signature exists, it will provide strong evidence for the AKLT chain state.

Finally, we should mention that the presence of itinerant fermions likely complicates this story as mentioned in the main text. The survival of the one magnon branch and or even the fundamental distinction between confined and deconfined spinons may disappear though arguments in the literature suggest an FL* state which preserves the fundamental distinction is possible<cit.>. 

