Inter-orbital nematicity and the origin of a single electron Fermi pocket in FeSe

The research of electronic properties of iron-based superconductors has proven to be a rich field that challenges our understanding of correlated multi-band systems. In this respect, particularly the superconducting iron-chalcogenides, FeSe and doped FeTe, have attracted considerable attention due to their unusual low-energy electronic states[1, 2, 3]. FeSe enters an orthorhombic phase near $T_{n}\sim 90$ K without concomitant static magnetic order at lower temperatures, unlike most of the iron-pnictide superconductors. There is strong evidence that the rotational symmetry breaking at $T_{n}$ is driven by the electronic degrees of freedom, and thus constitutes a rare example of electronic nematicity[4]. In particular, the spectrum of low-energy spin excitations, thought to drive electron pairing, is extremely anisotropic in untwinned crystals[5]. The superconducting phase sets in around $T_{c}\sim 9$ K and is therefore generated from an instability of the nematic “normal” state. Thus, it is not surprising that the superconducting properties also break rotational symmetry, as is indeed observed experimentally[6]. However, the origins of the extreme normal state nematicity, as well as the reasons for the absence of magnetism at ambient pressure, are still being debated[1, 2, 3].

The very high level of rotational anisotropy is evident, e.g., from theoretical modelling of experimental data measured at low temperatures in FeSe[3]. Quite generally, the anisotropy inherent to standard electronic bands in the orthorhombic state, is not enough to explain the much more extreme anisotropy observed experimentally[7, 8, 6, 9, 5, 10]. Thus, various interaction-driven feedback effects have been proposed that enhance the symmetry breaking[6, 11, 12, 13]. In this respect, several works have argued for orbital-dependent quasiparticle weights $Z_{\alpha}$ and shown how this effect may reconcile several experimental probes with simple models[14, 6, 9, 11, 15, 16, 17]. The existence of orbital-dependent $Z_{\alpha}$ has been explored extensively within dynamical mean-field theory (DMFT) and related methods, and arise naturally in multi-orbital systems with different degrees of orbital contributions to the Fermi level electronic states[18, 14, 19, 20, 21, 22]. For the specific case of FeSe, it was shown that $Z_{xy}$ must be strongly reduced from 1, but additionally that a large quasiparticle weight anisotropy of $Z_{xz}/Z_{yz}\sim 0.2-0.3$ was necessary in order to obtain agreement with experiments[6, 9, 23] within this approach. The origin of this rather small ratio remains unsettled, but may be related to strong feedback effects on the electronic states close to a stripe magnetic quantum critical point. One of the important implications of these extracted values for the orbital weights $Z_{\alpha}$ was that the electron pocket at $Y$ should be nearly completely incoherent, and thus difficult to observe with spectroscopic probes; indeed, no conclusive observation of this pocket has been reported by either angular-resolved photoemission spectroscopy (ARPES)[24, 25, 26] or scanning tunneling microscopy (STM)[9].

Subsequently, however, laser ARPES studies reported significant $d_{yz}$ content on the $\Gamma$-centered hole pocket, as well as a strong $d_{xz}$ component apparently inconsistent with the extreme “orbitally selective” scenario[27]. Furthermore, very recent synchrotron ARPES measurements on detwinned FeSe have proposed that the $Y$ pocket may be lifted entirely from the Fermi surface rather than unobservable due to incoherence[28, 29]. The explanation of the unusual spectroscopic observations of this pocket is therefore important not only to obtain the correct low-energy electronic band structure relevant for (detwinned) nematic FeSe, but also to ultimately understand the origin of the strong electronic nematicity in FeSe.

Of course, any theoretical approach for explaining the large electronic anisotropy of FeSe relies on having a reasonable starting point for inclusion of interactions in the band structure. It was initially proposed, based on comparison to ARPES data available at the time, that FeSe features nematicity in the $d_{xz}$ and $d_{yz}$ intra-orbital channels with form factors transforming as the irreducible representations (IRs) A_1g and B_1g of the given point group[30, 31, 32, 33, 34, 35, 36, 37]. More recently, however, several works have advocated for the relevance of (intra-orbital) nematicity involving also the $d_{xy}$ orbitals[38, 39, 40, 41, 42, 43]. Until recently, first principles approaches were unable to stabilize a nematic ground state without concomitant stripe magnetism. However a recent density functional theory DFT+U exploration of the energy landscape found a lowest energy nematic state transforming according to the E_u irreducible representation of the D_4h point group, containing significant inter-orbital nematic components[44]. This form of nematicity, which breaks also inversion symmetry, was shown to produce Fermi surfaces containing only a single electron pocket[44]. Whether this approach has material-specific predictive power is still an open question¹¹1 Applying a similar approach to LiFeAs, for example, produces a dramatically distorted low-energy band structure with missing $d_{xy}$ band at the Fermi level. R. Valenti, private communication., however, and in any case it remains important to identify the cause of nematic ordering using more transparent model-based methods.

