Non-local $d_{xy}$ nematicity and the missing electron pocket in FeSe

The origin of spontaneous electronic nematic ordering provides important information for understanding iron-based superconductors. Here, we analyze a scenario where the $d_{xy}$ orbital strongly contributes to nematic ordering in FeSe. We show that the addition of $d_{xy}$ nematicity to a pure $d_{xz}/d_{yz}$ order provides a natural explanation for the unusual Fermi surface and correctly reproduces the strongly anisotropic momentum dependence of the superconducting gap. We predict a Lifshitz transition of an electron pocket mediated by temperature and sulphur doping, whose signatures we discuss by analysing available experimental data. We present the variation of momentum dependence of the superconducting gap upon suppression of nematicity. Our quantitatively accurate model yields the transition from tetragonal to nematic FeSe and the FeSe$_{1-x}$S$_{x}$ series, and puts strong constraints on possible nematic mechanisms.


INTRODUCTION
The study of the iron-based superconductors ultimately revolves around understanding the interplay of superconductivity, magnetism and electronic nematicity. While the competition between magnetism and superconductivity has been extensively studied in other classes of superconductors, most notably the cuprates and heavy fermion compounds, the ironbased superconductors offer an ideal environment to study the interactions and consequences of the rotational symmetry breaking nematic state.
Currently, there is little consensus on the underlying mechanism that drives the formation of the electronic nematic state, with both charge-ordering and magnetic-ordering scenario's being proposed [1]. Nevertheless, it is often assumed that the leading consequence of the nematic mechanism on the ironbased superconductors is the lifting of the degeneracy of the d xz and d yz states of the Fe atoms. This has been thought to be responsible for many of the experimentally observed changes to certain physical properties [2]. This idea has additionally been held for the most well studied nematic iron-based superconductor, FeSe. This material exhibits an electronic nematic state at T s = 90 K as well as a superconducting state below T c = 8 K, without the formation of magnetic ordering [2][3][4]. Despite a relatively small tetragonal to orthorhombic lattice distortion to this material [5], particularly large anisotropic effects can be observed within measurements of the resistivity anisotropy [6], the spin susceptibility [7], the underlying electronic structure [8] and the momentum dependence of the superconducting gap [9][10][11][12][13][14].
However, it has become increasingly apparent that nematic order derived purely from a degeneracy breaking of the d xz and d yz states is unable to account for many of the experimental properties being observed within this material. Most notably, nematic ordering of the d xz and d yz orbitals predicts a Fermi surface topology which consists of one hole pocket and two electron pockets [15][16][17][18], whereas ARPES measurements of the nematic state of FeSe report the observation of one hole pocket and a single electron pocket [19][20][21][22][23][24].
This incorrect description of the Fermi surface has also made understanding the superconducting properties of FeSe challenging [10,11,25,26], without the addition of alternative anisotropic mechanisms, such as highly anisotropic orbital selective quasiparticle weights, [26,27] or orbital selective spin fluctuations [25].
Although nematicity is often thought of in terms of a lifting of the degeneracy of the d xz and d yz states, there is one additional form of orbital ordering, involving the relevant t 2g orbitals of the Fe atoms, which is consistent with the breaking of C 4 rotational symmetry. It is possible to have a non-local nematic ordering of the d xy states [28][29][30][31]. Although the d xy orbital does not conform to the B 1g irreducible representation of the point group associated with the tetragonal space group of many iron-based superconductors, an anisotropic hopping between the d xy orbitals of the nearest neighbour Fe-atoms, within a 2-Fe unit cell, does. This term is typically assumed to be negligible [28,31] and is often ignored within theories of the nematic state, however, there is no a-priori reason for this to be true. In fact, recent NMR [32] and ARPES [20,21] measurements have also suggested that the d xy orbital may be strongly affected by the onset of the nematic state.
With the goal to resolve the above mentioned inconsistencies between theory and experiment, we analyse in the following the question, what would be the consequence of nematicity, arising from the strongly correlated d xy orbital on the iron-based superconductors.
Focusing on the case of FeSe, we show from symmetry based arguments that if a sizeable non-local nematic ordering of the d xy orbital is the leading form of nematic instability, then many of the discrepancies between theoretical simulations and experimental measurements can be directly resolved. In particular, the missing electron pocket at the Fermi level and the highly anisotropic momentum dependence of the superconducting gap. We justify these claims by direct comparison of our theoretical simulations of the electronic structure with existing experimental measurements. As a crucial consequence of this d xy -dominated nematic scenario, we predict the existence of a Lifshitz transition of an electron pocket arXiv:2009.00507v2 [cond-mat.supr-con] 15 Mar 2021 Black dots denote Fe atoms. Se atoms located above (filled green) and below (empty green) the iron plane divide the system into sublattices A and B. b Fermi surface cut at kz = π for the tetragonal state, where the color describes the maximum orbital content of the band. The grey solid line describes the Brillouin zone boundary of the 2-Fe unit cell. Inset: Brillouin zones of the 2-Fe (solid grey) and 1-Fe (dashed black) unit cells corresponding to the "folded" and "unfolded" Brilloin zone, respectively. c Fermi surface map measured by ARPES of the tetragonal phase of FeSe at 100 K and d of the detwinned sample in the orthorhombic phase at 10 K taken from Ref. [19] under the Creative Commons Attribution 4.0 International License. as a function of temperature and sulphur doping, which can be verified experimentally.
This result not only provides a quantitatively accurate model to study the tetragonal to nematic phase transition in FeSe, but poses strong constraints on the possible mechanism associated with the nematic instability, and ultimately, on the mechanism of superconductivity in systems where nematic order and superconducting order appear in close proximity.

