Non-local dxy 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 dxy orbital strongly contributes to nematic ordering in FeSe. We show that the addition of dxy nematicity to a pure dxz/dyz 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 sulfur 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 FeSe1−xSx 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 iron-based 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 iron-based 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 twoelectron 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 as a function of temperature and sulfur 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.

RESULTS
The B 1g 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 in-plane on a square grid, and are connected by staggered out-of-plane pnictide or chalcogenide atoms, which gives rise to a two-Fe atom unit cell, as shown in Fig. 1a. This 2-Fe unit cell has primitive inplane 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. 1b. 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  and (R x , R y ), respectively. Black dots denote Fe atoms. Se atoms located above (filled green) and below (empty green) the iron plane divide the system into sub-lattices A and B. b Fermi surface cut at k z = π 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 .
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. 1b. 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. 1c. In Fig. 1d 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 onsite interactions on the Fe atoms or from interactions between two neighbouring Fe atoms on different sub-lattices. We list these four order 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 onsite d xz and d yz states. This form, often referred to as ferro-orbital ordering 35 , 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. 2a. This type of nematic order was originally discussed in the literature 15,[36][37][38] ; 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,36,[39][40][41][42] . 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,40,41 . 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,40,43,44 . 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 s-wave 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. 2b and 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,43 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. 2d.
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 Table 1. All B 1g nematic order parameters for FeSe.
All possible nematic order parameters within the B 1g irreducible representation of the P4/nmm space group involving interactions between d xz , d yz and d xy 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 Brillouin 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 d xy nematicity is a purely non-local phenomenon.
In Fig. 3a-c, we present the Fermi surface and band dispersion around the Ζ and A point for the tetragonal state, i.e. Φ h xz;yz ¼ Φ e xz;yz = Φ 4 = 0. In Fig. 3d-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 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,45 . In Fig. 3g-i we now present our scenario, where we use the same values of Φ h xz;yz and Φ e xz;yz as in Fig. 3d-f but now we set Φ 4 = 45 meV. We note that the Φ 4 term, as it is discussed in Table 1, 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 onsite d xy energy in the nematic state Δϵ xy d y xy;kþQY σ d xy;kþQY σ þd y xy;kþQX σ d xy;kþQX σ 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, while 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. 3i. 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 self consistently, 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 mean-fieldlike 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. 3c, 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 refs. 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,45,46 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,46,47 .