These recent developments raise a number of important questions related to the low-energy band structure and the associated superconducting gap structure of FeSe. For example, what underlying physical interactions naturally produce nematic order that generate single-electron-pocket bands[44, 43], and what are the consequences of this low-energy band structure for our broader understanding of this fascinating material?

Here, motivated by several earlier theoretical studies of the effects of longer-range Coulomb interactions $V$ on band structure renormalizations of iron-based superconductors[38, 46, 47, 48], we explore the nematic phase spontaneously generated by $V$. Ref. 38 originally investigated these effects for FeSe, and additionally found a dominant intra-orbital $d$-wave bond nematic order arising from nearest neighbor Coulomb repulsion. These results are, however, inconsistent with the current ARPES description of the FeSe Fermi surface and deserve to be re-examined. To this end, we have begun with the canonical DFT calculation of the FeSe band structure in the tetragonal phase[49], and introduced nematic order driven by longer-range Coulomb interactions systematically. We find that, in addition to the generation of small electron and hole pockets by $V$, it generates inter-orbital nematicity between $d_{yz}$ and $d_{xy}$, and $d_{xz}$ and $d_{xy}$ states in a large region of parameter space at low temperatures. This form of nematic order is distinct from those discussed previously in the literature[38, 44, 42, 43]. The inter-orbital nematic order components hybridize the low-energy bands near $X$ and $Y$, respectively, and thereby allow, depending on their amplitudes, for the lifting of one of the electron pockets. We explore this “$V$-scenario” for nematicity and Fermi surface anisotropy, and demonstrate how it naturally generates Fermi surfaces containing only a single electron pocket. We next discuss the comparison of the spin excitations obtained using this renormalized band structure with the highly anisotropic spectrum reported in the experimental literature[5]. Finally, we compute the resulting momentum-dependent superconducting gap structure, and discuss implications for other experimental probes of FeSe.

In order to explore how the electronic structure of FeSe is affected by longer-range Coulomb interactions, specifically nearest-neighbor (NN) density interactions, we apply the following many-body Hamiltonian

$\displaystyle\begin{split}H=&\,H_{\rm kin}+H_{\rm soc}+H_{\rm int}\\ =&-\sum_{i\,j}\sum_{\mu\nu}\sum_{\sigma}c^{\,\dagger}_{i\mu\sigma}\left(t^{\,\mu\nu}_{ij}+\mu_{0}\delta_{ij}\delta_{\mu\nu}\right)c^{\,\phantom{\dagger}}_{j\nu\sigma}\\ &+\lambda_{\rm SOC}\sum_{i}\sum_{\mu\nu}\sum_{\sigma\sigma^{\prime}}c^{\dagger}_{i\mu\sigma}\,\bm{L}^{\mu\nu}\cdot\bm{S}^{\sigma\sigma^{\prime}}c^{\phantom{\dagger}}_{i\nu\sigma^{\prime}}\\ &+\frac{V}{2}\sum_{\langle i,j\rangle}\sum_{\mu\nu}\sum_{\sigma\sigma^{\prime}}n^{\phantom{\dagger}}_{i\mu\sigma}n^{\phantom{\dagger}}_{j\nu\sigma^{\prime}},\end{split}$

(1)

where $\langle i,j\rangle$ indicates the set of NN sites. Here, the operator $c^{\phantom{\dagger}}_{i\mu\sigma}$ annihilates an electron in orbital $d_{\mu}$ with spin projection $\sigma$ at lattice site $\bm{R}_{i}$, and $n_{i\mu\sigma}$ represents the density operator. The kinetic part of the Hamiltonian[49] $H_{\rm kin}$ contains the hopping matrix elements and the chemical potential $\mu_{0}$, $H_{\rm soc}$ includes the atomic spin-orbit coupling (SOC) with strength $\lambda_{\rm SOC}$, and $H_{\rm int}$ describes the repulsive NN density interaction. For all results presented in the present work, we adjust $\mu_{0}$ to keep a fixed electron filling of $\langle n\rangle=6.0$. Note that for simplicity we completely discard the usual onsite Hubbard Coulomb repulsion $H_{U}$ in the model. As is well-known such interactions can lead to important band renormalization, magnetism, and nematicity[4, 50, 3, 51, 52], none of which lift the $Y$ electron pocket. Here, we assume that $H_{U}$ merely renormalizes the DFT band structure, though we do not include such effects explicitly, and explore the nematic instability generated solely from $H_{\rm int}$ containing only NN Coulomb repulsion.

