Mechanism for nematic superconductivity in FeSe

TL;DR

This paper proposes a microscopic model for nematic superconductivity in FeSe, explaining key experimental features such as neutron spectra, gap anisotropy, and pairing symmetry through a spin-fermion approach.

Contribution

It introduces a novel microscopic description of nematic quantum spin liquid and its role in superconductivity in FeSe, incorporating orbital-dependent Hund's coupling.

Findings

Reproduces neutron spectra with a broad continuum.

Predicts nodeless $s ext{-}±d$-wave pairing symmetry.

Explains gap anisotropy via orbital-dependent Hund's coupling.

Abstract

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\pm 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.

Figures

a of main text for definition of $J$’s). On a bipartite lattice, it is convenient to perform a unitary transformation by defining $a_{im}\equiv b_{im}$ on A sublattice, and $a_{j\uparrow}\equiv b_{j\downarrow}$, $a_{j\downarrow}\equiv-b_{j\uparrow}$ on B sublattice. The valence bond operator is then brought to the simpler form $A^{\dagger}_{ij}=\sum_{\sigma}a^{\dagger}_{i\sigma}a^{\dagger}_{j\sigma}$. Modular a constant, the spin Hamiltonian ${\cal H}_{S}=\sum_{ij}J_{ij}{\bm{S}}_{i}\cdot{\bm{S}}_{j}$ can be written in terms of the valence bond operators as