The B1g nematic state in FeSe
To understand the consequence of non-local d xy nematicity on the nematic state of FeSe we begin with a pedagogical discussion on the B 1g symmetric order parameters of the D 4h point group of the P4/nmm space group on the electronic structure of FeSe.
Within the P4/nmm crystal structure, the Fe atoms lie inplane on a square grid, and are connected by staggered outof-plane pnictide or chalcogenide atoms, which gives rise to a two-Fe atom unit cell, as shown in Fig. 1(a). This 2-Fe unit cell has primitive in-plane lattice vectors R 1 = R x − R y and R 2 = R x + R y which generate a Brillouin zone, defined by the coordinates k 1 = k x − k y , k 2 = k x + k y and k z , as shown in the inset of Fig. 1(b). In this 2-Fe Brillouin zone, the important high symmetry points are the k z = 0, Γ point and the k z = π, Z point at the zone center (k 1 , k 2 , k z ) = (0, 0, k z ) and the k z = 0, M point and k z = π A point at the zone corner (k 1 , k 2 , k z ) = (π, π, k z ). It is, however, useful to discuss the electronic structure of the P4/nmm lattice in terms of the lattice vectors R x and R y , and the momentum k x and k y , particularly to illustrate the effect of C 4 symmetry breaking. In this case, the Brillouin zone can be "unfolded" into an 45 degree rotated effective 1-Fe Brillouin zone, where the states around the M and A point are mapped to the X point at (k x , k y , k z ) = (π, 0, k z ) and the Y point at (k x , k y , k z ) = (0, π, k z ). Here, our model will retain the physical two-atom unit cell periodicity, however we will discuss symmetry breaking in terms of the coordinates k x and k y .
As a starting point, we employ the low-energy model of the d xz , d yz and d xy orbitals of the Fe atoms, following Ref. [33], and fit the parameters to quantitatively agree with ARPES data of the tetragonal state electronic structure of FeSe near the Γ (Z) point and M (A) point respectively [8,16,34], and Here the kinetic part H  0, this model describes a quasi-2D electronic structure which exhibits two corrugated d xz /d yz hole pockets at the Brillouin zone center with the inner hole pocket only present near k z = π, and two almost cylindrical d xz /d xy and d yz /d xy electron pockets at the Brillouin zone corners. The Fermi surface for this model is presented in Fig. 1(b). The states at the zone corner have majority d xz (d yz ) character along the length of the ellipse, but d xy at the outer tips. This model is in direct agreement with the low energy tetragonal electronic structure measured by ARPES, the Fermi surface of which, from Ref. [19], is presented in Fig. 1(c). In Fig. 1(d) we also show the experimentally measured Fermi surface of the nematic state of FeSe, which the nematic order parameter must ultimately reproduce.
If we focus on the d xz , d yz and d xy orbitals of the Fe atoms, within a single two atom unit cell, then we can write down four possible C 4 symmetry breaking nematic orders in the B 1g -channel, which arise either from on-site interactions on the Fe-atoms or from interactions between two neighbouring Fe-atoms on different sublattices. We list these four order nematic process real space momentum space expansion near expansion near xy yz xz TABLE I. All B 1g nematic order parameters for FeSe. All possible nematic order parameters within the B1g irreducible representation of the P4/nmm space group involving interactions between dxz, dyz and dxy orbitals within a single 2-Fe unit cell. Results presented in real and momentum space as well as Taylor-expanded in the vicinity of the high symmetry points (following the notations of Ref. [28]). Here we use coordinates of the unfolded Brilloiun zone. A sketch of of each nematic process is depicted in the first column where solid (dashed) arrows are associated with a positive (negative) amplitudes of the corresponding nematic order parameter. Note that dxy nematicity is a purely non-local phenomenon. parameters in Table 1 along with a schematic cartoon of the nematic process, mathematical formulation in real and momentum space and the effect each term has on the hole and electron pockets.
The first term, Φ 1 , describes a degeneracy breaking of the on-site d xz and d yz states. This form, often referred to as ferro-orbital ordering [36], generates an equivalent splitting of the d xz and d yz bands which affects both the hole and electron pockets equally, as shown in Fig. 2(a). This type of nematic order was originally discussed in the literature [15,[37][38][39], however, it later became apparent that whilst the evolution of the hole pocket is correctly captured by Φ 1 , the evolution of the electron pocket into a d xz dominated peanut, elongated along the k y axis, and a larger outer d xy dominated ellipse was not consistent with experiment [16,37,[40][41][42][43]. It was realised that a sign change in the nematic order parameter between the hole and electron pocket, involving the d xz and d yz orbitals was required [16,41,42]. Interestingly, this sign change of the nematic order between electron and hole pockets is analogous to the sign-changing s ± -wave superconducting gap and indicates that the same electronic interactions are most likely in play [29,30].
The second and third terms, Φ 2 and Φ 3 , describe non-local nematic ordering between nearest neighbour atoms involving the d xz and d yz states. Φ 2 describes an anisotropic hopping potential such that hopping is favoured in the x direction, but disfavoured in the y direction for both orbitals. This has been referred to as d-wave nematic order within the literature [17,18,41,44,45]. Conversely, Φ 3 describes a global increase in the hopping of the d xz orbital, and a global decrease in the hopping of the d yz orbital, often referred to as extended swave nematic order [16][17][18]. As shown in Table 1, Φ 2 only affects the electron pockets at the X/Y point, and Φ 3 only affects the hole pocket at the Γ point. We present the Fermi surface generated by the inclusion of Φ 2 or Φ 3 in Fig. 2(b) and 2(c) respectively.
The fourth term, Φ 4 is a purely non-local nematic ordering between the d xy states of nearest neighbour atoms, analogous to that of Φ 2 . This term, has been mentioned as a previously possible form of nematic order [28,29], however it was assumed to only produce a minor effect on the electronic band dispersion [31,44] and is often neglected entirely. If however, Φ 4 is made sufficiently large (> 50 meV) then it would be possible to generate a Lifshitz transition, which can induce a one-electron pocket Fermi surface. We show the generated Fermi surface obtained from the inclusion of Φ 4 > 50 meV in Fig. 2(d). As mentioned above, these are the only four symmetry allowed order parameters that can describe the evolution of the nematic state within the P4/nmm space group. Therefore, some combination of these four order parameters should be able to reproduce the Fermi surface of the nematic state of FeSe as measured by ARPES. As shown in Table 1, Φ 1 contributes equally to both the states around the Γ point and the X/Y point. However, Φ 2 and Φ 3 contribute solely to the X/Y and Γ point, respectively. It is therefore not possible to fit the magnitude of Φ 1 , Φ 2 and Φ 3 directly to experiment. This is because experimental measurements can only extract the total change in the band dispersion, which for the hole pocket is defined as Φ h xz,yz = Φ 1 + Φ 3 and for the electron pocket is Φ e xz,yz = Φ 1 + Φ 2 . In this model we therefore write the nematic order parameters as for the hole pocket, and for the electron pocket. We then fit the parameters Φ h xz,yz , Φ e xz,yz and Φ 4 directly to experimental values [19][20][21]. In Fig 3(a-c), we present the Fermi surface and band dispersion around the Γ and M point for the tetragonal state, i.e Φ h xz,yz =Φ e xz,yz =Φ 4 =0. In Fig, 3(d-f) we next present the Fermi surface and band dispersion's using Φ h xz,yz =15 meV, Φ e xz,yz = -26 meV and keep Φ 4 = 0 meV. This purely d xz /d yz nematic order correctly describes the 37 meV band separation of the two hole bands (which is obtained by a combination of the nematic splitting and spin orbit splitting where the spin orbit coupling parameter, λ h SOC = 23 meV [35]) and correctly describes the shape and orbital character of the d xz dominated elliptical hole pocket. Additionally, the experimentally observed peanut shaped electron pocket with d yz orbital content along the length of the peanut, and d xy orbital content at the tips, is also quantitatively captured. However, this scenario, with Φ 4 = 0 meV, also describes a second large electron pocket of predominantly d xy orbital character, which is not detected in experimental measurements. This type of order parameter is the main form considered within the literature [3,8,46]. In Fig. 3(g-i) we now present our scenario, where we use the same values of Φ h xz,yz and Φ e xz,yz as in Fig. 3(d-f) but now we set Φ 4 = 45 meV. We note that the Φ 4 term, as it is discussed in Table I, lifts the degeneracy of the d xy saddle point located at -50 meV at the M (A) point in the tetragonal state. This term alone however would ensure that the splitting of the saddle point is symmetric, such that the lower branch of the d xy band would also be pushed to a lower binding energy. Since a large shift in the lower d xy band is not observed in the experimental data [8] we add a finite energy shift to the on-site d xy energy in the nematic state which then generates an asymmetric nematic splitting. This term has the consequence that the lower branch of the d xy saddle point, remains in place, whilst the upper branch increases in energy by Φ 4 + ∆ xy where a value of ∆ xy = 40 meV was found to best fit the experimental data. The result of this, is that the upper d xy band is placed 36 meV above the Fermi level as shown in Fig. 3(i). Note that in general the onsite energy xy can change when entering the nematic state which corresponds to a Hartree shift of the d xy states. While we do not explicitly calculate this Hartree shifts of the orbitals selfconsistently, we do calculate the renormalization of the total chemical potential when entering into the nematic phase and show that it is indeed strongly affected by the transition. We explicitly assume that this Hartree shift has the same meanfield-like temperature dependence as the nematic order. It is important to note that this shift of the onsite energy only affects the precise size of the peanut shaped electron pocket and is not needed to obtain the correct Fermi surface topology of FeSe, which may be obtained if we set Φ 4 larger than 50 meV. It is however important to yield the correct position of the d xy dominated bands far below the Fermi energy, which is our reason for including it here. Note, as indicated by the vertical dashed line in Fig. 3(c,f,i), we find that the Fermi momentum (k F ) associated with d xy states of the peanut tips to be smaller in the tetragonal phase (c) than in the orthorhombic phase for our parametrization in (i) which is in good agreement with Ref. [19,20]. This is not the case if one neglects the d xy nematic order as seen in (f). Finally, we point out that the exact values for Φ 4 and ∆ xy are not particularly well constrained as ARPES data can not detect the actual position of the upper d xy band. However, as we will discuss below, a value of roughly this magnitude is required to correctly describe experimental properties as a function of temperature and sulphur doping.  The importance of a sizeable nematic order between d xy orbitals in the correct description of the electronic structure of FeSe is the main feature of our model of nematicty in FeSe. In contrast to the previous proposals [31,46,47] this scenario does not require strong incoherence of the d xy -dominated electron pocket. Furthermore, as we also discuss below, this nematic order between d xy orbitals is also required to resolve the multiple discrepancies between theory and experiment discussed in the literature [11,25,26,47,48].

Experimental evidence for non-local dxy nematicity
The main consequence of a sizeable non-local d xy nematic term is the generation of a Lifshitz transition at the corner of the Brillouin zone shortly after the onset of the nematic state. This Lifshitz transition is directly responsible for the evolution between a two-electron pocket Fermi surface in the tetragonal state into a one-electron pocket Fermi surface in the nematic state. Here we discuss available experimental data which can be interpreted in terms of this Lifshitz transition driven by nematic order and compare these experiments to our model calculations. In Fig. 4 we now compare our simulated electronic structure of the A point of FeSe, using the parameters discussed in Fig.  3(g-i), with the experimental ARPES data of detwinned crys-tals of FeSe from Ref. [19]. In Fig. 4(d) we present a close up of the Fermi surface predicted around the A point of our model which consists of a single peanut shaped pocket with dominant d yz orbital character along the length of the peanut and dominant d xy orbital character at the tips. The shape and size of this pocket is in direct agreement with the experimental measurement shown in Fig. 4(a) as well as the data published in Ref. [20,21,24]. In Fig. 4 (e,f) we present a cut of the band dispersion centred at the A point along the k x and k y axis, respectively. Our model exhibits a d yz saddle point located at -5 meV below E F , which can be observed along the k x axis as a hole like dispersion that almost reaches E F and along the k y axis as a very shallow electron like dispersion. A second saddle point of d xy orbital content is also present at -55 meV, which can equivalently be observed as an electron like dispersion along k x and a hole like dispersion along k y . Finally, a hole-like d xz band is also located around -60 meV along both the k x and k y axis. This simulated band dispersion is entirely consistent with the ARPES data presented in Fig. 4(b,c), as well as the data published on detwinned crystals of FeSe in Ref. [20,21,24]. However in the experimental data the d xy and d xz bands at ∼-50 meV can not be separately resolved due to self energy broadening effects [49,50]. This electronic structure is also fully consistent with ARPES measurements on "twinned" crystals, [11,12,16,39,41,42,51].
Despite the removal of an entire electron pocket within Comparison with ARPES data. a Experimental Fermi surface of the A point of FeSe measured by ARPES on a detwinned crystal taken from Ref. [19]. b,c Experimental band dispersion along the kx and ky axis respectively centered around the A point taken from Ref. [19]. d-f Corresponding simulated Fermi surface topology and band dispersions of the orbital projective model with the nematic order parameter.
g Change in the chemical potential as a function of temperature using a mean field temperature dependence Φi(T ) = Φi(0) 1 − ( T Ts ). h Corresponding change to the band positions at the A point as a function of temperature. i Experimental temperature dependence of the spectral function taken from Ref. [16]. j Simulated spectral function as a function of temperature calculated using Eq. (6) with a broadening parameter of Γ ∼ 10 meV. the nematic state, we find that the Luttinger theorem is conserved within this model. In Fig. 4(g) we plot the change in chemical potential µ(T ) self consistently calculated at a fixed number of particles and using a mean field temperature dependence of the nematic order parameter of the form , see Ref. [35] for details. We find that, within our parametrization, a shift of 13 meV is required to ensure particle number conservation between the tetragonal and nematic state, which is in good agreement with the 10 meV shift extracted from ARPES measurements [39,52]. Note that our parametrization accounts for charge neutrality, which is conserved throughout the whole doping range by the renormalization of the chemical potential µ. We emphasize that charge neutrality is also conserved at the Lifshitz transition of the d xy dominated electron pocket. Additionally, in Fig. 4(h) we present the evolution of the band positions at the A point calculated using the same temperature dependence and chemical potential shift. Under this assumption, we see that the Lifshitz transition occurs at T LT ∼70 K, shortly after the nematic transition at T s = 90 K.
Finally, in Fig. 4(j) we plot the temperature dependence of the spectral function, Here we have used an energy broadening parameter of Γ = 10 meV for all bands to simulate self energy broadening. We see that between 0 < T < 70 K the d xy band below the Fermi level resides very close in energy to the d xz band. Consequently, in the spectral function these two bands merge to a single broad peak at low temperatures, which together with the d yz excitation near E F suggests that two peaks would be observed by the spectral function, which is in good agreement with the experimental measurements reported in Ref. [16,20,41,42]. We show the experimental measurement of Ref. [16], in Fig. 4(i) for comparison.
We note that there has been some ambiguity as to the experimental interpretation of this spectral function, particularly regarding the orbital content of the lower peak at -50 meV. With some authors claiming that the lack of a clearly observable splitting of the peak closest to E F in the spectral function at 100 K indicates that the d xz and d yz states must remain degenerate [11,16,42,52]. However, in our simulation, we find that the direct observation of this splitting would be challenging, particularly with large self energy broadening effects that are present at the A point in the tetragonal state [16,41,42,50]. Coupled with the observation of d xz spectral weight of the band at -50 meV [20,21,41], as well as the symmetry allowed possibility of an orbital transmutation of d xz and d xy weight as a function of temperature [31], we find that our interpretation is entirely consistent with the measured ARPES data.
Further evidence of a Lifshitz transition of the d xy orbital can be found in Ref. [20]. There, Yi. et. al. reported temperature dependent ARPES measurements of the A point of detwinned crystals of FeSe, where the orthorhombic domains were aligned along a given axis by uniaxial strain. When studying the temperature dependence of the spectral function around the A point at +10 meV above the Fermi level, they report the observation of a shrinking inner electron band which disappears at this particular energy at approximately T = 60 K, as shown in 5(a). Our calculation of the momenta at which the electronic dispersion crosses the energy of +10 meV above the Fermi level is presented in Fig. 5(b). Indeed, it yields the same result as shown in Fig. 5(a) with excellent agreement on temperature dependence and band position. This agreement between theory and experiment can be interpreted as a smoking gun for the d xy Lifshitz transition, however, we caution that there are some experimental technicalities, which can complicate the interpretation. For example, changing the temperature of a uniaxially strained system, as was performed in Ref. [20], may alter the relative ratio of the orthorhombic do-main orientations, due to additional thermal contraction. This could imply that the disappearing intensity may instead arise from a change in population of orthorhombic domains, rather than a band shifting above E F . Nevertheless, if confirmed by other experimental measurements, this could be the most direct evidence for a d xy Lifshitz transition. A similar observation has also been made for the sulphur doped FeSe 0.9 S 0.1 system [22,23]. Alternate evidence for a Lifshitz transition at the zone corner, which does not require uniaxial strain, can be found by studying Quantum Oscillation measurements of the sulphur doped, FeSe 1−x S x , systems. Isovalent sulphur doping is known to suppress the nematic state, which, despite slightly changing the anion height of the lattice, leaves the underlying electronic structure effectively equivalent to that of the tetragonal state of FeSe [8,54,55]. Thus, by varying the sulfur content in FeSe 1−x S x it is possible to indirectly study the temperature dependence of the nematic order. Recently, Shubnikovde Haas oscillation measurements of the sulphur doped series have been reported by Coldea et. al. [53]. These experiments detect extremal areas of Fermi surface pockets, which in FeSe corresponds to the 2D area of the hole and electron pockets at both k z = 0 and k z = π respectively.
Previous interpretation of this experimental data was based on the assumption that two electron pockets are present in the nematic state. However, we find that our scenario, where there is only one electron pocket, is also qualitatively consistent with the quantum oscillation data if one assigns the Z and A pocket (labelled as 1 and 2 in Fig. 6 (a,c)) to the δ and γ frequencies of Ref. [53] and the M and Γ pockets (labelled as 3,4 in Fig. 6 (e,g)) to the α and β frequencies. Our scenario does not predict an additional frequency, which is however not detected in the experiment for the entire doping range.
As a function of sulphur doping, a Lifshitz transition associated with the d xy band should lead to a sharp sudden change in extremal areas measured by this experimental technique. We illustrate this in Fig. 6(i). By again assuming a mean field dependence of the nematic order parameter , where x 0 = 18%) we find that the four unique extremal areas present in the nematic state remain relatively consistent between 0% and 13% sulphur doping, with the inner hole pocket of the Z point crossing the Fermi level at 10%. However, at the Lifshitz transition of the d xy band, we find that all extremal areas sharply increase in magnitude and two additional extremal areas are generated. Bringing the total number of unique extremal areas in the tetragonal state to seven. The reason for this sharp increase is twofold, the first is the sudden change of chemical potential as discussed in Fig. 4(g), which mainly affects the hole pockets, and the second is the sudden change of the electron pockets from a two-fold symmetric, peanut-like, dispersion to a four-fold symmetric, petal-like, dispersion. Whilst the two largest areas corresponding to the k z = π hole pocket and the outer electron pocket are well separated from the k z = 0 pockets, the other five extremal areas are likely to have very Corresponding experimental data taken from Ref. [53]. We have included a dashed circle in i,j to highlight the evidence for a re-entering dxy band. Note that we use a different color code in a-i than the authors of Ref. [53] and that we assign our pockets differently to the data and trend lines in j.

similar (and very small) values.
Although a Lifshitz transition of the electron band was not the interpretation of Ref. [53], we observe that the data points measured for x > 15% are in fact consistent with a re-entering d xy pocket at the A and M point, we display the experimental data of Ref. [53] in Fig. 6(j) for comparison. Very similar results can also be observed by suppressing nematic order under application of physical pressure to FeSe [56]. We note that the conversion of these simulated areas into units of Tesla using the Onsanger relation (F = A 2πe where A is the area) indicate that the simulated areas are approximately 70-80% of that measured by experiment for all doping values, nevertheless the qualitative trends are consistent.

Momentum Dependence of the superconducting gap
We now consider how the superconducting state will be affected by the presence of non-local d xy nematicity, and in particular how the gap structure will evolve from the nematic to tetragonal state.
Superconductivity in the iron-based superconductors is believed to be driven by the interband repulsive interactions enhanced by spin fluctuations, and indeed the numerous comparisons between theory and experiment have shown that this assumption can correctly capture the symmetry of the superconducting gap [57][58][59]. The gap symmetry of the nematic state of FeSe has been a subject of intense experimental investigation, with multiple ARPES and scanning tunnelling microscopy (STM) measurements reporting a highly anisotropic momentum dependence of the superconducting gap [9][10][11][12][13][14]. These measurements revealed regions in momentum space which are either nodal or posses a very small superconducting gap. Previously, it was shown that the solution of the linearised superconducting gap equation assuming a one-electron pocket Fermi surface and random phase approximation (RPA) correctly reproduced the experimental data in the nematic state [11], however this was achieved by phenomenologically removing the second electron pocket that was still present in the underlying band structure. In the scenario we have considered here, with sizeable d xy nematicity, the second electron pocket is shifted well above the Fermi level and we are now in a position to study how the momentum dependence of the superconducting gap evolves as nematicity is suppressed, and a d xy Lifshitz transition occurs, which can be experimentally probed by measuring the momentum dependence of the superconducting gap as a function of sulphur doping.
We utilize our orbital projected band model which allows for an analytical treatment of the pairing problem by solving mean-field BCS gap equations self-consistently at zero temperature. Within the projected band model we adopt the interaction Hamiltonian of Ref. [60] which combines Kohn-Luttinger and spin-fluctuation approaches. This leads to a projection on those pairing channels which contain intra-orbital U he and inter-orbital J he and J ee pairings that results in interband Cooper-pair hopping between hole and electron pockets, corresponding to (π, 0), (0, π) spin fluctuations, and between both electron pockets, corresponding to (π, π) fluctuations, respectively. The interaction Hamiltonian then reads where µ, ν ∈ {xz, yz}. Even though our conclusions in the  previous section suggest a sizeable contribution of d xy orbitals to nematicity, for superconducting pairing the d xy orbital is not so important. Superconductivity is dominated by interband, and intra-orbital, Cooper pair scattering. As there is no d xy weight on the hole pocket, the d xy orbital can only contribute to the superconducting pairing via scattering between two electron pockets (J ee ), which is considerably weakened if only one electron pocket is present. In what follows we neglect the contribution of the d xy orbital to superconducting pairing. Later, in Sec. D, we will confirm this assumption by solving the full linerised gap equation for a 10-orbital tight binding model, and including the d xy orbital in the pairing.
We define superconductivity in orbital space, perform a mean-field decoupling and solve the BCS-gap equations together with the particle number equation self-consistently for the constant orbital pairing gaps assuming U he J he = J ee . In the gap equations we focus on the large hole pocket, the peanut shaped electron pocket and the incipient electron band above the Fermi level. We further neglect contributions from inter-band Cooper-pairing and integrate within an energy cutoff of 50 meV around the Fermi level. Further details are given in Ref. [35]. The superconducting gap on the outer hole pocket is given by with the constant orbital gaps ∆ h xz/yz and the k-dependent orbital dressing factors |a h xz/yz (k)| 2 = | d xz/yz |h out | 2 , which correspond to d xz/yz orbital content of the hole pockets, respectively. Assuming small SOC at the A point the band-gap on peanut shaped pocket is given by in good approximation by ∆ e out (k) ≈ |a e yz (k)| 2 ∆ e yz + |a e xy (k)| 2 ∆ e xy ≈ |a e yz (k)| 2 ∆ e yz (9) , with |a e yz/xy (k)| 2 = | d yz/xy |e X | 2 the d yz/xy orbital weight on the peanut [35].
In Fig. 7(a) we show the orbital weight |a h xz/yz (k)| 2 on the hole pocket as function of the Fermi angle measured with respect to k x direction. In experiment, the gap at the hole pocket exhibits deep minima (and possibly accidental nodes) at θ(k F ) = ±π/2 and has a maximum at θ(k F ) = 0 resembling the d yz orbital weight. According to Eq. (8) this behaviour requires ∆ h xz ∆ h yz . In Fig. 7(b) we present the solution of the BCS gap equations as function of U he , which yield the large ratio ∆ h yz /∆ h xz necessary to reproduce the correct angular dependence of ∆ h out (k). A similar result for the hole pocket was also obtained in Ref. [44] where nematicity was solely given by Φ e xz,yz . Within their model, however, two electron pockets were present at the Fermi level and they were not able to correctly reproduce the angular dependence of the gap at the peanut shaped electron pocket. Ref. [25] also calculated Eq. (8) and Eq. (9) for a model where nematicity was given solely by Φ e xz,yz . However here they added an additional strong anisotropy to the superconducting pairing such that the second electron pocket did not contribute to superconductivity, but was still present at the Fermi level, in disagreement with ARPES measurements as discussed in Section B.
The weight |a e yz/xy (k)| 2 of the d xz and d xy orbitals on the peanut shaped electron pocket is shown in Fig. 7(b). Experimental results suggest that the angular dependence of the gap follows the d yz orbital weight of the peanut which has a minimum at φ(k F ) = 0 and is almost isotropic elsewhere. This behaviour is captured by Eq. (9).
Furthermore, in Fig. 7(c) we show that |∆ h yz | > |∆ e yz |, which leads to max(|∆ h out (k F )|) > max(|∆ e out (k F )|) as it is seen in experiment. In Fig. 7(e) and (f) we show the order parameters ∆ h out (k F ) and ∆ e out (k F ) projected onto the Fermi surface and as a function of the Fermi angle and find excellent agreement with experimental data [9][10][11][12][13][14], of which we show data extracted from STM measurements of Ref. [10], depicted as red and blue circles in Fig 7(e). For completeness, in Fig. 7(d) we present the total density of states (DOS) ρ(ω) = − 1 π k Tr τ0+τ3 2 ImĜ(ω, k) calculated for the gap values that correspond to Fig. 7(d). The DOS is nearly Vshaped, dedicated to the almost nodal gaps, and exhibits two peaks at energies ∼ 1.5 meV and ∼ 2.2 meV, which correspond to the gap maximum of the electron and hole gap, respectively.
Finally, in Fig. 8 we make a prediction about how the Fermi surface gap structure in the orthorhombic state of FeSe 1−x S x will evolve upon decreasing the nematic order parameter, for example as a function of increasing sulphur content. A schematic phase diagram highlighting the doping induced Lifshitz transition at the zone corner is shown in Fig. 8(d). Here, we model the doping dependence of the nematic order param- , where x 0 = 18%, and self consistently calculate the pairing gaps as well as the chemical potential. Note, however, that at this stage we keep the interactions momentum-independent and do not alter its value as the doping increases and thus do not account for altering nesting conditions which will affect the pairing strength and gap size but not the momentum dependence. We observe that throughout the entire doping range the overall pairing symmetry exhibits a sign-change between the hole and electron pockets.
At the same time, we find that within a scenario of the d xydominant non-local nematicity there is a clear impact of the d xy Lifshitz transition on the angular dependence of the gap on the electron pockets as can been seen in Fig. 8 (b,c). In particular, additional maxima/minima or accidental nodes for x = 0.17 of the gap on the inner and outer electron pockets at θ = ±π/2 develop near the Lifshitz transition presumably due to enhanced (π, π) interaction between them. This could be an explanation for the two distinct superconducting gaps detected as a function of sulphur doping [61,62]. This specific angular dependence may be verified by ARPES/STM measurements and is one of the predictions of our d xy -nematic scenario. Observe also that the angular dependence of the superconducting gaps in the tetragonal phase, x > 0.18, is probably more complex due to recent observation of the timereversal symmetry breaking [4].

Comparison with 10-orbital tight binding model
Throughout this manuscript we have used a minimal orbital projective model to describe the electronic dispersion around the Fermi level. To confirm the generality of our conclusions and to correctly account for d xy contribution to superconductivity we have also studied a 10-orbital tight binding model that has additionally been fit to ARPES data in the tetragonal state [34]. We then use the same form of nematic order parameter which for the hole pockets can be written as and for the electron pockets can be written as where d † xy A/B (d xy A/B ) is the creation (annihilation) operator for the d xy orbital on sublattice A or B in the two-atom unit cell respectively. Here we use the values Φ h xz,yz =7 meV, Φ e xz,yz = -14 meV and keep Φ 4 = 28 meV.
To account for the Hartree term introduced in Eq. (5) we additionally include the two equations This term ensures the asymmetric splitting of the two d xy bands around the M (A) point are captured whilst leaving the additional d xy band, located around -50 meV at the Γ (Z) point, unaffected. As with the case of the orbital projective model, Eq. (12) does not break any additional symmetries, or change the Fermi surface topology, however it is required to reproduce the fine details of the band dispersions below the Fermi level. We set ∆ xy = 28 meV within this formalism. We present the Fermi surface and band dispersion for this model in Ref. [35] showing nearly identical band dispersions to the orbital projective model of Fig. 1. We now study the momentum dependence of the superconducting gap predicted from this model. To do this we solve the standard linearized superconducting gap equation [63], Here, we integrate the pairing vertex, Γ µν (k, k ) over all k states on the Fermi surface for each band, ν. The eigenvectors corresponding to the largest eigenvalue then describe the symmetry of the leading superconducting instability suggested by the pairing vertex. Importantly, within this framework we treat all orbitals, including the d xy orbital, as equally coherent [64].
We use a pairing vertex that describes spin-fluctuation mediated superconductivity utilising the random phase approximation (RPA) in the presence of spin orbit coupling [11,65]. The exact form of the pairing vertex is rather long and left to the appendix, however the important properties that control this pairing vertex are i) the orbital content present at the Fermi level and ii) the form of the spin susceptibility used to generate pairing.
One has to bear in mind, however, that the RPA approximation is only valid in the regime where ωsf EF < 1, where ω sf is the spin fluctuation frequency and E F is the Fermi energy [66], which is not true for FeSe where the Fermi energy is <10 meV [39]. However, the exact details of the spin susceptibility, and pairing vertex, are less important to the superconducting properties of FeSe than the available states present around the Fermi level. To emphasise this, we solve the linearised superconducting gap equation using two distinct forms of the spin susceptibility. In the first we employ the RPA pairing vertex, calculated explicitly from the itinerant electronic structure in the tetragonal and nematic state, where, a s µ (k) describe the eigenvector of orbital s and band µ at momentum k, E µ (k) is the corresponding eigenvalue and f (E) is the Fermi function calculated at T=10 K in the nematic state and 100 K in the tetragonal state.
In the second form, we employ a greatly simplified phenomenological form of the spin susceptibility, such as has been used to study the cuprates and other ironbased superconductors [67,68]. Here, Q i are the antiferromagnetic wavevectors (π, 0) and (0, π) and X Qi is a constant proportional to the real part of the spin response at the corresponding Q i , with units eV −1 , which we set to 1 for each Q i for illustrative purposes. We additionally set ξ (a measure of the correlation length of the spin fluctuations) to unity and neglect inter-orbital contributions in this approximation. The spin susceptibility calculated using Eq. (14) and Eq. (15) are presented in Fig. 9(a) and (d) respectively. We observe that the itinerant form of the spin susceptibility predicts maximum spin fluctuations to occur at (π, π), in agreement with previous models of the tetragonal state [15,27,45,47], and that the nematic order parameter only induces a slight increase (decrease) to the spin susceptibility at (π, 0) ((0, π)). Although this form of spin susceptibility does not correctly describe the fully dynamics of the spin response as seen in inelastic neutron scattering experiments [7,69] we find that remarkably, the correct highly anisotropic momentum dependence of the superconducting gap is still obtained as the leading solution to the linearised gap equation of the nematic state, as shown in Fig. 9(b). This is equivalently the case for the phenomenological spin susceptibility in Fig. 9(e). The reason for this agreement is that, due to the Fermi surface topology, superconducting pairing in the nematic state is controlled by approximately only two scattering vectors Q ≈ (0, 0) and Q ≈ (π, 0), such that the rest of the spin susceptibility structure does not impact the superconducting properties.
However, in the nominally tetragonal state, with two electron pockets, we find that the form of spin susceptibility does have an effect on the leading pairing symmetry. The itinerant pairing vertex, with maximum spin fluctuations at (π, π), suggests d-wave superconductivity will be the leading instability, as shown in Fig. 9(c). However, if (0, π) and (π, 0) fluctuations are assumed to be dominant then s ± becomes the leading instability ( Fig. 9(f)).
Additionally, we find that the only contribution of the d xy orbital in superconducting pairing is the possible generation of accidental nodes at the tips of the electron pockets, which verifies that the assumption made in Sec. C is valid. This observation emphasizes the robustness of our prediction for the impact of the Lifshitz transition of one electron band due to d xy −nematicity on the superconducting gap near the nematic to tetragonal transition.  9. Comparison of the superconducting gap simulations with a 10-orbital tight binding model. a Spin susceptibility calculated using the itinerant spin susceptibility defined in Eq. (14). b,c Eigenvectors of the leading eigenvalue from Eq. (13) calculated using the itinerant spin susceptibility in b the nematic state and c The tetragonal state. The inset defines the angular dependence of the gap for each pocket in arbitrary units, taken anticlockwise from the kx-axis. ex and ey define the momentum dependence of the gap at (π, 0) and (0, π) respectively, whilst ho and hi describe the inner and outer hole pocket. d-f Equivalent results obtained from the phenomenological spin susceptibility defined in Eq. (15). Here we set U=0.5 eV, U = U 6 ,J = 0.1U and J = J throughout and use the unfolded 1-Fe unit cell [34]. Further details can be found in the supplemental material.