The symmetry properties of FeSe, and thereby $H$, are governed by the non-symmorphic space group $P4/nmm$, which consists of eight point group elements $\{g\,|\,(0,0,0)\}$ for $g\in$ D_2d, and $\{gI\,|\,(\nicefrac{{1}}{{2}},\nicefrac{{1}}{{2}},0)\}$ with $I$ denoting inversion. In order for the above to constitute a closed group, the product of space group elements is defined modulo integer lattice translations [53]. For both analytical tractability and numerical simplicity we perform our study in the 1-Fe unit cell, for which the space group $P4/nmm$ is lowered to the symmorphic subgroup D_2d. Although one thereby lacks certain symmetry elements of the original system, it is still possible to capture the violation of rotational symmetry and the emergence of nematic order. In particular, we need to focus only on the generators of the $\rm D_{2d}$ group $g=\{E,\,S_{4},\,C_{2},\,C^{\prime}_{2},\,\sigma_{d}\}$, and the associated five IRs $\Gamma=\{\rm A_{1},\,A_{2},\,B_{1},\,B_{2},\,E\}$. Here, the generator $S_{4}$ refers to the combined operation of $C_{4}\sigma_{xy}$.

We incorporate the effects of NN Coulomb interactions by performing a Hartree-Fock mean-field (MF) decoupling of $H_{\rm int}$, and introduce the following homogeneous bond-order fields in wave-vector space [38, 46]

$\displaystyle N^{\mu\nu}_{\bf{k}\sigma}=-\frac{2}{\mathcal{N}}\sum_{\bf{k}^{\prime}}f_{\bf{k}-\bf{k}^{\prime}}^{\rm A_{1}}\langle c^{\dagger}_{\bf{k}^{\prime}\nu\sigma}c^{\phantom{\dagger}}_{\bf{k}^{\prime}\mu\sigma}\rangle,$

(2)

with $\mathcal{N}$ being the number of lattice sites, and $f_{\bf{k}}^{\rm A_{1}}$ being the form factor of the interaction, which transforms as the IR $\rm A_{1}$. Specifically for the nearest-neighbor interaction studied here $f_{\bf{k}}^{\rm A_{1}}\propto\sum_{{\rm d}}\cos(k_{\rm d}\,a_{\rm d})$, where $\rm d$ labels the set of primitive lattice vectors, and $a_{\rm d}$ is the lattice constant in the $\rm d$-direction. Note that in the decoupling of $H_{\rm int}$, the Hartree terms contribute only to a constant energy shift, which can be readily absorbed in the chemical potential $\mu_{0}\mapsto\mu_{0}-4V\sum_{\bf{k}^{\prime}\mu\sigma}\langle n_{\bf{k}^{\prime}\mu\sigma}\rangle/\mathcal{N}$.

Upon reaching the nematic transition temperature $T_{n}$, the system will spontaneously violate $S_{4}$ symmetry, and we therefore seek to split the bond order fields into $S_{4}$ symmetry-preserving and symmetry-breaking terms, in order to identify the origin of nematicity driven by NN interactions. The symmetry-preserving terms, denoted by $N_{{\bf k}\sigma,{\rm br}}^{\mu\nu}$, will in general lead to band renormalizing (br) terms, and is obtained by averaging over the four cycles of $S_{4}$

$\displaystyle N_{\bf{k}\sigma,{\rm br}}^{\mu\nu}=\frac{1}{4}\sum_{\ell=0}^{3}\sum_{\mu^{\prime}\nu^{\prime}}\big{[}D(S^{\ell}_{4})\big{]}^{\mu\mu^{\prime}}\,N^{\mu^{\prime}\nu^{\prime}}_{S_{4}^{-\ell}\bf{k}\,\sigma}\,\big{[}D^{\dagger}(S^{\ell}_{4})\big{]}^{\nu^{\prime}\nu},$

(3)

where $D(S_{4}^{\ell})$ is the representation of $S_{4}^{\ell}$ in orbital space. As a consequence of SOC, one should in general also perform the $S_{4}$ averaging in spin space. However, by exploiting that the fields $N_{{\bf k}\sigma}^{\mu\nu}$ are diagonal in spin space and the fact that the action of $S_{4}$ is equivalent to a $\pi/2$-rotation about the spin $z$-axis, we find that $S_{4}$ leaves the spin projection $\sigma$ of $N_{{\bf k}\sigma}^{\mu\nu}$ invariant. Momentum dependent band renormalizing terms, such as $N_{{\bf k}\sigma,{\rm br}}^{\mu\nu}$, have been previously put forward as potential candidates for explaining the band structure renormalizations of iron-pnictides and iron-chalcogenides[38, 46, 51].