a of main text. The corresponding parameters are: $\varepsilon_{\Gamma}=14$, $\delta\varepsilon_{\Gamma}=11$, $\frac{1}{2m_{\Gamma}}=-350$, $b=-70$, $c=120$, $\lambda=9$, $\varepsilon_{1}=-20$, $\varepsilon_{3}=-60$ $\frac{1}{2m_{1}}=75$, $\frac{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

Equations

$${\mathcal{H}}_S=\sum _{ij}{J}_{ij}{\mathbf{S}}_i\cdot {\mathbf{S}}_j,$$

$${\mathcal{S}}_{\mathbf{q},\omega }\sim \sum _{\mathbf{k}}\left\{\cosh \left[2\left({\theta }_{\mathbf{k}}+{\theta }_{\mathbf{k}+\tilde{\mathbf{q}}}\right)\right]-1\right\}\delta \left({\omega }_{\mathbf{k}}+{\omega }_{\mathbf{k}+\tilde{\mathbf{q}}}-|\omega |\right),$$

$${\mathcal{H}}_c=\sum _{\mathbf{k},\alpha \beta ,\nu }{ϵ}_{\alpha \beta }^{\mu \nu }(\mathbf{k})c_{\alpha \mu }^{†}(\mathbf{k})c_{\beta \nu }(\mathbf{k}),$$

$${\mathcal{H}}_{\mathrm{int}}=-\sum _{i,\alpha ,\mu \nu }{J}_{\alpha }{\mathbf{S}}_i\cdot c_{i\alpha \mu }^{†}{\mathbf{\sigma }}_{\mu \nu }{c}_{i\alpha \nu },$$

$${\mathcal{H}}_S=-\frac{1}{2}\sum _{ij}{J}_{ij}{A}_{ij}^{†}{A}_{ij}.$$

$${\mathcal{H}}_S^{(\mathrm{MF})}=\lambda \sum _{i\sigma }{a}_{i\sigma }^{†}{a}_{i\sigma }+\frac{1}{2}\sum _{i\delta \sigma }{Q}_{\delta }\left(a_{i\sigma }^{†}{a}_{i+\delta ,\sigma }^{†}+a_{i\sigma }{a}_{i+\delta ,\sigma }\right).$$

$${\alpha }_{\mathbf{k}\sigma }=\cosh {\theta }_{\mathbf{k}}{a}_{\mathbf{k}\sigma }-\sinh {\theta }_{\mathbf{k}}{a}_{-\mathbf{k}\sigma }^{†},$$

$${\mathcal{H}}_S^{(\mathrm{MF})}=\sum _{\mathbf{k}\sigma }{\omega }_{\mathbf{k}}\left({\alpha }_{\mathbf{k}\sigma }^{†}{\alpha }_{\mathbf{k}\sigma }+\frac{1}{2}\right),$$

$$F=\sum _{\delta }\frac{|Q_{\delta }{|}^2}{2J_{\delta }}-\frac{1}{2}(2S+1)\lambda +\frac{1}{\beta }\int \frac{d^2\mathbf{k}}{(2\pi {)}^2}\ln \left[2\sinh \left(\frac{1}{2}\beta {\omega }_{\mathbf{k}}\right)\right],$$

$$Q{\gamma }_{\mathbf{k}}=2Q_x\cos k_x.$$

$$F^{(1)}=\frac{Q_x^2}{J_1}-\frac{1}{2}(2S+1)\lambda +\frac{1}{\beta }\int \frac{d^2\mathbf{k}}{(2\pi {)}^2}\ln \left[2\sinh \left(\frac{1}{2}\beta {\omega }_{\mathbf{k}}\right)\right],$$

$$S+\frac{1}{2}$$

$$\frac{Q_x}{J_1}$$

$$Q{\gamma }_{\mathbf{k}}=2Q_x\cos k_x+4Q_{x+y}\cos k_x\cos k_y+4Q_{x+2y}\cos k_x\cos (2k_y).$$

$$F^{(2)}=\frac{Q_x^2}{J_1}+\frac{2Q_{x+y}^2}{J_2}+\frac{2Q_{x+2y}^2}{J_4}-\frac{1}{2}(2S+1)\lambda +\frac{1}{\beta }\int \frac{d^2\mathbf{k}}{(2\pi {)}^2}\ln \left[2\sinh \left(\frac{1}{2}\beta {\omega }_{\mathbf{k}}\right)\right],$$

$$S+\frac{1}{2}$$

$$\frac{Q_x}{J_1}$$

$$\frac{Q_{x+y}}{J_2}$$

$$\frac{Q_{x+2y}}{J_4}$$

$${\mathcal{H}}_{0,\Gamma }=\sum _{\mathbf{k}}{\psi }_{\gamma ,\mathbf{k}}^{†}{h}_{\Gamma }(\mathbf{k}){\psi }_{\Gamma ,\mathbf{k}}.$$

$$h_{\Gamma }(\mathbf{k})=\left({\epsilon }_{\Gamma }+\frac{k^2}{2m_{\Gamma }}\right){\tau }^0\otimes {\sigma }^0+\left[\delta {\epsilon }_{\Gamma }+b{k}^2\cos (2{\theta }_k)\right]{\tau }^3\otimes {\sigma }^0+c{k}^2\sin (2{\theta }_k){\tau }^1\otimes {\sigma }^0+\lambda {\tau }^2\otimes {\sigma }^3,$$

$${\mathcal{H}}_{0,X}=\sum _{\mathbf{k}}{\psi }_{X,\mathbf{k}}^{†}{h}_X(\mathbf{k}){\psi }_{X,\mathbf{k}},$$

$$h_X(\mathbf{k})=\left[{\epsilon }_1+\frac{k^2}{2m_1}-a_1{k}^2\cos (2{\theta }_k)\right]\frac{{\tau }^0+{\tau }^3}{2}\otimes {\sigma }^0+\left[{\epsilon }_3+\frac{k^2}{2m_3}-a_3{k}^2\cos (2{\theta }_k)\right]\frac{{\tau }^0-{\tau }^3}{2}\otimes {\sigma }^0+2vk\sin {\theta }_k{\tau }^2\otimes {\sigma }^0.$$

$${\mathcal{H}}_{0,Y}=\sum _{\mathbf{k}}{\psi }_{Y,\mathbf{k}}^{†}{h}_Y(\mathbf{k}){\psi }_{Y,\mathbf{k}},$$

$$h_Y(\mathbf{k})=\left[{\epsilon }_1+\frac{k^2}{2m_1}+a_1{k}^2\cos (2{\theta }_k)\right]\frac{{\tau }^0+{\tau }^3}{2}\otimes {\sigma }^0+\left[{\epsilon }_3+\frac{k^2}{2m_3}+a_3{k}^2\cos (2{\theta }_k)\right]\frac{{\tau }^0-{\tau }^3}{2}\otimes {\sigma }^0+2vk\cos {\theta }_k{\tau }^2\otimes {\sigma }^0.$$

$${\mathcal{D}}^{(\pi ,0)}(\mathbf{q},\Omega )\sim \sum _{\mathbf{k},\omega }{J}_{yz}^2{G}^{(yz)}(\mathbf{k},\omega )G^{(yz)}(\mathbf{k}+\mathbf{q},\omega +\Omega ),$$

$${\mathcal{D}}^{(\pi ,\pi )}(\mathbf{q},\Omega )\sim \sum _{\mathbf{k},\omega }{J}_{xy}^2{G}^{(xy)}(\mathbf{k},\omega )G^{(xy)}(\mathbf{k}+\mathbf{q},\omega +\Omega ),$$

$${\mathcal{S}}_{\mathrm{int}}=-\frac{1}{2}{\int }_0^{\beta }d\tau \sum _{\alpha ,{\alpha }^{\prime }}{J}_{\alpha }{J}_{{\alpha }^{\prime }}{\chi }_{ij}(\tau ){\mathbf{s}}_{i\alpha }(\tau )\cdot {\mathbf{s}}_{j{\alpha }^{\prime }}(0),$$

$${\mathcal{H}}_{\mathrm{pair}}=\sum _{\mathbf{k}{\mathbf{k}}^{\prime }ab}{\Gamma }_{ab}(\mathbf{k},{\mathbf{k}}^{\prime })h_a^{†}(\mathbf{k})h_b({\mathbf{k}}^{\prime }),$$

$$-\sum _j{\oint }_{{\mathrm{FS}}_j}\frac{d{\mathbf{k}}_{\parallel }^{\prime }}{2\pi v_F({\mathbf{k}}^{\prime })}{\Gamma }_{ij}(\mathbf{k},{\mathbf{k}}^{\prime })g_j({\mathbf{k}}^{\prime })=\lambda g_i(\mathbf{k}),$$

$${\chi }^{\prime \prime }(\mathbf{q},\omega )\sim \sum _{\mathbf{k}}\left(1-\frac{{\Delta }_{\mathbf{k}}{\Delta }_{\mathbf{k}+\mathbf{q}}}{E_{\mathbf{k}}{E}_{\mathbf{k}+\mathbf{q}}}\right)\left[1-f(E_{\mathbf{k}})-f(E_{\mathbf{k}+\mathbf{q}})\right]\delta (\omega -E_{\mathbf{k}}-E_{\mathbf{k}+\mathbf{q}}).$$

$${\mathcal{H}}_S=J\sum _i\left[\cos \theta \left({\mathbf{S}}_i\cdot {\mathbf{S}}_{i+1}\right)+\sin \theta {\left({\mathbf{S}}_i\cdot {\mathbf{S}}_{i+1}\right)}^2\right],$$

Full text

FeSe presents new challenges to both perspectives, including explaining its nematic order McQueen et al. (2009)(see Fig. 1(a)), its absence of magnetism, its gapped but active spin fluctuations at $(\pi,\pi)$ in addition to $(\pi,0)$ Wang et al. (2016a) 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 Mukherjee et al. (2015); Kreisel et al. (2015). Renormalization group approaches found the effective interactions promoting spin density wave to be also promoting orbital order Chubukov et al. (2016); Xing et al. (2016); Classen et al. (2016). Approaches focusing on sizable local moments Gretarsson et al. (2011) led to proposals of quadrupolar order accompanying nematic order Yu and Si (2015); Wang et al. (2016b) and the proposal of a quasi-one dimensional quantum paramagnet stateWang et al. (2015) of AKLT (Affleck-Kennedy-Lieb-Tasaki)Affleck et al. (1987) 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 frustrationGlasbrenner et al. (2015); Wang et al. (2015). Indeed FeSe is close to a classic situation for frustrated magnets in the much studied $J_{1}$-$J_{2}$ modelChandra et al. (1990); Misguich and Lhuillier (2004)(see Fig. 1(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 Fernandes et al. (2010). 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) respectivelyHaldane (1988); Read and Sachdev (1990). Moreover DMRG studies on $J_{1}$-$J_{2}$ model noted an intermediate paramagnetic phase between stripe order and Neel order state Jiang et al. (2009, 2012). 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 statesGong et al. (2016). 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 Gretarsson et al. (2011) coexist with itinerant carriers of all three $t_{2g}$ orbitalsNakayama et al. (2014); Watson et al. (2015); Suzuki et al. (2015). In order to capture the dual character Moon et al. (2010) we turn to a spin-fermion model Kou et al. (2009); Lv et al. (2010); Yin et al. (2010); Stanev and Littlewood (2013); Liang et al. (2014); Lee (2015): ${\cal H}={\cal H}_{c}+{\cal H}_{S}+{\cal H}_{\rm int}$, where ${\cal H}_{c}$ and ${\cal H}_{S}$ describe the itinerant carriers and local moments respectively that are coupled through ${\cal H}_{\rm int}$. For the spin model

[TABLE]

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 Fang et al. (2008); Xu et al. (2008). But due to the near itinerancy of the core electrons, longer range terms are also expected Glasbrenner et al. (2015).

The $J_{1}$-$J_{2}$ model has been extensively studied both classically and quantum mechanically (see Refs.[Chandra et al., 1990; Misguich and Lhuillier, 2004; Jiang et al., 2009, 2012]). Within classical models the role of frustration is clear from the fact that the model can be recast as ${\cal H}_{S}=J_{2}\sum({\bm{S}}_{1}+{\bm{S}}_{2}+{\bm{S}}_{3}+{\bm{S}}_{4})^{2}$ up to a constant at $J_{2}=J_{1}/2$ point, where ${\bm{S}}_{1-4}$ are the four spins on each plaquette $\langle 1234\rangle$ 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 boundaryMisguich and Lhuillier (2004) 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 Wang et al. (2016a), and (2) the ordering tendencies obey the classical condition of ${\bm{S}}_{1}+{\bm{S}}_{2}+{\bm{S}}_{3}+{\bm{S}}_{4}=0$ on a plaquette.

A prominent feature of the INS data Wang et al. (2016a) is its broad and gapped continuum of spectral weight (Fig.2a) 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 Balents (2010). Indeed it is a common feature of slave-particle mean field theories. So we will choose Schwinger boson mean field theory (SBMFT) Arovas and Auerbach (1988) as our mean field theory approach. Additional features of Fig. 2 (a-c) we aim to capture include:

• •

The simultaneous presence of both $(\pi,\pi)$ spin fluctuations and $(\pi,0)$, $(0,\pi)$ spin fluctuations.

• •

The quasi-one-dimensional dispersion $\omega\sim\sin k_{y}$ Müller et al. (1981); Lake et al. (2013); Vlijm and Caux (2014) found in the shape of the upper and lower bounds.

• •

The observed cross-shaped spectrum around $(\pi,\pi)$.

To find these features in a SBMFT, we turn to the known Sachdev and Read (1991) SBMFT phase diagram of the $J_{1}$-$J_{2}$ model (Fig. 3). Note that the Neel and stripe long range order for small $J_{2}/J_{1}$ and large $J_{2}/J_{1}$ are expectedRead and Sachdev (1990) to melt into $C_{4}$ symmetric and $C_{2}$ symmetric QSL’s respectively (see Fig. 3). 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. 3 with additional neighbor couplings.

To construct the ansatz, we now turn briefly to the specifics of SBMFT. In Schwinger boson representation, each spin ${\bm{S}}_{\bf r}$ is represented by two bosonic operators $b_{{\bf r}\sigma}$, $\sigma=\uparrow,\downarrow$ with the constraint $\sum_{\sigma}b^{\dagger}_{{\bf r}\sigma}b_{{\bf r}\sigma}=2S$. The spin operator is then ${\bm{S}}_{\bf r}=\frac{1}{2}\sum_{\sigma\sigma^{\prime}}b^{\dagger}_{{\bf r}\sigma}{\bm{\sigma}}_{\sigma\sigma^{\prime}}b_{{\bf r}\sigma^{\prime}}$, with $\bm{\sigma}$ the Pauli matrices. We can then expand $H_{{\bf r},{\bf r^{\prime}}}\equiv J_{{\bf r},{\bf r^{\prime}}}{\bm{S}}_{\bf r}\cdot{\bm{S}}_{\bf r^{\prime}}$ in terms of the spin singlet operator $A^{\dagger}_{{\bf r},{\bf r^{\prime}}}=b^{\dagger}_{{\bf r}\uparrow}b^{\dagger}_{{\bf r^{\prime}}\downarrow}-b^{\dagger}_{{\bf r}\downarrow}b^{\dagger}_{{\bf r^{\prime}}\uparrow}$ to obtain $H_{{\bf r},{\bf r^{\prime}}}=-J_{{\bf r},{\bf r^{\prime}}}\frac{1}{2}A^{\dagger}_{{\bf r},{\bf r^{\prime}}}A_{{\bf r},{\bf r^{\prime}}}+S^{2}$. Finally, we mean-field decompose $H_{{\bf r},{\bf r^{\prime}}}$ and introduce mean fields $\langle A_{{\bm{r}},{\bm{r}^{\prime}}}\rangle$ using $A^{\dagger}_{{\bf r},{\bf r^{\prime}}}A_{{\bf r},{\bf r^{\prime}}}=\langle A^{\dagger}_{{\bf r},{\bf r^{\prime}}}\rangle A_{{\bf r},{\bf r^{\prime}}}+A^{\dagger}_{{\bf r},{\bf r^{\prime}}}\langle A_{{\bf r},{\bf r^{\prime}}}\rangle-\langle A^{\dagger}_{{\bf r},{\bf r^{\prime}}}\rangle\langle A_{{\bf r},{\bf r^{\prime}}}\rangle$. We will further assume the bosons do not condense for we are interested in the spin liquid phase.

Defining $A_{\hat{\mu}}\equiv\langle A_{{\bm{r}},{\bm{r}}+{\hat{\mu}}}\rangle$, we keep $A_{\hat{x}}\neq 0$ and the diagonals $A_{\hat{x}\pm\hat{y}}\neq 0$ and $A_{{\hat{x}}\pm 2{\hat{y}}}\neq 0$ for states in the shaded region of Fig. 3. 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 $(\pi,\pi)$ Néel state and the $(\pi,0)$/$(0,\pi)$ 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. 3. It is a “zero flux state” in that the smallest loop has zero “flux” obeying the so-called flux expulsion principle Tchernyshyov et al. (2006) 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 ${\cal S}_{{\bm{q}}\omega}\equiv{\rm Im}\langle S^{z}({\bm{q}},\omega)S^{z}(-{\bm{q}},\omega)\rangle$ associated with our ansatz. At $T=0$, it is of the form Auerbach and Arovas (1988)

[TABLE]

where $\theta_{\bm{k}}$ is the angle in the Bogoliubov transformation of SBMFT (see SM1 for explicit expression), and ${\tilde{\bm{q}}}={\bm{q}}-(\pi,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. 2(d-f). They capture the basic features of the neutron spectra: (1) The spectrum is gapped (Fig. 2d), as a result of the absence of long range magnetic ordering. (2) Both $(\pi,\pi)$ and $(\pi,0)$/$(0,\pi)$ spin fluctuations are present (Fig. 2d, e). (3) The spectrum displays the novel feature of continuum with the bounds exhibiting quasi-one-dimensional dispersion (Fig. 2d).

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 $(\pi,\pi)$ and $(\pi,0)$, but not at $(0,\pi)$. 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 $(\pi,\pi)$, $(0,\pi)$ and $(\pi,0)$ of roughly the same magnitude even in a single nematic domainMukherjee et al. (2015); Kreisel et al. (2015). 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 ${\bm{q}}$-wavevectors compete, leading to nodes Mukherjee et al. (2015); Kreisel et al. (2015). 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:

[TABLE]

where $c^{\dagger}_{\alpha\mu}({\bm{k}})$ creates an itinerant electron with momentum ${\bm{k}}$, spin $\mu$ and orbital index $\alpha$. The Fermi surface of FeSe consists of two electron pockets around the M points and one hole pocket around the $\Gamma$ point Nakayama et al. (2014); Watson et al. (2015); Suzuki et al. (2015). Following Cvetkovic and Vafek (2013); Fernandes and Chubukov (2016), we take a simple symmetry based approach of expanding the dispersion $\epsilon_{\alpha\beta}^{\mu\nu}({\bm{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 $\Gamma$ point, from $d_{yz}$ and $d_{xy}$ around $(\pi,0)$, from $d_{zx}$ and $d_{xy}$ around $(0,\pi)$. 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 $\Gamma$. The resulting simplified Fermi surface is shown in Fig.4a, see SM2 for explicit parameters. 111Our 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 Lv et al. (2010):

[TABLE]

where ${\bm{\sigma}}$ represents the vector of Pauli matrices, and $J_{\alpha}>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. 222Coupling 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 $\langle A_{{\bm{r}},{\bm{r}}\pm{\hat{x}}}\rangle$ in the nematic QSL state generates an interaction among conduction electrons along the $x$-direction, which drives bond-centered nematic order with $\varphi_{c}\equiv\langle c^{\dagger}_{{\bm{r}}+{\hat{x}},\alpha}c_{{\bm{r}},\alpha}-c^{\dagger}_{{\bm{r}}+{\hat{y}},\alpha}c_{{\bm{r}},\alpha}\rangle\neq 0$ below the temperature at which the nematic QSL develops. The observed nematic transition at $T_{s}\sim 90K$ McQueen et al. (2009) is consistent with this picture. Furthermore, $\varphi_{c}$ linearly couples to $\varphi_{o}\equiv\frac{n_{zx}-n_{yz}}{n_{zx}+n_{yz}}$, where $n_{zx,yz}$ denote occupation of $zx$ and $yz$ orbitals, and $\varphi_{s}\equiv M_{x}^{2}-M_{y}^{2}$, where ${\bm{M}}$ represents the magnetic moment. These different measures of nematicity are consistent with orbital imbalance observed in ARPES Nakayama et al. (2014); Watson et al. (2015); Suzuki et al. (2015) ($\varphi_{o}\neq 0$) and the observed NMR resonance line splitting Baek et al. (2015) ($\varphi_{s}\neq 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. 4b,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,\pi)$ for one nematic domain renders the determination of gap sign in different pockets unfrustrated. In contrast, in the itinerant model where $(\pi,\pi)$, $(\pi,0)$ and $(0,\pi)$ 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 Song et al. (2011); Kasahara et al. (2014a), penetration depth Kasahara et al. (2014a), thermal conductivity measurements Bourgeois-Hope et al. (2016), observing low energy excitations Song et al. (2011); Kasahara et al. (2014a) despite the evidence of a full gap Bourgeois-Hope et al. (2016); Sprau et al. (2016). (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 Wang et al. (2016). More specifically, our gap function is a combination of $d$-wave as induced by $(\pi,\pi)$ spin fluctuations and $s_{\pm}$ as induced by $(\pi,0)$ spin fluctuations. We consider a single nematic domain, where the pairing interaction concentrates around $(\pi,q_{y})$. Two examples of the gap function (in arbitrary units) are shown in Fig.4b,c, where $J_{xy}=J_{yz}=J_{zx}$ and $J_{xy}=0.4J_{yz}=0.4J_{zx}$ respectively.

Now we turn to the question of orbital dependence of the Hund’s coupling. Fig. 4b,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. 4b, $d$- and $s$-wave are at par in Fig. 4c). 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 Yi et al. (2015) implies $J_{xy}\ll 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 measurementsSprau et al. (2016) (see Fig. 4c,d). In Ref.  Sprau et al. (2016) 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 Yi et al. (2015). This orbital dependent Hund’s coupling amplify the role of $(\pi,0)$ spin fluctuation in pairing despite larger spectral weight at $(\pi,\pi)$, which is consistent with the observation of sharp spin resonance at $(\pi,0)$ Wang et al. (2016) (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.Wang et al., 2016a 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 $(\pi,\pi)$ 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 Yi et al. (2015), it is conceivable these systems make better use of already more prominent $(\pi,\pi)$ fluctuation to achieve higher $T_{c}$. We note here that the nematic QSL state we propose is distinct from the proposal of Ref. Wang et al., 2015 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. Senthil et al., 2003, 2004 that has recently been revisited using DMRG Lee et al. (2016).

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.

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}\equiv b_{im}$ on A sublattice, and $a_{j\uparrow}\equiv b_{j\downarrow}$, $a_{j\downarrow}\equiv-b_{j\uparrow}$ on B sublattice. The valence bond operator is then brought to the simpler form $A^{\dagger}_{ij}=\sum_{\sigma}a^{\dagger}_{i\sigma}a^{\dagger}_{j\sigma}$. Modular a constant, the spin Hamiltonian ${\cal H}_{S}=\sum_{ij}J_{ij}{\bm{S}}_{i}\cdot{\bm{S}}_{j}$ can be written in terms of the valence bond operators as

[TABLE]

We then apply mean field theory to the bosonic Hamiltonian Arovas and Auerbach (1988). Defining $Q_{ij}=J_{ij}\langle A_{ij}\rangle\equiv Q_{\delta}$, the quadratic part of the mean field Hamiltonian reads:

[TABLE]

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

[TABLE]

with $\tanh\left(2\theta_{\bm{k}}\right)=-Q\gamma_{\bm{k}}/\lambda$. Here $Q\gamma_{\bm{k}}$ denotes the Fourier transform of $Q_{\delta}$: $Q\gamma_{\bm{k}}=\sum_{\delta}Q_{\delta}e^{-i{\bm{k}}\cdot{\bm{\delta}}}$. The resulting Hamiltonian reads

[TABLE]

with the dispersion $\omega_{\bm{k}}=\sqrt{\lambda^{2}-(Q\gamma_{\bm{k}})^{2}}$. Integrating out the bosonic fields, one obtains the free energy

[TABLE]

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}\neq 0$, and

[TABLE]

Its free energy is

[TABLE]

and the self-consistency equations are

[TABLE]

This ansatz state is basically determined by the dimensionless parameters: $Q_{x}/\lambda$. We have plotted the spin structure factor taking $Q_{x}/\lambda=0.498$ (see Fig. 5a,b,c). Here $S=0.677$.

We consider then our ansatz state, namely the quantum melted stripe state, where $Q_{x}\neq 0$, $Q_{x+y}\neq 0$, $Q_{x+2y}\neq 0$, and

[TABLE]

Its free energy is

[TABLE]

and the self-consistency equations are

[TABLE]

This ansatz state is basically determined by the three dimensionless parameters: $Q_{x}/\lambda$, $Q_{x+y}/Q_{x}$ and $Q_{x+2y}/Q_{x}$. We have plotted the spin structure factor taking $Q_{x}/\lambda=0.398$, $Q_{x+y}/Q_{x}=0.025$ and $Q_{x+2y}/Q_{x}=0.1$ (see Fig. 5d,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 Cvetkovic and Vafek (2013); Fernandes and Chubukov (2016). Consider first the Fermi pocket near the $\Gamma$ point, a single hole pocket has been detected, with the $d_{yz}$ and $d_{zx}$ orbitals dominating the spectral weight. We introduce a spinor $\psi^{T}_{\gamma,{\bm{k}}}=(c_{yz,{\bm{k}}\uparrow},-c_{zx,{\bm{k}}\uparrow},c_{yz,{\bm{k}}\downarrow},-c_{zx,{\bm{k}}\downarrow})$, and the kinetic energy part of the Hamiltonian is of the form

[TABLE]

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 Cvetkovic and Vafek (2013). The result reads

[TABLE]

where $\tau$ and $\sigma$ are Pauli matrices in orbital and spin space respectively, and ${\bm{k}}=(k_{x},k_{y})=k(\cos\theta_{k},\sin\theta_{k})$.

For the electron pocket near $(\pi,0)$, the $d_{yz}$ and $d_{xy}$ orbitals dominate the spectral weight. We introduce a spinor $\psi^{T}_{X,{\bm{k}}}=(c_{yz,{\bm{k}}\uparrow},c_{xy,{\bm{k}}\uparrow},c_{yz,{\bm{k}}\downarrow},c_{xy,{\bm{k}}\downarrow})$, and the kinetic energy part of the Hamiltonian is of the form

[TABLE]

with

[TABLE]

Here ${\bm{k}}$ is measured from $(\pi,0)$.

For the electron pocket near $(0,\pi)$, the $d_{zx}$ and $d_{xy}$ orbitals dominate the spectral weight. We introduce a spinor $\psi^{T}_{Y,{\bm{k}}}=(c_{zx,{\bm{k}}\uparrow},c_{xy,{\bm{k}}\uparrow},c_{zx,{\bm{k}}\downarrow},c_{xy,{\bm{k}}\downarrow})$, and the kinetic energy part of the Hamiltonian is of the form

[TABLE]

with

[TABLE]

Here ${\bm{k}}$ is measured from $(0,\pi)$.

With a proper choice of the parameters, we can obtain a single hole pocket around $\Gamma$, a single electron pocket around $(\pi,0)$, and a single electron pocket around $(0,\pi)$ as shown in Fig.4a of main text. The corresponding parameters are: $\varepsilon_{\Gamma}=14$, $\delta\varepsilon_{\Gamma}=11$, $\frac{1}{2m_{\Gamma}}=-350$, $b=-70$, $c=120$, $\lambda=9$, $\varepsilon_{1}=-20$, $\varepsilon_{3}=-60$ $\frac{1}{2m_{1}}=75$, $\frac{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 $\Gamma$- and M-points, the induced self-energy will be predominantly near ${\bm{q}}=(\pi,\pi)$, $(\pi,0)$ and $(0,\pi)$ (see Fig.6). We expect the neutron spectrum near these points to be smeared. For ${\bm{q}}\sim(\pi,0)$, the self-energy is predominantly from $d_{yz}$ orbitals,

[TABLE]

with fermion Green’s function $G({\bm{k}},\omega)$. Here ${\bm{k}}$ is at the $\Gamma$ pocket, and ${\bm{k}}+{\bm{q}}$ at the $M_{x}$ pocket. For ${\bm{q}}\sim(\pi,\pi)$, the self-energy is predominantly from $d_{xy}$ orbitals,

[TABLE]

where ${\bm{k}}$ and ${\bm{k}}+{\bm{q}}$ are at the two $M$-pockets. With the suppression of $J_{xy}$, we expect the Landau damping effect to be weaker near ${\bm{q}}=(\pi,\pi)$.

We can estimate the strength of the Hund’s coupling from the resulting superconducting transition temperature: $T_{c}\sim E_{F}e^{-1/\lambda}$, with the dimensionless coupling $\lambda=N_{0}V\sim J_{H}^{2}/(E_{F}J_{\rm ex})$. The energy scales involved are: (1) Fermi energy $E_{F}\sim 10$ meV Kasahara et al. (2014b), (2) the exchange interaction $J_{\rm ex}\sim 100$ meV Glasbrenner et al. (2015), and (3) $T_{c}\sim 1$ meV. From these, we obtain Hund’s coupling $J_{H}\sim 20$ meV, which is much smaller than the exchange interaction $J_{\rm 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

[TABLE]

with the itinerant electron spin density ${\bm{s}}_{i\alpha}\equiv\sum_{\mu\nu}c^{\dagger}_{i\alpha\mu}{\bm{\sigma}}_{\mu\nu}c_{i\alpha\nu}$, and the local moment spin correlation function $\chi_{ij}(\tau)\equiv\langle T_{\tau}S_{i}^{a}(\tau)S_{j}^{a}(0)\rangle$. The induced interaction is highly anisotropic, and the dominant interaction term is the nearest-neighbor interaction (say along the $x$-direction): $J_{H}^{2}\chi c^{\dagger}_{{\bm{r}}\alpha}\sigma^{a}_{\alpha\beta}c_{{\bm{r}}\beta}c^{\dagger}_{{\bm{r}}+{\hat{x}},\alpha^{\prime}}\sigma^{a}_{\alpha^{\prime}\beta^{\prime}}c_{{\bm{r}}+{\hat{x}},\beta^{\prime}}$. This interaction results in a phase transition to a nematic state with order parameter $\langle c^{\dagger}_{{\bm{r}}+{\hat{x}},\alpha}c_{{\bm{r}},\alpha}\rangle\neq 0$, or more generally, $\varphi_{c}\equiv\langle c^{\dagger}_{{\bm{r}}+{\hat{x}},\alpha}c_{{\bm{r}},\alpha}-c^{\dagger}_{{\bm{r}}+{\hat{y}},\alpha}c_{{\bm{r}},\alpha}\rangle\neq 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^{\dagger}_{\alpha\alpha^{\prime}}({\bm{k}})=\frac{1}{\sqrt{2}}\left(c^{\dagger}_{{\bm{k}}\alpha\uparrow}c^{\dagger}_{-{\bm{k}}\alpha^{\prime}\downarrow}-c^{\dagger}_{{\bm{k}}\alpha\downarrow}c^{\dagger}_{-{\bm{k}}\alpha^{\prime}\uparrow}\right)$. 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 $\chi({\bm{q}})\equiv\chi({\bm{q}},\Omega_{n}=0)$, which can be obtained from the dynamic spin structure factor via $\chi({\bm{q}})=-\int d\omega{\cal S}_{{\bm{q}},\omega}/\omega$. The special form of $\chi({\bm{q}})$ in FeSe (see Fig. 7) 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 ${\cal H}_{\rm pair}\sim J_{\alpha}^{2}\chi({\bm{k}}-{\bm{k}}^{\prime})h^{\dagger}_{\alpha\alpha}({\bm{k}})h_{\alpha\alpha}({\bm{k}}^{\prime})$.

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^{\dagger}_{{\bm{k}}\alpha\mu}=\sum_{a}\eta^{*}_{\alpha a\mu}({\bm{k}})d^{\dagger}_{{\bm{k}}a\mu}$ with band index $a$. Note that since the spin-orbit coupling here is in the $\sigma^{3}$ channel, different spins do not mix. The pair operator in the band basis is $h^{\dagger}_{a}({\bm{k}})=\frac{1}{\sqrt{2}}\left(d^{\dagger}_{{\bm{k}}a\uparrow}d^{\dagger}_{-{\bm{k}}a\downarrow}-d^{\dagger}_{{\bm{k}}a\downarrow}d^{\dagger}_{-{\bm{k}}a\uparrow}\right)$. Omitting the frequency dependence, the pairing Hamiltonian is of the form

[TABLE]

with the projected pairing interaction $\Gamma_{ab}({\bm{k}},{\bm{k}}^{\prime})=\frac{1}{2}\sum_{\alpha}J_{\alpha}^{2}\chi({\bm{k}}-{\bm{k}}^{\prime})M^{*}_{\alpha a}({\bm{k}})M_{\alpha b}({\bm{k}}^{\prime})$. The orbital content is encoded in the form factor $M_{\alpha a}({\bm{k}})=\eta_{\alpha a\uparrow}({\bm{k}})\eta_{\alpha a\downarrow}(-{\bm{k}})$. The gap function is then defined as $\Delta_{a}({\bm{k}})=\sum_{{\bm{k}}^{\prime}b}\Gamma_{ab}({\bm{k}},{\bm{k}}^{\prime})\langle h_{b}({\bm{k}}^{\prime})\rangle$. The gap symmetry function $g_{i}({\bm{k}})\propto\Delta_{a}({\bm{k}})$ on the Fermi surface is determined by the eigen equation

[TABLE]

where ${\bm{k}}_{\parallel}$ denotes momentum along the Fermi surface ${\rm FS}_{j}$, and $v_{F}({\bm{k}})=|\nabla E_{a}({\bm{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 $\epsilon_{\alpha\beta}^{\mu\nu}({\bm{k}})$ (2) local moment spin susceptibility $\chi({\bm{q}})$, and (3) Hund’s couplings $J_{\alpha}$.

We show here more concretely how $(\pi,\pi)$ spin fluctuation mediated pairing enhances the superconducting $T_{c}$. We first estimate the value of $\lambda$ 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}\sim E_{F}e^{-1/\lambda}$ Kirzhnits et al. (1973); Klimin et al. (2014). With $T_{c}\sim 8$ K, $E_{F}\sim 10$ meV Kasahara et al. (2014b), we obtain $\lambda\sim 0.37$. In bulk FeSe, pairing occurs predominantly among $d_{yz}$ orbitals as mediated by $(\pi,0)$ spin fluctuations, while $(\pi,\pi)$ spin fluctuation mediated pairing among $d_{xy}$ orbitals is largely suppressed. This corresponds to taking $(J_{xy},J_{zx},J_{yz})\sim(0,1,1)$. (Note that due to the near absence of $(0,\pi)$ 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})\sim(1,1,0)$. We have obtained the resulting eigenvalue $\lambda^{\prime}=2.98\lambda$, which gives $T_{c}\sim 47$ K. Hence $(\pi,\pi)$ spin fluctuation mediated pairing is indeed able to account for the much higher $T_{c}$ in heavily doped FeSe ($T_{c}\sim 48$ K Miyata et al. (2015)), and a large part of the $T_{c}$ increase in monolayer FeSe ($T_{c}\sim 50-64$ K Lee et al. (2014)).

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 formBulut et al. (1992)

[TABLE]

The gap function changes sign between the $\Gamma$ pocket and $(\pi,0)$ pocket, and between the $(\pi,0)$ pocket and $(0,\pi)$ pocket. Since the sizes of the Fermi pockets are small, one expects the resonances to be localized in momentum space around ${\bm{q}}=(\pi,0)$ and ${\bm{q}}=(\pi,\pi)$.

Due to lack of full knowledge of the gap function in the whole Brillouin zone, we approximate the neutron intensity at ${\bm{q}}=(\pi,0)$ and ${\bm{q}}=(\pi,\pi)$ 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 ${\bm{q}}$, i.e. we replace $\sum_{\bm{k}}$ by $\oint_{\bm{k}}\oint_{{\bm{k}}^{\prime}}e^{-|{\bm{k}}-{\bm{k}}^{\prime}-{\bm{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.8(a) (with $q_{0}=0.2$). One can see that the resonance at ${\bm{q}}=(\pi,\pi)$ is suppressed due to worse nesting of the corresponding Fermi surface and gap functions.

Actually one expects further suppression of $(\pi,\pi)$ resonance with more realistic band structures. In particular, it has been found in quantum oscillation measurements Terashima et al. (2014) that the Fermi pocket at $(0,\pi)$ (with $k_{F}\simeq 0.13{\rm\AA}^{-1}$) is much larger than the Fermi pocket at $(\pi,0)$ (with $k_{F}\simeq 0.043{\rm\AA}^{-1}$). Indeed a simple check by enlarging the Fermi pocket at $(0,\pi)$ by a factor of two while keeping the rest of the band structure and gap functions fixed suppreses the resonance at ${\bm{q}}=(\pi,\pi)$ to be almost vanishing (see Fig.8(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 Read and Sachdev (1990) 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. Wang et al., 2015 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

[TABLE]

with $-\pi\leq\theta\leq\pi$. As shown in Fig.9, its phase diagram contains (see Läuchli et al. (2006); Manmana et al. (2011) and references therein): (1) a ferromagnetic phase for $-\pi<\theta<-3\pi/4$ and $\pi/2<\theta\leq\pi$, (2) a dimerized phase for $-3\pi/4<\theta<-\pi/4$, (3) a gapped and topologically ordered Haldane phase for $-\pi/4<\theta<\pi/4$ (the AKLT state corresponds to $\theta=\arctan\frac{1}{3}\simeq 0.1024\pi$) and (4) a gapless phase with antiferroquadrupolar (AFQ) correlations for $\pi/4<\theta<\pi/2$ Läuchli et al. (2006). The dimerized phase and the Haldane phase are separated by a critical point, the Takhtajan-Babujian (TB) point Takhtajan (1982); Babujian (1982). At the TB point, the system possesses gapless spinon excitations. The spinon continuum is manifest in the dynamical spin structure factor (Fig.10) as obtained using algebraic Bethe ansatz-based method Vlijm and Caux (2014). 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.11) Schmitt et al. (1998). Returning back to the critical point, the two-magnon excitations merge with the one magnon excitations and only a continuum of spinon excitations remainSchmitt et al. (1998).

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. 11. 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 possibleSenthil et al. (2003, 2004); Lee et al. (2016).