DISCUSSION
The d xy orbital has been largely underestimated within theories of the nematic state for the iron-based superconductors. With most theories focusing solely on the degeneracy lifting of the d xz and d yz orbitals. Although, the possibility of anisotropic hopping between the d xy orbitals of nearest neighbour atoms have been also considered [29,31,70], it was thought to have a minor role for the electronic structure modification in the nematic state. The resulting contradiction of theoretically predicted Fermi surface topology in the nematic state with that experimentally measured by ARPES had previously been attributed to the incoherence of the d xydominated electron pocket. This was indirectly supported by DMFT calculations, arguing about the stronger correla-tions effect on the d xy orbital [71], yet we are not aware of DMFT calculations that would correctly reproduce the measured one-electron pocket Fermi surfaces of the nematic state of FeSe. In order to generate the experimentally observed Fermi surface, we adopted another scenario where the sizeable d xy nematic order parameter is larger than the d xz /d yz term and find that it allows consistent explanation of several existing experimental controversies. In addition, this scenario predicts a Lifshitz transition of one electron band in the sulphur doped FeSe 1−x S x systems, which should occur prior to the orthorhombic to tetragonal transition. This may explain recent resistivity measurements in high magnetic field which have reported quantum-critical behaviour at a sulphur doping value of x = 0.16 [72], two percentage points lower than the disappearance of the nematic state which occurs at x = 0.18 [2], as well as the sudden change in superconducting gap observed in specific heat and STM measurements [61,62]. Assuming that the nominal doping values are correct between different experiments, we expect that the observed quantum critical behaviour is associated with the Lifshitz transition predicted here.
Apart from the interest on its own, the sizeable d xynematicity poses further constraints for microscopic theories of nematicity. The anisotropic hopping pathway of the d xy orbital, coupled with the experimental observation of a momentum dependent d xz /d yz rearrangement suggest that nematicity is a strongly non-local phenomena (likely of magnetic rather than purely orbital origin) and may be governed by nearest neighbour interactions [73]. Observe also that this order needs to be re-examined within the RG weak-coupling analysis performed previously [29,30,74]. It is also interesting to note that local ab-initio calculations, such as LDA, LDA + DMFT, and even QSGW+DMFT are unable to account for such a prominent evolution of the bands upon entering the nematic state [42,50,75]. Although a recent DFT calculation proposed an E u symmetry nematic order parameter, using symmetry preconditioned wavefunctions [76], that generates an electronic structure where one of the electron pockets is above the Fermi level. The lack of a consistent ab-initio description of the nematic state, coupled with our observation of a maximum nematic order parameter on the most correlated d xy band strongly suggests that nematicity arises from a nonlocal strongly-correlated phenomenon. It is also interesting to note that Te doping induces an orbital selective Mott transition of the d xy orbital, and suppresses the nematicity [77,78].
The scenario of the d xy dominant nematicity has also important consequences for alternative theoretical explanations for experimental properties of FeSe, such as the need for highly anisotropic orbital-selective quasiparticle weights or spin fluctuations to describe the momentum dependence of the superconducting gap [10,[25][26][27]79]. These theories were born from the observation that a two-electron pocket Fermi surface could not correctly describe experimental measurements, and that some form of anisotropy was required, either within the quasiparticle weight of the d xz and d yz orbitals, or in the direction of the spin fluctuations. Importantly, both ideas have the consequence that they suppress any contribution of a second electron pocket from superconducting pairing. What we have shown here is that the discrepancies between theory and experiment can be simply explained as a consequence of starting with an incorrect Fermi surface topology. This point has been also argued in several publications [11,48], however, in these papers a justification for the removal of the second electron pocket was lacking. Here, we have presented a scenario which can justify the removal of the second electron pocket, by assuming large nematic ordering in the d xy orbital.
We want to point out that the addition of a Hartree shift alongside the nematic order parameters is necessary in our model to reproduce the asymmetric nematic splitting of the d xy band, suggested by experiments [20], which constrains the position of the d xy dominated bands far below the Fermi energy. Future work is required to elucidate the potential link of this term with nematicity.
To conclude, we have studied a model of FeSe and FeSe 1−x S x which exhibits dominant nematic order in the d xy channel and have shown that this assumption can naturally explain a multitude of unusual experimental properties observed within this system. This finding greatly constrains the possible mechanisms that can be attributed to the phenomena of nematicity.
Note added: During the submission process of our manuscript we became aware of a theoretical paper by Steffensen et. al. [80] also addressing the question of a single electron pocket Fermi surface of FeSe.
Acknowledgements We thank Andrey Chubukov, Yuji Matsuda, Taka Shibauchi, Peter Wahl, Amalia Coldea, Steffen Boetzel, Matthew Watson and Timur Kim for insightful conversations. LCR acknowledges funding from the Royal Commission for the Exhibition of 1851. JB, MAM, and IME. gratefully acknowledge financial support by the German Research Foundation within the DFG Project ER463/14-1.
Competing Interests The Authors declare no Competing Financial or Non-Financial Interests.
Data Availability The authors declare that all data supporting the findings of this study are available within the paper and its supplementary information files.
Author Contribution LCR and JB conceptualised the project. JB, MAM and IE performed calculations using the orbital projective model. LCR and ME performed calculations using the 10-orbital tight binding model. All authors contributed to writing the manuscript.