Here, the main focus is on the $S_{4}$ symmetry breaking (sb) terms leading to nematic order, which are simply the remainder of the total field $N_{\bf{k}\sigma,{\rm sb}}^{\mu\nu}=N_{\bf{k}\sigma}^{\mu\nu}-N_{\bf{k}\sigma,{\rm br}}^{\mu\nu}$. In order to extract the form factors entering in the symmetry-preserving and symmetry-breaking fields, we perform the following projection onto the normalized lattice versions of the basis functions $f_{\bf{k}}^{\gamma}$ belonging to the group $\rm D_{2d}$

$\displaystyle N_{\bf{k}\sigma,{\rm br/sb}}^{\mu\nu}=\sum_{\gamma}f^{\gamma}_{\bf{k}}N_{\gamma\sigma,{\rm br/sb}}^{\mu\nu}.$

(4)

As a consequence of the NN interaction, the sum over $\gamma$ is restricted to include only the lowest order lattice harmonics, i.e. only linear combinations of $\cos(k_{\rm d}a_{\rm d})$ and $\sin(k_{\rm d}a_{\rm d})$ enter in $f_{\bf{k}}^{\gamma}$, while a longer-ranging interaction in general allows for higher order terms. Thus, we arrive at the following MF decoupled interaction

$\displaystyle H^{\rm MF}_{\rm int}=V\sum_{{\bf{k}},\gamma}\sum_{\mu\nu}\sum_{\sigma}c^{\,\dagger}_{\bf{k}\mu\sigma}f_{\bf{k}}^{\gamma}\left(N_{\gamma\sigma,{\rm br}}^{\mu\nu}+\alpha N_{\gamma\sigma,{\rm sb}}^{\mu\nu}\right)c^{\,\phantom{\dagger}}_{\bf{k}\nu\sigma}.$

(5)

Note that we here introduced the enhancement factor $\alpha$, since the interaction strength $V$ in the symmetry-preserving and symmetry-breaking channels in general will be different, since the fields belong to different IRs of the point group, and can potentially lead to distinct couplings in a renormalization group procedure. Thus we allow $\alpha\neq 0$ throughout, similar to earlier works[38, 46].

Lastly, in App. A we elaborate further on the construction of $N_{\gamma\sigma,{\rm br/sb}}^{\mu\nu}$, display their explicit matrix structure, and list the basis functions $f_{\bf{k}}^{\gamma}$ used in the upcoming sections.

Equipped with the above model Hamiltonian and symmetry considerations, we are now able to calculate the bond order fields $N_{\gamma\sigma,{\rm br/sb}}^{\mu\nu}$ self-consistently for a fixed electron filling $\langle n\rangle=6.0$, (for details see App. A). In our calculations we adopt the kinetic Hamiltonian for FeSe given in Ref. [49], after including SOC in the spin $z$-direction of bare strength $\lambda_{\rm SOC}=30\,$meV. However, the exact value of SOC is not qualitatively important for the emergence of inter-orbital nematic components and the concomitant disappearance of one of the electron pockets, as found further below. However, SOC is important, within the current model, for the FS to contain only a single hole pocket at $\Gamma$. Due to weak inter-layer coupling in the $c$ direction, we restrict our study to the quasi-2D lattice of a single layer of Fe and Se atoms. The resulting orbitally-resolved Fermi surface (FS) and band structure in the absence of NN interactions are displayed in Fig. 1(a) and 1(c), respectively.

The effects of NN interactions, captured by the terms entering in Eq. (5), give rise to the two aforementioned effects, namely band renormalization and $S_{4}$ symmetry breaking. Let us momentarily discuss the former, $N_{\gamma\sigma,{\rm br}}^{\mu\nu}$, by focusing on the low-energy $t_{2g}$ orbitals. For this case, the fields lead to the following two effects: $i)$ collective down-shift of the hole pockets at $\Gamma$ and $M$, and $ii)$ up-shift of the Dirac points at $X$ and $Y$. The down-shift of the hole pockets is readily attributed to the self-consistent fields in Fig. 2(a), which couple to the form factor $f_{\bf{k}}^{\rm\,s}=\cos k_{x}+\cos k_{y}$. Specifically the fields $N_{{\rm s}\sigma,{\rm br}}^{xz\,xz}=N_{{\rm s}\sigma,{\rm br}}^{yz\,yz}=-0.22$ lead to a down-shift of $d_{xz}$,$d_{yz}$-dominated bands at $\Gamma$, while $N^{xy\,xy}_{{\rm s}\sigma,{\rm br}}=0.24$ ensures the down-shift of $d_{xy}$-dominated bands at $M$. The latter occurs since the form factor has a relative sign on the two pockets $f_{\Gamma}^{\,\rm s}=-f_{M}^{\,\rm s}$. Additionally we note that the effects of the $\rm s$-wave fields on the electron pockets are minimal since $f^{\,\rm s}_{X}=f^{\,\rm s}_{Y}=0$.