Experimental evidence for non-local d xy 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 twoelectron pocket Fermi surface in the tetragonal state into a oneelectron 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. 3g-i, with the experimental ARPES data of detwinned crystals of FeSe from ref. 19 . In Fig. 4d 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. 4a as well as the data published in ref. 20,21,24 . In Fig. 4e, 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. 4b, c, as well as the data published on detwinned crystals of FeSe in refs. 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 48,49 . This electronic structure is also fully consistent with ARPES measurements on "twinned" crystals 11,12,16,38,40,41,50 .
Despite the removal of an entire electron pocket within the nematic state, we find that the Luttinger theorem is conserved within this model. In Fig. 4g 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 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 38,51 . 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. 4h 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~7 0 K, shortly after the nematic transition at T s = 90 K.
Finally, in Fig. 4j we plot the temperature dependence of the spectral function, Tr f Im Gðk; ωÞgfðω; TÞ; at the A point, where G(k, ω) = [ω−(H(k)−μ)+iΓ] −1 and f(ω, T) is the Fermi 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 refs. 16,20,40,41 . We show the experimental measurement of ref. 16 , in Fig. 4i 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,41,51 . 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,40,41,49 . Coupled with the observation of d xz spectral weight of the band at −50 meV 20,21,40 , 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 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. 5b. Indeed, it yields the same result as shown in Fig. 5a with excellent agreement on temperature dependence and band position. This agreement  19 . b, c Experimental band dispersion along the k x and k y 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 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.
L.C. Rhodes et al. 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 domain 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,52,53 . 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, Shubnikov-de Haas oscillation measurements of the sulphur doped series have been reported by Coldea et al. 54 . 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. 6a, c) to the δ and γ frequencies of ref. 54 and the M and Γ pockets (labelled as 3, 4 in Fig. 6e, 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. 6i. 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% sulfur 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  20 , Spectral intensity divided by the Fermi function taken at +10 meV above E F as a function of temperature. The disappearance of a band can be inferred at T = 60 K. b Simulated evolution of the momentum of states near the A point, along the k x axis k = (k x , π, π), which have an energy of E =+10 meV as a function of temperature. j Corresponding experimental data taken from ref. 54 . We have included a dashed circle in i, j to highlight the evidence for a re-entering d xy band. Note that we use a different color code in a-i than the authors of ref. 54 and that we assign our pockets differently to the data and trend lines in j.
L.C. Rhodes et al. extremal areas in the tetragonal state to seven. The reason for this sharp increase is two-fold, the first is the sudden change of chemical potential as discussed in Fig. 4g, 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. While 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 similar (and very small) values. Although a Lifshitz transition of the electron band was not the interpretation of ref. 54 , 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. 54 in Fig. 6j for comparison. Very similar results can also be observed by suppressing nematic order under application of physical pressure to FeSe 55 . 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 inter-band 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 [56][57][58] . 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-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. 59 , 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 0 he and J 0 ee pairings that results in inter-band 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 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 inter-band, 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 0 ee ), which is considerably weakened if only oneelectron 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 meanfield 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 0 he ¼ J 0 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 Cooperpairing and integrate within an energy cutoff of 50 meV around the Fermi level. Further details are given in Supplementary note 3. 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 ja h xz=yz ðkÞj 2 ¼ jhd xz=yz jh out ij 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 with ja e yz=xy ðkÞj 2 ¼ jhd yz=xy je X ij 2 the d yz/xy orbital weight on the peanut.
In Fig. 7a we show the orbital weight ja h xz=yz ðkÞj 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. 7b 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. 43 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 ja e yz=xy ðkÞj 2 of the d yz and d xy orbitals on the peanut shaped electron pocket is shown in Fig. 7b. 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. 7c we show that jΔ h yz j>jΔ e yz j, which leads to max ðjΔ h out ðk F ÞjÞ> max ðjΔ e out ðk F ÞjÞ as it is seen in experiment. In Fig. 7e 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-14 , of which we show data extracted from STM measurements of ref. 10 , depicted as red and blue circles in Fig. 7e. For completeness, in Fig. 7d we present the total density of states (DOS) ρðωÞ ¼ À 1 π P k Tr τ0þτ3 2 ImĜðω; kÞ calculated for the gap values that correspond to Fig. 7d. The DOS is nearly V-shaped, 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. 8d. Here, we model the doping dependence of the nematic order parameters , 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. 8b, 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 sulfur doping 60,61 . 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 y 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 equation This term ensures the asymmetric splitting of the two d xy bands around the M (A) point are captured while 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 Supplementary Note 4, 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 62 , Here, we integrate the pairing vertex, Γ μν ðk; k 0 Þ overall k 0 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 (Although in principle the quasiparticle weight of each orbital can be different, here we assume that the difference is small enough that it can be considered negligible).
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,63 . The exact form of the pairing vertex is described in Supplementary Note 5; 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 spinfluctuation frequency and E F is the Fermi energy 64 , which is not true for FeSe where the Fermi energy is <10 meV 38 .
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 iron-based superconductors 65,66 . 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. 9a 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,44,46 , 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,67 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. 9b. This is equivalently the case for the phenomenological spin susceptibility in Fig. 9e. 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. 9c. However, if (0, π) and (π, 0) fluctuations are assumed to be dominant then s ± becomes the leading instability (Fig. 9f).
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.

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,68 , 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 xy -dominated electron pocket. This was indirectly supported by DMFT calculations, arguing about the stronger correlations effect on the d xy orbital 69 , yet we are not aware of DMFT calculations that would correctly reproduce the measured oneelectron 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 sulfur 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 70 , 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 Fig. 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 k x -axis. e x and e y define the momentum dependence of the gap at (π, 0) and (0, π) respectively, while h o and h i 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 0 ¼ U 6 , J = 0.1U and J 0 ¼ J throughout and use the unfolded 1-Fe unit cell 34 . Further details can be found in Supplementary note 5. specific heat and STM measurements 60,61 . 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 xy -nematicity 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 71 . Observe also that this order needs to be re-examined within the RG weak-coupling analysis performed previously 29,30,72 . 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 41,49,73 . Although a recent DFT calculation proposed an E u symmetry nematic order parameter, using symmetry preconditioned wavefunctions 74 , 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 non-local 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 75,76 .
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]77 . 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,47 ; 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. 78 also addressing the question of a single electron pocket Fermi surface of FeSe.

Superconducting gap simulations for the orbital projective model
Our simulations of the superconducting gap for the orbital projective model follow the methodology introduced in ref. 43 details of how this calculation was performed, as well as the parameters used in this model can be found in Supplementary Material 1-3.
Superconducting gap simulations for the 10-orbital tight binding model using the 10-orbital tight binding model. All authors contributed to writing the manuscript.

FUNDING
Open Access funding enabled and organized by Projekt DEAL.