SUPPLEMENTARY NOTE 1
In this note, we discuss the parametrization of states near Γ and Z point. We model the fermionic states at the center of the Brillouin zone using the effective low energy model introduced in Ref. [33]. In the tetragonal as well as in the orthorhombic phase the states near the Fermi level at the Γ and Z point can be described in terms of a two component spinor Ψ † Γ/Z,σ (k) = (d † xz,σ (k), d † yzσ (k)) with momentum k measured with respect to (0, 0, 0) at Γ and (0, 0, π) at Z point,respectively The nematic order parameter, which we discuss in the main text, and spin-orbit coupling enter the parametrization via Written in coordinates of the unfolded Brilloin zone the Hamiltonian is given by where we used c = − b 2 and µ is the chemical potential that fixes the number of particles accounting for Luttinger theorem. Exact diagonalization leads the band dispersions corresponding to the larger outer and smaller inner hole pocket h in and h out , respectively. Note that the inner hole pocket sinks below the Fermi level once h xz,yz − µ < ( Within our model we define superconductivity as momentum independent mean field order parameters in orbital space. The order parameters for d yz and d xz intra-orbital onsite pairing are given by which in Nambu notation Ψ † k = (Ψ † Γ/Z,↑ (k), Ψ T Γ/Z,↓ (−k)) yields the BdG Hamiltonian H with the pairing matrix∆ h = τ0+τ3 We perform a unitary transformation that diagonalizes H Γ/Z,↑ (k) from orbital into band space and obtain the BdG-Hamiltonian band basiŝ in which the intra-and inter-band gaps are given by respectively, and and In absence of SOC, diagonalization of Eq. (33) yields the band dispersion ξ reg,inc X/Y = A X/Y ± B 2 X2/Y2 + B 2 X3/Y3 with the "+" and "-" solution corresponding to a regular and incipient electron band at X and Y point, respectively. We define superconductivity as mean field order parameters in orbital space. The order parameters for d yz and d xz intra-orbital onsite pairing are given by Here we neglect the contribution of d xy orbitals to the pairing problem as discussed in the main text. In the Nambu basis with the pairing matrices∆ X (k) = τ0+τ3 2 ∆ e yz and∆ Y (k) = τ0+τ3 2 ∆ e xz . A unitary transformation U † k = α e xy X /xy Y (k) α e xz/yz (k) −α e * xz/yz (k) α e * xy X /xy Y (k) transforms H X/Y from orbital into band space with two eigenvalues of which ξ reg X/Y = A X/Y + B 2 X2/Y2 + B 2 X3/Y3 describe the regular bands with the corresponding eigenvectors |X reg,σ and |Y reg,σ . Following Ref. [33] we define the projector U F S = diag(|X reg,↑ , |Y reg,↓ , |X reg,↓ , |Y reg,↑ ) and treating SOC as a small perturbation we can transform Eq. (40) to the reduced basis which involves states near the Fermi level only.
where the intra-band gaps and the orbital content are given by and |a e xy X /xy Y (k)| 2 = sin 2 (φ |a e xz/yz (k)| 2 = cos 2 (φ The parameter determines the SOC induced splitting between inner and outer electron pocket. By performing another unitary transformation we bring Eq. (42) into band basis desribing the states at M and A point in presence of weak SOC.
The normal state band dispersion of Eq. (47) is given by ξ e in/out (k) = 2 |κ|) 2 and the superconducting intra-and inter-band order parameters can be parameterized as where For the gap at the peanut shaped pocket we find for small SOC ∆ e out (k) ≈ ∆ reg X (k). The fitting parameters for the dispersion near the A point listed in Tab In this note, we discuss the pairing interaction. We assume superconductivity to be driven by repulsive onsite interactions involving only intra-orbital interactions that translate mostly into pure inter-band pair-hopping interaction between hole and electron pockets and between both electron pockets corresponding to (π, 0), (0, π) and (π, π) spin-fluctuation mediated pairing, We treat Eq. (52) using a mean-field decomposition into the pairing terms introduced above: From now on we neglect contributions from inter-band pairing gaps (such as Eq. (25) and (50)) and write the free energy as where α ∈ {h in , h out , e in , e out , e inc1 , e inc2 } labels the bands in presence of SOC and∆ = ∆ h xz , ∆ h yz , ∆ e yz , ∆ e xz T the gaps in orbital space and E α (k) = ξ α (k) 2 + |∆ α (k)| 2 . The pairing matrix is given bŷ Minimizing the free energy with respect to∆ † and focusing on the larger hole pocket and the two electron bands near the Fermi surface (i.e {h in , e in , e out }) yields the BCS gap-equations Note that the system of equations (56) is not a closed one, as the chemical potential µ depends on the total number of particles (N) which we fix at 10 K in the non-superconducting orthohombic state as N (∆(T ), µ(T )) = N (0, µ(10K)) = − ∂F ∂µ = α,k A self-consistent treatment of the pairing problem requires solving the coupled set of equations (56) and (57).

SUPPLEMENTARY NOTE 4
In this note, we discuss the 10-orbital tight binding model used in part D of the results section in the main text. The tight binding model uses the form and notation of Ref.
[? ] and the hopping parameters, in units of eV, are 2D parameters Where V c/s = 1 2 U c/s χ c/s U c/s , χ c/s = χ 0 [1 ± U c/s χ 0 ] −1 is the RPA enhanced spin susceptibility in the charge (c) or spin (s) channel and In these calculations we assume spin rotational invariance, U = U − 2J, J = U 6 and J = J. U is set to 0.5 eV.