By contrast to the above, the fields presented in Fig. 2(b) are dominant at the $X$ and $Y$ points due to the $\rm d$-wave form factor $f_{\bf{k}}^{\,\rm d}=\cos k_{x}-\cos k_{y}$. In fact, the field $N_{{\rm d}\sigma,{\rm br}}^{xz\,xz}=0.09$ ($N_{{\rm d}\sigma,{\rm br}}^{yz\,yz}=-0.09$) shifts $d_{xz}$($d_{yz}$)-dominated bands at $Y$ ($X$), pushing the Dirac points up closer to the Fermi level. Similar to the $\rm s$-wave fields, also here the symmetry of the form factor, $f^{\,\rm d}_{X}=-f^{\,\rm d}_{Y}$, compensates the relative sign between the two fields $N^{xz\,xz}_{\rm d\sigma,br}/N^{yz\,yz}_{\rm d\sigma,br}=-1$. Furthermore, we point out that the electron pockets acquire a peanut-like shape for appropriate values of $V$, since the fields $N_{\rm d\sigma,br}^{\mu\nu}$ do not affect the $d_{xy}$ orbitals. Lastly note that the down- (up-) shift of the hole pockets (Dirac points) leads to smaller $\Gamma$ and $M$ ($X$ and $Y$) pockets. In fact, it was previously shown that $N_{\bf{k}\sigma,{\rm br}}^{\mu\nu}$ can completely remove the $M$ pocket in FeSe, while for LiFeAs it was the $\Gamma$ pockets which were pushed away from the Fermi level [51].

Turning our attention to the symmetry breaking fields, $N_{\gamma\sigma,{\rm sb}}^{\mu\nu}$, we find that only fields that are coupled to the $\rm d$-wave form factor play an essential role, see Fig. 2(c) and 2(d). By applying the same logic as for the band renormalizing terms, we observe that the field $N_{\rm d\sigma,sb}^{xz\,xz}$ ($N_{\rm d\sigma,sb}^{yz\,yz}$) leads to an up- (down-) shift of $d_{xz}$($d_{yz}$)-dominated bands at $Y$ ($X$). While this effect appears similar to the one encountered for $N_{\rm d\sigma,br}^{\mu\nu}$, we stress that here the Dirac points are shifted in opposite directions due to the violation of $S_{4}$-symmetry. This reduction in symmetry also allows for anisotropic inter-orbital $d_{xz}-d_{xy}$ and $d_{yz}-d_{xy}$ hybridization terms, namely $\alpha N^{xz\,xy}_{\rm d\sigma,sb}\approx 0.10$ and a vanishing coupling $\alpha N^{yz\,xy}_{\rm d\sigma,sb}$. It is in fact these particular couplings which lead to the distinct hybridization gaps at $X$ and $Y$ evident from Fig. 1(d), and they will, in synergy with all the effects discussed above, result in the highly anisotropic FS shown in Fig. 1(b), featuring only a single electron pocket at the $X$ point. Thus, as a function of temperature, a Lifshitz transition necessarily takes place such that a shrinking $Y$-pocket eventually disappears at a certain temperature below $T_{n}$. A detailed discussion of the experimental evidence for such a temperature-induced Lifshitz transition can be found in Ref. 43.

To gain deeper insight into the emergent nematic order, we now proceed and study more carefully the IRs of $N_{{\bf k}\sigma,\rm br/sb}^{\mu\nu}$. By construction, the $S_{4}$ symmetry preserving fields can only consist of terms transforming as A₁ and A₂, and, in fact, we find through our self-consistent calculations that only A₁ terms are non-zero in $N_{{\bf k}\sigma,{\rm br}}^{\mu\nu}$, see App. B for further details.

In contrast, $N_{{\bf k}\sigma,{\rm sb}}^{\mu\nu}$ can contain terms transforming as the remaining three IRs of the group, i.e. $\rm B_{1}$, $\rm B_{2}$ and $\rm E$, allowing for the nematic order to arise in any of these three channels. As we explicitly show in App. B, we find non-zero $S_{4}$ symmetry breaking terms belonging to two distinct IRs, namely B₁ and E. This implies that the system undergoes two consecutive phase transitions upon lowering of the temperature. In order to quantify this, we focus on the $t_{2g}$ orbitals, and define the following leading order parameters

$\displaystyle\begin{split}N_{\rm B_{1}}&=(N_{\rm d\sigma,sb}^{xz\,xz}+N_{\rm d\bar{\sigma},sb}^{xz\,xz}+N_{\rm d\sigma,sb}^{yz\,yz}+N_{\rm d\bar{\sigma},sb}^{yz\,yz})/4,\\ N_{\rm E}&=\left([N_{{\rm d}\sigma,{\rm sb}}^{xz\,xy}+N_{{\rm d}\bar{\sigma},{\rm sb}}^{xz\,xy}]/2,[N_{{\rm d}\sigma,{\rm sb}}^{yz\,xy}+N_{{\rm d}\bar{\sigma},{\rm sb}}^{yz\,xy}]/2\right)\\ &\equiv(N_{{\rm E}_{x}},N_{{\rm E}_{y}}),\end{split}$

(6)

where $\bar{\sigma}$ is the opposite spin projection of $\sigma$. For details see App. B. In Fig. 3 we show the values of these order parameters for various temperatures, and indeed find that $N_{\rm B_{1}}$ becomes non-zero at $T_{n}$, while only the single component $N_{{\rm E}_{x}}$ condenses at lower temperatures $T^{\prime}_{n}$. An intra-orbital nematic order thus arises at $T_{n}$, lowering the point group symmetry from $\rm D_{\rm 2d}$ to $\rm D_{2}$, and leads to the up- (down-) shift of $d_{xz}$($d_{yz}$)-dominated bands at $Y$ ($X$) discussed in the previous paragraph. The second transition at $T_{n}^{\prime}$ further reduces the point group to $\rm C_{2}$ and allows for anisotropic inter-orbital $d_{xz}-d_{xy}$ or $d_{yz}-d_{xy}$ hybridization terms, i.e. the effects ultimately leading to the highly anisotropic FS, shown in Fig. 1(b).

We stress that the effects of the $e_{g}$ orbitals and remaining fields $N^{\mu\nu}_{\gamma\sigma,{\rm br/sb}}$ not displayed in Fig. 2, are all incorporated in our calculations, but do not lead to qualitative changes of the low-energy band structure and FS, and are therefore not explicitly mentioned in the above discussion. For a complete overview of all fields $N_{\gamma\sigma,{\rm br/sb}}^{\mu\nu}$, see App. A. Furthermore, we note that the final nematic band structure, as seen e.g. in Fig. 1(d), obviously depends on the starting point, i.e. the tetragonal DFT band structure. Thus, the final quantitative energy scales, i.e. the hybridization gap at $Y$ and the required amplitudes of $V$ and $\alpha$ for generating a FS similar to that shown in Fig. 1(b), depend on the initial bare band structure. Lastly, we note that the purely intra-orbital nematic order generated from NN Coulomb repulsion found earlier [38, 46], can be reproduced here when applying the same band structure as in Ref. 38. We have not located the exact band property that leads to the additional inter-orbital nematic fields from the DFT bulk FeSe band model used in the current study[49], but note that the solution presented in Fig. 1 is quite generic at low temperatures and exists for a wide parameter range.

The low-energy electronic structure established above has important consequences for spin excitations that should be compared with experiments. In this respect, a series of inelastic neutron scattering experiments[54, 55, 5] have established a remarkable set of magnetic phenomena in FeSe. At low temperatures in the nematic phase, but high energies $\sim$ 100 meV, strong $(\pi,\pi)$ (Néel) fluctuations dominate the spectrum, with weaker but still prominent $(\pi,0)$ and $(0,\pi)$ (stripe-like) fluctuations. As the energy is lowered to a few tens of meV, a spin gap develops in the $(\pi,\pi)$ spectrum, but $(\pi,0)$ spin fluctuations strengthen. Notably, the only measurement of detwinned FeSe crystal finds that at low energies the intensity of the $(0,\pi)$ fluctuations essentially vanishes[5]. This extraordinary result should emerge from a proper theory of low energy spin and orbital degrees of freedom. We show now that, counter-intuitively, the current proposal for low-energy nematic electronic structure is not sufficient to explain the above-mentioned magnetic properties, and requires the additional physics of orbital selective correlations.

In Fig. 4(a), we first illustrate the bare magnetic susceptibility Re$\chi_{0}({\bf q},\omega=0)$, together with the enhanced susceptibility obtained in the random phase approximation $\chi^{RPA}$. The spectrum is clearly dominated by intense $(\pi,\pi)$ fluctuations arising from scattering between the $d_{xy}$ states, which are not observed in experiment. Furthermore, $(\pi,0)$ and $(0,\pi)$ states are nearly degenerate, in contradiction to the results of Ref. 5. Note that these issues are also common to conventional spin-fluctuation theories of FeSe including a $Y$ pocket[23, 56], or other novel schemes to lift the $Y$ pocket[43] via nematic order. In Ref. 23, it was proposed that they could be resolved by assuming orbitally selective incoherence of $d_{xy}$, $d_{xz}$, and $d_{yz}$ states, such as should take place according to previous theory[22] and as observed in some experiments[20]. With a phenomenological insertion of orbitally dependent quasiparticle weights $Z_{\alpha}$, with $\alpha$ an orbital index, Kreisel et al.[23] could fit inelastic neutron data on FeSe using very strong suppression of $d_{xy}$ weights and somewhat smaller suppression of $d_{xz}$ and $d_{yz}$ weights, assuming also a large ratio of at least 1.7 for the ratio $Z_{yz}/Z_{xz}$.

We therefore explore how orbital incoherence could improve the agreement of the susceptibility calculated from the current highly nematic electronic structure with experiment. In Fig. 4(b), we present the results for weak correlations, and the evolution of the susceptibility as correlations are enhanced by a substantial suppression of the $d_{xy}$ quasiparticle weight (see Appendix C for a brief discussion of the Ansatz of Kreisel et al.[23]). As anticipated from the contribution of the orbitally resolved susceptibility of the $d_{xy}$ orbital (see Fig. 4(a)), with reduction of $Z_{xy}$ comes the suppression of the Néel peak and concomitant moderate enhancement of the $(\pi,0)$ stripe peak. On the other hand the $(\pi,0)/(0,\pi)$ anisotropy is still much weaker than that reported in Ref. 5. In Fig. 4(c), we therefore show the effect of additionally increasing the $Z_{yz}/Z_{xz}$ anisotropy. It is easy to see that much larger $(\pi,0)/(0,\pi)$ anisotropies are obtained in this case, but also that substantially smaller quasiparticle weight ratios, of order $\sim 1.3$, are required compared to Ref. 23, due to the effect of inter-orbital nematic order introduced here.

On the other hand, these various hypothetical renormalizations do not lead to substantial changes in the gap structure obtained for FeSe within the corresponding spin fluctuation pairing theory. This is simply because in a multiband system even at $\omega=0$ contributions to $\chi(\bf q)$ arise from states tens or even hundreds of meV from the Fermi level[57]. By contrast, Fermi surface states determine the anisotropy of the effective pairing interaction completely. This effect can be seen easily in Fig. 5, where we plot the leading eigenvector of the linearized gap equation (see Appendix C) for the same three cases described in Fig. 4. The overall structure of the gaps and the density of states are seen to be virtually identical.

In this work we have pursued a scenario where nematicity originates entirely from NN Coulomb repulsion, even though it is well-known that onsite interactions alone may also drive a nematic instability as a precursor to stripe magnetism[4]. Interestingly, orbitally resolved studies within the spin-nematic onsite interaction-only scenario, do find sizable inter-orbital nematic susceptibilities[58]. Such studies, however, have only been performed in the tetragonal paramagnetic phase, and it remains to be seen whether spontaneous breaking in the inter-orbital channel, similar to the current proposal, can also arise from higher order processes solely from onsite interactions. However, the absence of magnetism in FeSe and the overall agreement of the here-presented results to the experimental facts, raises the question of whether indeed NN-repulsion is the generator of nematicity for FeSe. From cRPA[59] and DMFT[60] calculations it is known that interactions are generally larger in FeSe than any of the other iron-based superconducting systems, presumably because a lack of intervening spacer layers reduces the screening. Onsite interactions $U$, $J$ are known to lead to band renormalizations, in particular Fermi surface pocket shrinkage[47, 48] and orbital-dependent band narrowing. It is tempting to speculate that this, in turn, effectively enhances the importance of the unusually large $V$ in FeSe, thereby boosting instabilities driven by this channel.

Recently, another theoretical study investigated the possibility of nematic order lifting one of the electron pockets (above the Fermi level) in FeSe[43]. In agreement with the present work, a Lifshitz transition necessarily takes place as a function of temperature, and a superconducting gap structure consistent with experiments follows directly from the resulting FS with the missing $Y$-pocket. A main differences between Ref. 43 and the current approach is the nature of the nematic order; whereas Rhodes et al.[43] begin with a phenomenological $\bf k\cdot\bf p$ expansion around $\Gamma$, $X$ and $Y$, and explore the role of a large intra-orbital B_1g $d_{xy}$-nematic order parameter imposed by “hand” on the band structure, together with additional onsite Hartree shifts to this orbital, we have started from a nearest-neighbor Coulomb interaction, and shown that (self-consistently generated) inter-orbital $d_{xy}-d_{xz/yz}$ nematic components lift the $Y$-pocket. Thus, while a B_1g intra-orbital $d_{xy}$ nematic component is also present in our approach as seen from Fig. 2(d), the most important nematic components for generating a FS containing no electron $Y$-pocket are the distinct inter-orbital components.

As mentioned in the “Results” section, the nematic lowest-temperature phase advocated in this work, resides in the $\rm C_{2}$ (monoclinic) group, containing only a $C_{2}$ rotation around the $x$-axis as the remaining symmetry operation. This constitutes a clear prediction within the current scenario; the electronic sector should exhibit a double transition as the temperature is lowered. We are not aware of such evidence, e.g. from specific heat data[61, 62, 63], which may be because of the very small electronic entropy change at the second transition, or simply because the two transitions are accidentally close (in temperature) (see App. B for more details) in FeSe. In addition, if the nematic order couples strongly enough to the atomic lattice, such symmetry lowering could be tested by experiments sensitive to the overall crystal point group. In this regard, however, we note that recent experimental and theoretical studies have advocated for a more complex static crystal microstructure in FeSe than previously thought. Locally, FeSe appears to accommodate a myriad of inhomogeneous nanoscale lattice distortions, where only the spatially averaged structure complies with the standard tetragonal (orthorhombic) crystal symmetries at high (low) temperatures[64, 65, 44, 66]. This may in fact be the result of two nematic channels very close in in energy. Lastly, it is tempting to speculate that the inter-orbital nematic order discussed here, may be relevant for the recent experimental studies of structural transitions in Bi_1-xSr_xNi₂As₂, which exhibit transitions from a high-temperature tetragonal phase to a triclinic low-temperature phase[67].

Recent theoretical works explored the possibility of pinned local nematic order[68, 69, 70]. In particular, Ref. 70 explored the local disorder-induced nematicity from nonmagnetic impurities in the tetragonal phase at $T>T_{n}$. Several experimental works have reported evidence for such local nematic order above $T_{n}$[64, 65]. The results of Ref. 70 were obtained by applying a one-band model and relied on interactions leading to a single-component $\rm B_{\rm 1g}$ nematic order, and hence the question arises how local pinned nematic order gets affected by the presence of the substantial inter-orbital components found in this work? We have answered this question by performing a real-space calculation similar to that of Ref. 70, but generalized to the multi-orbital case. For the intra-orbital $\rm B_{1}$ order parameter $N_{\rm B_{1}}$ studied here, we observe that it enters the Landau free energy expansion in the exact same manner as the order parameter studied in Ref. 70. Therefore the local impurity-induced order exhibits a ”flower-shape” pinned by nonmagnetic impurities[70]. However, for the lower temperature phase, the spatially dependent local nematic impurity-induced order is distinct from that of Ref. 70, since the order parameter $N_{E}$ couples differently to the nonmagnetic impurity. It remains an interesting future study to investigate the detailed experimental consequences of these various local nematic orders.

In summary, we have discovered an inter-orbital nematic order generated from nearest-neighbor Coulomb repulsion for electronic models relevant for FeSe in particular, and perhaps the iron-based systems more broadly. A natural property of this kind of nematic order is the generation of highly anisotropic Fermi surfaces, featuring in some cases, only a single hole and electron pocket. We have shown how, for FeSe, this explains photoemission data and the experimentally extracted very anisotropic superconducting gap structure. However, a consistent picture of the neutron scattering data, cannot be straightforwardly obtained within this nematic scenario, without invoking additional self-energy effects, which in the simplest case, involves sizable corrections in the form of reduced orbitally dependent quasiparticle weights.

We thank M. H. Christensen and I. Eremin for insightful conversations. D.S. and B.M.A. acknowledge support from the Carlsberg Foundation. B. M. A. acknowledges support from the Independent Research Fund Denmark grant number 8021-00047B. P. J. H. was supported by the U.S. Department of Energy under Grant No. DE-FG02-05ER46236.