Fermi surface pockets in electron-doped iron superconductor by Lifshitz transition

The Fermi surface pockets that lie at the corner of the two-iron Brillouin zone in heavily electron-doped iron selenide superconductors are accounted for by an extended Hubbard model over the square lattice of iron atoms that includes the principal 3d xz and 3d yz orbitals. At half filling, and in the absence of intra-orbital next-nearest neighbor hopping, perfect nesting between electron-type and hole-type Fermi surfaces at the the center and at the corner of the one-iron Brillouin zone is revealed. It results in hidden magnetic order in the presence of magnetic frustration within mean field theory. An Eliashberg-type calculation that includes spin-fluctuation exchange finds that the Fermi surfaces undergo a Lifshitz transition to electron/hole Fermi surface pockets centered at the corner of the two-iron Brillouin zone as on-site repulsion grows strong. In agreement with angle-resolved photoemission spectroscopy on iron selenide high-temperature superconductors, only the two electron-type Fermi surface pockets remain after a rigid shift in energy of the renormalized band structure by strong enough electron doping. At the limit of strong on-site repulsion, a spin-wave analysis of the hidden-magnetic-order state finds a"floating ring"of low-energy spin excitations centered at the checkerboard wavenumber (pi,pi). This prediction compares favorably with recent observations of low-energy spin resonances around (pi,pi) in intercalated iron selenide by inelastic neutron scattering.


Introduction
The family of iron-based superconductors first discovered ten years ago has established a new route to high critical temperatures [1]. In particular, stoichiometric FeSe becomes superconducting below a modest critical temperature of 8 K. Electron doping of FeSe raises the critical temperature dramatically into the range 40-110 K, however [2,3]. The latter has been achieved by laying a monolayer of FeSe over a substrate [4,5,6,7], by intercalating layers of FeSe with organic compounds [8,9,10], by dosing thin films of FeSe with alkali metals [11,12,13], and by applying a gate voltage to thin films of FeSe [14,15]. Angle-resolved photoemission spectroscopy (ARPES) on such heavily electron-doped FeSe reveals two electron-type Fermi surface pockets at the corner of the two-iron Brillouin zone [5,6]. It also reveals buried hole bands at the center of the twoiron Brillouin zone that lie 60 meV below the bottom of the former electron bands [8,16]. At low temperature, ARPES also finds an isotropic gap at the electron Fermi surface pockets consistent with a superconducting state [8,17]. The energy gap at the Fermi level is confirmed by scanning tunneling microscopy on heavily electron-doped FeSe [7,10]. By comparison with stoichiometric FeSe, which has a relatively low critical temperature, it has been argued that the phenomenon of high-temperature superconductivity in heavily electron-doped FeSe is due to the appearance of a new electronic groundstate [18,19].
In addition, recent inelastic neutron scattering experiments on intercalated FeSe find a ring of low-energy magnetic excitations centered at the wave number (π/a, π/a) associated with Néel order over the square lattice of iron atoms in a single layer [20,21,22]. Here, a denotes the lattice constant. It has been suggested recently by one of the authors that this ring of low-energy magnetic excitations is a result of hidden Néel order among the iron 3d orbitals [18,23]. (See Fig. 6.) Such hidden magnetic order can emerge from frustration among local magnetic moments over the square lattice of iron atoms in FeSe [24,25]. One of the authors has shown that adding electrons to the local magnetic moments at half filling results in electron-type Fermi surface pockets at the corner of the two-iron Brillouin zone, but with no Fermi surface at the center [18]. It is important to mention, at this stage, that conventional band structure calculations for electron-doped FeSe typically predict an additional hole-type Fermi surface pocket at the center of the two-iron Brillouin zone, in marked contrast to ARPES measurements [8,17,26] In this paper, we shall introduce an extended Hubbard model over the square lattice of iron atoms in heavily electron-doped FeSe that harbors hidden Néel order among the 3d xz /3d yz orbitals because of perfect nesting of electron-type and of hole-type Fermi surfaces at the center and at the corner of the one-iron Brillouin zone, with nesting wavenumber (π/a, π/a). True Néel order is suppressed because of magnetic frustration from super-exchange interactions across the Se atoms. At half filling, mean field theory similar to that for the one-orbital Hubbard model over the square lattice [27] finds a stable hidden spin-density wave (hSDW) at the same wavenumber, but with a nodal D xy gap in the quasi-particle spectrum at the Fermi level. Also in analogy with the one-orbital Hubbard model over the square lattice [28,29,30], we identify two collective modes of the mean field theory that represent hidden spin-wave excitations. They vanish at the Néel wave number (π/a, π/a) with an acoustic dispersion in frequency. We argue in the Discussion section that degeneracy of these hidden Goldstone modes with true spin excitations results in a ring of low-energy magnetic excitations similar to what is observed by inelastic neutron scattering in intercalated FeSe [20,21,22].
Next, we shall formulate an Eliashberg theory [31,32,33] for the extended Hubbard model of a single layer of heavily electron-doped FeSe that is based on exchange of the above hidden spin-wave excitations by electrons and by holes (particle-hole channel). A solution of the associated Eliashberg equations [34] finds a Lifshitz transition of the electron/hole Fermi surfaces to pockets centered at the corner of the two-iron Brillouin zone at moderate to strong on-site Coulomb repulsion. This result is consistent with similar results obtained by one of the authors at the limit of strong on-site Coulomb repulsion [18], which predict electron-type Fermi surface pockets alone at the corner of the two-iron Brillouin zone at any level of electron doping. At strong but finite on-site Coulomb repulsion, the present Eliashberg-type calculation finds a threshold electron doping beyond which electron-type Fermi surface pockets appear alone. Below, we introduce the two-orbital Hubbard model for heavily electron-doped FeSe.

Perfect nesting of Fermi surfaces
We retain only the 3d xz /3d yz orbitals of the iron atoms in the following description for a single layer of heavily electron-doped FeSe. In particular, let us work in the isotropic basis of orbitals d− = (d xz − id yz )/ √ 2 and d+ = (d xz + id yz )/ √ 2. The kinetic energy is governed by the hopping Hamiltonian where repeated indices α and β are summed over the d− and d+ orbitals, where repeated index s sums over electron spin, and where i, j and i, j represent nearest neighbor (1) and next-nearest neighbor (2) links on the square lattice of iron atoms. Above, c i,α,s and c † i,α,s denote annihilation and creation operators for an electron of spin s in orbital α at site i. The reflection symmetries shown by a single layer of FeSe imply that the above intra-orbital and inter-orbital hopping matrix elements show s-wave and dwave symmetry, respectively [35,36,37]. In particular, nearest neighbor hopping matrix elements satisfy with real t 1 and t ⊥ 1 , while next-nearest neighbor hopping matrix elements satisfy with real t 2 and pure-imaginary t ⊥ 2 .
The above hopping Hamiltonian is easily diagonalized by plane waves of d x(δ)z and id y(δ)z orbitals that are rotated with respect to the principal axis by a phase shift δ(k): where N = 2N Fe is the number of iron site-orbitals. Their energy eigenvalues are respectively given by ε are diagonal and off-diagonal matrix elements, with It is notably singular at k = 0 and Q AF , where the matrix element ε ⊥ (k) vanishes.
Let us now turn off next-nearest neighbor intra-orbital hopping: t 2 = 0. Notice, then, that the above energy bands satisfy the perfect nesting condition where Q AF = (π/a, π/a) is the Néel ordering vector on the square lattice of iron atoms. The relationship (7) is an expression of a particle-hole symmetry that the hopping Hamiltonian (1) exhibits at t 2 = 0. (See Appendix A.) As a result, it can easily be shown that the Fermi level at half filling of the bands lies at ε F = 0. Figure 1 shows such perfectly nested electron-type and hole-type Fermi surfaces for hopping parameters t 1 = 200 meV, t ⊥ 1 = 500 meV, t 2 = 0 and t ⊥ 2 = 100 i meV. We shall now demonstrate how the perfectly nested Fermi surfaces shown by Fig. 1 can result in an instability to long-range hidden Néel order. It is useful to first write the creation operators for the eigenstates (4) of the electron hopping Hamiltonian, H hop : where α = 0 and 1 index the d− and d+ orbitals, and where n = 1 and 2 index the anti-bonding and bonding orbitals id y(δ)z and d x(δ)z . The inverse of the above is then It is then straight-forward to show that the spin magnetization for true (m = 0) or for hidden (m = 1) magnetic order, Figure 1. Band structure with perfectly nested Fermi surfaces at half filling, ε + (k) = 0 and ε − (k) = 0, with hopping matrix elements t 1 = 200 meV, t ⊥ 1 = 500 meV, t 2 = 0, and t ⊥ 2 = 100 i meV. Point nodes of quasi-particle gap are marked on Fermi surfaces.
The contribution to the static spin susceptibility from inter-band scattering that corresponds to true (m = 0) or to hidden (m = 1) Néel order is then given by the Lindhard function where n F is the Fermi-Dirac distribution. Next, application of the perfect-nesting condition (7) yields a more compact expression for the inter-band contribution to the static spin susceptibility (13): We conclude that the static susceptibilities for true and for hidden Néel order diverge logarithmically as χ inter (0, Q AF ) ∝ lim ǫ→0 c 2 D + (0) ln(W bottom /ǫ) and χ inter (1, Q AF ) ∝ . The one-iron Brillouin zone is divided into a 10, 000 × 10, 000 grid, while the δ-function is approximated by (4k B T 0 ) −1 sech 2 (ε/2k B T 0 ). Here, k B T 0 is 3 parts in 10, 000 of the bandwidth. The peak is a logarithmic van Hove singularity at ε + (π/a, π/a). lim ǫ→0 s 2 D + (0) ln(W bottom /ǫ), with corresponding density of states weighted by the magnitude square of the matrix element (12): Above, W bottom = −ε + (0, 0). The weighted densities of states (15) and (16) are of comparable strength at the Fermi level, ε = 0. For example, numerical calculations that are described in the caption to Fig. 2 yield the values c 2 D + (0) = 0.135 a −2 eV −1 and s 2 D + (0) = 0.072 a −2 eV −1 for these quantities. Here, hopping matrix elements coincide with those listed in the caption to Fig. 1. Hidden magnetic order is therefore possible.

Hidden magnetic order and excitations
We have just seen how the perfect nesting of the Fermi surfaces shown by Fig. 1 results in an instability towards long-range Néel order per d+ and d− orbital. Below, we shall introduce an extended Hubbard model over the square lattice that includes these orbitals alone. Within mean field theory, we shall see that long-range hidden Néel order exists at half filling because of magnetic frustration by super-exchange interactions.

Extended Hubbard model
We shall now add on-site interactions due to Coulomb repulsion and super-exchange interaction via the Se atoms to the hopping Hamiltonian (1). The Hamiltonian then has three parts: H = H hop + H U + H sprx . On-site Coulomb repulsion is counted by the second term in the Hamiltonian [38], where n i,α,s = c † i,α,s c i,α,s is the occupation operator, where S i,α = s,s ′ c † i,α,s σ s,s ′ c i,α,s ′ is the spin operator, and where n i,α = n i,α,↑ + n i,α,↓ . Above, U 0 > 0 denotes the intraorbital on-site Coulomb repulsion energy, while U ′ 0 > 0 denotes the inter-orbital one. Also, J 0 < 0 is the Hund's Rule exchange coupling constant, which is ferromagnetic, while J ′ 0 denotes the matrix element for on-site-orbital Josephson tunneling. The third and last term in the Hamiltonian represents super-exchange interactions among the iron spins via the selenium atoms: Above, J are positive super-exchange coupling constants over nearest neighbor and next-nearest neighbor iron sites. We shall assume henceforth that magnetic frustration is moderate to strong: J . In isolation, the above term of the Hamiltonian then favors "stripe" spin-density wave order at half filling over conventional Néel order.

Mean field theory
Following the mean-field treatment of antiferromagnetism in the conventional Hubbard model over the square lattice at half filling [27,28,29,30], assume that the expectation value of the magnetic moment per site, per orbital, shows hidden Néel order: where m i,α = 1 2 n i,α,↑ − 1 2 n i,α,↓ . The super-exchange terms, H sprx , make no contribution within the mean-field approximation, since the net magnetic moment per iron atom is null in the hidden-order Néel state. And we shall neglect the onsite Josephson tunneling term in (17) H U . This is valid at the strong-coupling limit, U 0 → ∞, where the formation of a spin singlet per iron-site-orbital is suppressed. We are then left with the two on-iron-site repulsion terms and the Hund's Rule term in H U .
The mean-field replacement of the intra-orbital on-site term (U 0 ) is the standard one [27]. In particular, replace The first term above can be absorbed into the chemical potential, the last term above is a constant energy shift, leaving a mean-field contribution to the Hamiltonian . The mean-field replacement of the inter-orbital oniron-site repulsion term (U ′ 0 ) in H U is also standard: The first two terms above can be absorbed into a shift of the chemical potential because n i,d− = n i,d+ , while the third and last term above is again a constant energy shift. The inter-orbital repulsion term, hence, makes no contribution to the Hamiltonian within the mean-field approximation. Last, we make the same type of mean-field replacement for the Hund's Rule term (J 0 ) in H U : Again, the last term above is just a constant energy shift. The first two terms, however, contribute to the mean-field Hamiltonian: in the case of hidden magnetic order (19).
The net contribution to the mean-field Hamiltonian from interactions in the present two-orbital Hubbard model is then Notice that the last sum above is simply twice the hidden (m = 1) ordered moment S z (m, Q AF ) defined by (10). Inspection of (11) then yields that the mean-field Hamiltonian for the present two-orbital Hubbard model takes the form with a gap function wherek = k + Q AF , and where Here, we have used the result (12) for the matrix element in the case of hidden magnetic order (m = 1). Here also, intra-band scattering (n ′ = n) has been neglected because it shows no nesting. After shifting the sum in momentum of the first term in (21) by Q AF for the anti-bonding band, n = 1, we arrive at the final form of the mean-field Hamiltonian: Above, we have set the ± sign in the matrix element (12) to minus for convenience. The mean-field Hamiltonian (24) is diagonalized in the standard way by writing the electron in terms of new quasi-particle excitations [28,29,30]: Above, u(k) and v(k) are coherence factors with square magnitudes where The mean-field Hamiltonian can then be expressed in terms of the occupation of quasiparticles by The quasi-particle excitation energies are then E(k) for particles and E(k) for holes. Notice that the gap (22) in the excitation spectrum has D xy symmetry. (See Fig. 1.) At half filling then, the energy band −E(k) is filled, while the energy band +E(k) is empty. Last, inverting (25) yields As expected, the quasiparticles are a coherent superposition of an electron of momentum k in the bonding band 2 with an electron of momentum k + Q AF in the anti-bonding band 1. Finally, to obtain the gap equation, we exploit the pattern of hidden Néel order (19), and equivalently write the gap maximum (23) as Using expression (11) and the result (12) in the case of hidden magnetic order (m = 1) yields the relationship wheren = 1 + (n mod 2). Again, we have neglected intra-band scattering. Also, notice that the sum over the bands (4) yields which is orbitally isotropic. Substituting in (25) and the conjugate annihilation operators, and recalling that the n = 1 quasi-particle band is filled in the groundstate, while the n = 2 quasi-particle band is empty, yields c † s (n,k)c s (n, k) = −(sgn s)u(k)v(k) for the expectation value. We thereby obtain the relationship or equivalently, the gap equation In the limit ∆ 0 → ∞, we then have In the thermodynamic limit, N Fe → ∞, integration along one of the principal axes followed by a series expansion in turn yields m 0,0 = 1/4 for the previous expression. In the limit |t ⊥ 1 | → 0, on the other hand, we have | sin 2δ(k)| → 1 by (6b), which yields an ordered magnetic moment m 0,0 = 1/2. Figure 3 shows the ordered magnetic moment at U 0 → ∞ versus hybridization of the 3d xz and 3d yz orbitals.
The above mean field theory predicts quasiparticles with excitation energies that disperse as E + (k) and E − (k), where E ± (k) = (ε 2 ± (k) + ∆ 2 0 [sin 2δ(k)] 2 ) 1/2 . They reach zero at point nodes located where the Fermi surface crosses a principal axis, at which the phase shift δ(k) is a multiple of π/2. These point nodes are shown by Fig. 1. The quasi-particle energy spectra disperse about the nodes in a Dirac-cone fashion: where v F is the Fermi velocity of ε ± (k) at the node, and where δ ′ D is the gradient of the phase shift δ(k) at the node. Here (k , k ⊥ ) denote the momentum coordinates about a zero point node in the directions parallel and perpendicular to the Fermi surface there. Last, notice that combining the spectra E + (k) and E − (k) in the folded Brillouin zone results in four Dirac cones at the Fermi level.

Low-energy collective modes
The groundstate of the above mean field theory is the filled energy band −E(k): Inspection of (28) yields that it can also be expressed as where |1 = k,s c † s (1, k)|0 is the filled anti-bonding band of non-interacting electrons. Next, observe that each pair of factors above per momentum over spin ↑ and ↓ can be expressed as u 2 (k) exp(−[v(k)/u(k)]n · s 1 ,s 2 c † s 1 (2,k)σ s 1 ,s 2 c s 2 (1, k)), with unit vector n along with the z axis. Now define a new spin quantization axis z ′ = y, along with the remaining axes x ′ = z and y ′ = x. If, more generally, we let the axis of the sub-lattice The one-iron Brillouin zone is divided into a 10, 000 × 10, 000 grid.
magnetizationn lie in the z-x plane, then the spin operator in the argument of the exponential above becomes Here, c ′ s and c ′ † s are the electron annihilation and creation operators in the new quantization axes. Here also, φ is the angle thatn makes with the z axis. Re-expanding the exponential operator above then yields the groundstate (31) in the new quantization axes: It has indefinite spin S y along the new quantization axis: where |Ψ (m) 0 are projections of the groundstate that have spin S y equal to m . Then because S y = i ∂ ∂φ , we have that the macroscopic phase angle and the macroscopic spin along the y axis satisfy the commutation relationship [40] [φ, S y ] = i .
Define, next, the macroscopic magnetization, M y = S y /V , where V is the area, and let it and the phase angle φ vary slowly over the bulk. Their dynamics is then governed by the hydrodynamic Hamiltonian [41,42] Above, χ ⊥ and ρ s denote, respectively, the transverse spin susceptibility and the spin stiffness of the present hidden spin-density wave state. Given the commutation relationship [φ(r), M y (r ′ )] = i δ (2) (r−r ′ ), we obtain the following dynamical equations: The magnetization thus satisfies the wave equationM y = c 2 0 ∇ 2 M y , with propagation velocity c 0 = (ρ s /χ ⊥ ) 1/2 . We conclude that the present hidden spin-density wave state supports antiferromagnetic spin-wave excitations that disperse acoustically in frequency: ω(k) = c 0 |k|. And since the above dynamics can be rotated by 90 degrees about the z axis, there then exist two acoustic spin-wave excitations per momentum.

Transverse spin susceptibility and spin rigidity
To compute the transverse spin susceptibility, we apply an external magnetic field along the y axis by adding the term −h i α S (y) i,α to the Hamiltonian H hop +H U +H sprx . The on-site-orbital repulsion terms, the Hund's Rule coupling terms, and the super-exchange terms can then be replaced by the isotropic mean-field approximations Yet the external magnetic field cants the antiferromagnetically aligned moments per orbital along the y axis by the transverse magnetization per orbital, m ⊥ . It makes a contribution to the above mean-field replacements that can be accounted for by making the replacement h → h + 2U(0) m ⊥ in the paramagnetic term that we added above to the Hamiltonian, where We thereby arrive at the formula for the transverse spin susceptibility per iron atom, where χ ⊥ is the naive transverse spin susceptibility that neglects the effect of canting.
The formula for the naive transverse spin susceptibility is well known [43], and it is derived in Appendix C. It reads At the weak-coupling limit, U(0), U(π) → 0, the transverse spin susceptibility (37) is given by χ (0) ⊥ , and the quotient in the sum over momentum above is equal to 2 δ[ε + (k)]. We thereby obtain the Pauli paramagnetic susceptibility at weak-coupling, where D + (ε) is the density of states of the bonding band ε + (k). (See Fig. 2.) At the strong-coupling limit, U 0 → ∞, it's useful to re-write the formula (37) as [43] Observe, next, the following identity for the quotient in (38): Applying the gap equation (30) then yields the result U(π)χ (0) In the limit U 0 → ∞, it can be shown that I 1 = I 2 as the pure-imaginary hopping matrix element t ⊥ 2 tends to zero. (See Appendix C.) In such case, U(π)χ (0) ⊥ = 1, and we thereby achieve the result It coincides precisely with the corresponding result for the Heisenberg model (67) [25] after making the assignments J 1 = J (sprx) 1 for two of the exchange coupling constants. Here J n and J ⊥ n represent intra-orbital and inter-orbital Heisenberg exchange coupling constants among the d+ and d− orbitals.
Next, to compute the spin stiffness ρ s at half filling, we follow the calculation of the same quantity in the case of the conventional Hubbard model over the square lattice [28,29,30,43,44]. At zero temperature, the spin rigidity saturates the f sum rule for the spin current. We therefore arrive at the expression for it. In the weak-coupling limit, where the gap function vanishes, we therefore get lim U (π)→0 ρ s = N −1 Fe ′ k t 1 (cos k x a + cos k y a), where the prime notation indicates the condition that ε + (k) < 0. In such case, the sum over momenta lies inside the Fermi surface that is centered at the center of the Brillouin zone. (See Fig. 1.) At strong coupling U(π) → ∞, on the other hand, it is useful to rewrite the above expression (43) as Here, we have substituted in the expression for the coherence factor (26). Approximate now all dispersions in energy about the Dirac nodes at the Fermi surface ε + (k) = 0: e.g.; , and sin 2δ(k) ∼ = 2δ ′ D k y , where the coordinates of the Dirac node are (k D , 0). After taking the thermodynamic limit, N Fe → ∞, and after cutting off the resulting integrals in momentum by k 1 ∼ k D /2 in both the x and in the y directions, we obtain the following result in the limit of strong coupling: In this limit, the spin stiffness at half filling therefore scales as ρ s ∼ [t 2 /U]ln[U/t] with the scale of the hopping matrix elements t, and with the scale of the on-site repulsion energy U. It is useful to compare the latter result for the rigidity of hidden magnetic order at strong on-site-orbital repulsion (45) with that obtained from the corresponding two-orbital Heisenberg model (67) [25]: for the remaining two exchange coupling constants.

Eliashberg theory
The previous mean-field approximation of the extended two-orbital Hubbard model for a single layer of heavily electron-doped iron-selenide predicts Dirac quasi-particle excitations at nodes where the Fermi surface crosses a principal axis. (See Fig. 1.) Below, we shall demonstrate how the Fermi surface at weak coupling experiences a Lifshitz transition to Fermi-surface pockets at the corner of the two-iron Brillouin zone as the on-site-orbital repulsion grows strong. We will achieve this by first formulating an Eliashberg theory for the extended Hubbard model in the electron-hole channel.

Hidden spinwaves and interaction with electrons
It was revealed in section 3.3 that the above mean field theory for the hidden Néel state of the Hubbard model over the square lattice harbors spin-wave excitations that collapse to zero energy at the Néel wavenumber Q AF . Also, the hidden Néel state of the corresponding Heisenberg model over the square lattice exhibits the very same hidden spin-wave excitations [25]. Consider then the propagator for hidden spinwaves: iD(q, ω) = 1 √ 2 m + (π) 1 √ 2 m − (π) | q,ω , where m ± (π) = m x (π) ± i m y (π). Here, m(π) = m d− − m d+ is the hidden magnetic moment. The propagator takes the form in the case of the above mean-field theory, as well as in the case of the linear spinwave approximation of the Heisenberg model (67) [25]. It shows a pole in frequency that disperses acoustically as ω b (q) = c 0 |q| about Q AF , where c 0 = (ρ s /χ ⊥ ) 1/2 is the hidden-spin-wave velocity, and whereq = q + Q AF . In the former case, s 1 is equal to the sub-lattice magnetization, m 0,0 , while s 1 is given by the electron spin in the latter case. Last, above, χ ⊥ is the transverse spin susceptibility of the hidden Néel state. Before continuing, it is important to mention that the Eliashberg theory to be developed below shows renormalized Fermi surfaces that continue to exhibit perfect nesting: ε + (k) = ν and ε − (k) = −ν. (See Fig. 4.) Here, ν is a relative energy shift between the bonding and anti-bonding bands that emerges naturally. Given the self-consistent nature of the Eliashberg theory [33,34], it is then natural to take the transverse spin susceptibility χ ⊥ and the spin rigidity ρ s above from the meanfield theory developed in the previous section, but with the previous relative shift in the chemical potentials included. That is accomplished by making the replacement ε + (k) → ε + (k) − ν in the gap equation (30) and in expressions (38) and (44) for χ (0) ⊥ and for ρ s .
The previous mean-field theory implies that the hidden spinwaves in question interact with independent electrons governed by the hopping Hamiltonian, H hop . This is evident from the mean-field form of the interaction (20) in isotropic form: − i α U(π)m i,α · 2S i,α . Keeping only the transverse contributions yields the interaction − i α U(π)(m + i,α S − i,α + m − i,α S + i,α ). Plugging in expression (9) and its conjugate for the electron creation and destruction operators yields the following interaction contribution to the Hamiltonian with hidden spinwaves: where q = k −k ′ is the momentum transfer, and where the matrix element above is the prior one for hidden order (m = 1). (See Appendix B.) As required, notice that H e−hsw is invariant under the particle-hole symmetry described in Appendix A: c s (n, k) ↔ ±ic † s (n,k) and m − (π, q) ↔ −m + (π, q). Because we shall use Nambu-Gorkov formalism [33,45,46] below, it is useful to write the above electron-hidden-spinwave interaction in terms of spinors: with spinor Above, τ 1 is the Pauli matrix along the x axis. Also, the explicit matrix element M n,k;n,k ′ has been substituted in. (See Appendix B.)

Electron propagator and Eliashberg equations
Let C s (k, t) denote the time evolution of the destruction operators (49) C s (k), and let C † s (k, t) denote the time evolution for the conjugate creation operators C † s (k). The Nambu-Gorkov electron propagator is then the Fourier transform iG(k, ω) = dt 1,2 e iωt 1,2 T [C s (k, t 2 )C † s (k, t 1 )] , where t 1,2 = t 2 −t 1 , and where T is the time-ordering operator. It is a 2 × 2 matrix. In the absence of interactions, its matrix inverse is then given by where τ 0 is the 2 × 2 identity matrix, and where τ 3 is the Pauli matrix along the z axis. Guided by the previous mean field theory, let us next assume that the matrix inverse of the Nambu-Gorkov Greens function takes the form Here, Z(k, ω) is the wavefunction renormalization, ∆(k) is the quasi-particle gap (22), and ν is a relative energy shift of the bands that preserves perfect nesting. In particular, the form (51) of the Nambu-Gorkov Greens function is consistent with the perfect nesting that is equivalent to (7). Matrix inversion of (51) yields the Nambu-Gorkov Greens function [33,45,46] and G 2 = 0. Above, the excitation energy is To obtain the Eliashberg equations, recall first the definition of the self-energy correction: G −1 = G (0)−1 − Σ. Comparison of the inverse Greens functions (50) and (51) then yields [33,34] for it. Next, we neglect vertex corrections from the electron-hidden-spinwave interaction (48). This approximation will be justified a posteriori at the end of the next subsection. The self-energy correction is then approximated by with iω m = iω n − iω n ′ , and with q = k −k ′ . Here, we have Wick rotated to pure imaginary Matsubara frequencies at non-zero temperature T . Observe, finally, that τ 1 τ µ τ 1 = sgn µ τ µ , where sgn 0 = +1 = sgn 1 , and where sgn 2 = −1 = sgn 3 . Identifying expressions (54) and (55) for the self-energy corrections then yields the following selfconsistent Eliashberg equations at non-zero temperature: The Greens functions above are listed in (52). Last, the above Eliashberg equations can be expressed at real frequency. In particular, it becomes useful to write the propagator for hidden spinwaves (46) as A series of decompositions into partial fractions followed by summations of Matsubara frequencies yields Eliashberg equations in terms of Fermi-Dirac and Bose-Einstein distribution functions at real frequency. They are listed in Appendix D. At zero temperature, these reduce to Below, we find solutions to the above equations.

Fermi-surface pockets at corner of Brillouin zone
Let us now work within the approximation that the gap function ∆(k) shows no frequency dependence. The negative sign in the gap equation (58) implies then that ∆(k) is null. Henceforth, we shall then set ∆(k) = 0. Furthermore, let us neglect any angular dependence acquired either by the wavefunction renormalization, Z(k, ω), or by the relative energy shift of the bands, ν, on momentum around the Fermi surface: ε + (k) = ν. This is exact for ν near the upper band edge of ε + (k) in the absence of nearest-neighbor intra-orbital hopping, t 1 = 0, in which case the Fermi surface pockets at (π/a, 0) and at (0, π/a) are circular. Following the standard procedure [34], we then multiply both sides of the remaining two Eliashberg equations (58) by δ[ε + (k)−ν]/D + [ν] and integrate in momentum over the first Brillouin zone. The previous Eliashberg equations (58) thereby reduce to where and where the wavefunction renormalization is averaged over the new Fermi surface: . Above, we have approximated the function U 2 F 0 (Ω; ν, ε ′ ) of ε ′ by its value at the renormalized chemical potential, U 2 F 0 (Ω; ν, ν). It is also understood in (60) that the limit implicit in the last δ-function factor is taken last.
The effective spectral weight of the hidden spinwaves, U 2 F 0 (Ω; ν, ν), can be evaluated by choosing coordinates for the momentum of the electron, (k , k ⊥ ), that are respectively parallel and perpendicular to the Fermi surface of the bonding band (FS + ): ν = ε + (k). (See Figs. 1 and 4.) This yields the intermediate result where v = ∂ε + /∂k is the group velocity. Yet the dispersion of the spectrum of hidden spinwaves ω b (q) is acoustic at low energy. This then yields the limiting dependence on frequency for their effective spectral weight, U 2 F 0 (Ω; ν, ν) = ǫ E (ν)/Ω as Ω → 0, with a constant pre-factor Above, c 0 is the velocity of hidden spinwaves at Q AF .
We can now find solutions to the remaining Eliashberg equations (59a) and (59b). In particular, assume that the relative energy shift ν lies near the upper edge of the bonding band ε + (k) at (π/a, 0) and at (0, π/a). Figure 4 displays the Fermi surfaces in such case. Substituting in the simple pole in frequency above for U 2 F 0 (Ω; ν, ν) yields the first Eliashberg equation: Here, we have reversed the order of integration. Also above, [ν − W, ν] is the range of integration over ε ′ in (59a), where W = W bottom + W top is the bandwidth of ε + (k), while ω uv is an ultra-violet cutoff in frequency for the hidden spinwaves. Assuming W/Z ≫ Ω + ω yields the equation in the low-frequency limit, where x = Ω/ω above. The final result for the wavefunction renormalization at the Fermi surface is then Z = (π 2 /4)(ǫ E /ω) as ω → 0. By (52), the spectral weight of quasi-particle excitations is 1/Z. It therefore vanishes at the Fermi level, ω = 0. This result is then consistent with the characterization of the hSDW state as a Mott insulator. The second Eliashberg equation (59b) can be evaluated in a similar way. After substituting in the simple pole in frequency for U 2 F 0 (Ω; ν, ν), integrating first over ε ′ yields the equation Assuming, once again, the inequalities W /Z ≫ Ω and W/Z ≫ ω then yields the following equation in the low-frequency limit: Here, we have expanded the previous argument of the logarithm in powers of Ω. The final result for the second Eliashberg equation is then where ω ir ∼ c 0 /L is an infra-red cutoff in frequency. Above, the previous result from the first Eliashberg equation has been substituted in. We can now check the previous ε ± (k x ,k y ) = ± ν , d xz ORBITAL inequalities that were assumed. The second Eliashberg equation (65) implies that the energy scale (62) is of order ǫ E ∼ ε + (π/a, 0)/ln(ω uv /ω ir ) for ν near the upper edge of the bonding band ε + (k). The ratio W/ωZ ∼ W/ǫ E is therefore of order ln(ω uv /ω ir ), which diverges logarithmically as the infra-red cutoff in frequency ω ir tends to zero.
We shall finally estimate the Eliashberg energy scale (62) ǫ E . For simplicity, assume small circular Fermi surface pockets (t 1 = 0) of Fermi radius k F , which is related to the concentration x 0 of electron/holes per pocket by k F a = (2πx 0 ) 1/2 . The Fermi velocity is then v F = 2t ⊥ 1 k F a 2 . Also, the phase shift (6b) is approximately sin 2δ(k) = [(t ⊥ 2 /i)/2t ⊥ 1 ](k F a) 2 sin 2φ, where φ is the angle that k makes about the center of the Fermi surface pocket at (π/a, 0) or at (0, π/a). We then get the expression for the Eliashberg energy scale (62). Comparing this estimate with the second Eliashberg equation (65), while fixing ν to the upper edge of the band ε + (k), then yields that the area of the electron/hole Fermi surface pockets shown in Fig. 4 vanishes logarithmically with the size L of the system for any positive U(π). Note, however, that such behavior is effectively ruled out at the weak-coupling limit, U(π) → 0, because the ordered magnetic moment s 1 vanishes exponentially in such case. On the other hand, if instead the scale L of the system is fixed, then the previous comparison of (65) and (66) yields U(π) ∝ x −3/4 0 . By the previous estimate for the phase shift at the new Fermi surface pockets, the electron-hidden-spin-wave interaction (48) then scales as U(π) −1/3 . This justifies our neglect of vertex corrections at the limit of strong on-site repulsion, U 0 → ∞.
Last, a self-consistent solution to the above Eliashberg equations at the limit of large on-iron-site-orbital Coulomb repulsion, U 0 , also exists at a relative shift of the bands ν near the bottom edge of the bonding band, ε + (k), instead. In particular, [ν, ν + W ] is now the range of integration over ε ′ in the Eliashberg equations (59a) and (59b). The previous results for the wavefunction renormalization (64) and for the relative energy shift between the two bands (65) hold after making the replacement ν = W top with −ν = W bottom in the latter. It is important, now, to observe that the density of states of the bonding band at the upper band edge is larger than the density of states at the bottom edge by Fig. 2. The condensation energy is of order −D + (ν)∆ 2 0 , however. By the definition (23) for ∆ 0 and by Fig. 3 for the ordered magnetic moment, the condensation energy dominates the kinetic (hopping) energy at strong on-site repulsion U 0 compared to the bandwidth. This argues in favor of the former solution in such case, with ν at the upper edge of the band ε + (k).

Discussion
The previous mean field theory analysis of the extended two-orbital Hubbard model for heavily electron-doped FeSe finds that hidden Néel antiferromagnetic order is expected at perfect nesting (7) t 2 = 0 when true Néel order is suppressed by magnetic frustration. (See Figs. 1 and 4.) Below, we compare the observable consequences that have been listed above with analogous theoretical results at the strong-coupling limit [18], U 0 → ∞, and with recent experimental evidence for such hidden magnetic order in the superconducting state of intercalated FeSe [20,21,22].

Comparison of weak coupling and strong coupling
In subsections 3.3 and 3.4, we showed how mean field theory for the hidden Néel state of the two-orbital Hubbard model that describes heavily electron-doped FeSe agrees both qualitatively and quantitatively with the corresponding Heisenberg model at large on-iron-site-orbital repulsion. In particular, a hydrodynamical analysis (35) predicts two acoustically dispersing spin-wave excitations per momentum at the "checkerboard" wavenumber Q AF . This agrees with the large-s 0 analysis of the corresponding Heisenberg model [25]. Second, the transverse spin susceptibility of the hSDW state (37) was calculated above as well. After tuning the inter-orbital next-nearest neighbor hopping t ⊥ 2 , which is pure imaginary, towards the limit of weak hybridization between the 3d xz and 3d yz orbitals, the transverse susceptibility of the hSDW state at the limit of strong on-iron-site-orbital Coulomb repulsion (42) is found to agree with the same quantity calculated from the corresponding Heisenberg model in the large-s 0 limit [25]. Also, the spin rigidity (43) of the hSDW state was computed above at the limit of strong on-iron-site-orbital repulsion. A comparison with the same results for the corresponding Heisenberg model yields exchange coupling constants that are consistent with hidden Néel order. Figure 4 is the central result of the paper, however. It shows the Fermi surfaces of the extended Hubbard model in the hSDW state at half-filling and at strong on-ironsite-orbital Coulomb repulsion, as predicted by Eliashberg theory in the particle-hole channel. The rigid-band approximation, in turn, predicts electron-type Fermi surface pockets at wavenumbers (π/a, 0) and (0, π/a) upon electron doping at concentrations x > x 0 . Here, x 0 denotes the concentration of electrons/holes inside the Fermi surface pockets shown in Fig. 4. It vanishes as U 0 diverges. This result agrees with Schwinger-boson-slave-fermion mean field theory of the corresponding local-moment (t-J) model at electron doping [18], in which case U 0 → ∞ and x 0 → 0. It also notably agrees with ARPES on heavily electron-doped FeSe [5,6,8,9]. In particular, x 0 may represent a threshold concentration of electron doping at which hSDW order gives way to superconductivity.

Comparison of hidden magnetic order with experiment
The Heisenberg model that corresponds to the present extended Hubbard model for heavily electron-doped FeSe reads [25] where the α = d− or d+. The results obtained in subsection 3.4 for the transverse spin susceptibility and for the spin rigidity of the hSDW state are consistent with Heisenberg exchange coupling constants that satisfy J 1 > J ⊥ 1 and J 2 = J ⊥ 2 . At this strong-coupling limit, U 0 → ∞, adding electrons with only inter-orbital nearest neighbor hopping, t ⊥ 1 > 0, yields a local-moment model that has been recently analyzed within the Schwinger-boson-slave-fermion formulation [18]. The mean field theory of the corresponding hSDW state is well behaved. As mentioned previously, it shows two electron-type Fermi surface pockets at the corner of the two-iron Brillouin zone, with 3d xz and 3d yz orbital character respectively. (Cf. Fig. 4.) Schwinger-bosonslave-fermion mean field theory also finds two branches of spin-wave excitations that correspond to true and to hidden magnetic moments, S i,d− + S i,d+ and S i,d− − S i,d+ , respectively. They are governed by the Heisenberg model (67) in the large-s 0 limit, but with the replacement J ⊥ Here, s 0 is the electron spin, and x is the electron doping concentration per site-orbital. Figure 5 shows the spin-wave spectra from such a large-s 0 approximation for the hSDW state of the local-moment model near a critical Hund's Rule coupling J 0c where the spectrum softens completely at "stripe" SDW wavenumbers (π/a, 0) and (0, π/a) [18]; i.e., ∆ cSDW → 0.
The results shown by Fig. 5 for the spin-excitation spectrum of the hSDW state are obtained from the local-moment model for heavily electron-doped FeSe that includes only inter-orbital nearest neighbor hoping, t ⊥ 1 > 0. The 3d xz and 3d yz orbitals are good quantum numbers in such case. Likewise, the true and hidden magnetic moments just cited are also good quantum numbers. Unfortunately, unlike the previous analysis of the extended Hubbard model, the Schwinger-boson-slave-fermion mean field theory that such results are based on is not well behaved when mixing between the two orbitals (pure imaginary t ⊥ 2 ) is turned on [18]. Figure 6 shows, however, the points in momentum and energy at which the two branches of the spin-excitation spectrum are degenerate. It reveals a "floating ring" of low-energy magnetic excitations about the Néel wave number Q AF . Observable spin excitations should exist along the floating ring at weak mixing between the 3d xz and 3d yz orbitals. Similar low-energy spin resonances around Q AF have been observed recently in the superconducting phase of intercalated FeSe by inelastic neutron scattering [20,21,22]. Such experiments are then consistent with the hSDW state proposed here for heavily electron-doped FeSe.

Summary and Conclusions
Understanding the mechanism behind the high-temperature superconductivity displayed by heavily electron-doped iron selenide remains elusive. In an attempt to solve this mystery, we have shown how the electron-type Fermi surface pockets that exist at the corner of the two-iron Brillouin zone in heavily electron-doped iron-selenide can emerge from an extended Hubbard model over the square lattice of iron atoms that includes only the 3d xz and 3d yz orbitals. At half-filling, and in the absence of next-nearest neighbor intra-orbital hopping, perfect nesting exists between hole-type and electrontype Fermi surfaces displayed by Fig. 1. The nesting wavenumber is (π/a, π/a), which corresponds to checkerboard (Néel) order. It notably differs from parent compounds to iron-pnictide high-temperature superconductors, which display "stripe" spin-density order, with nesting vector (π/a, 0). The former checkerboard nesting leads to hidden magnetic order at that ordering wavenumber when true Néel order is suppressed by magnetic frustration. An extended Hartree-Fock calculation of the Eliashberg type reveals that hole and electron Fermi surfaces become centered at the corner of the two-iron Brillouin zone at moderate to strong on-site Coulomb repulsion because of the exchange of antiferromagnetic spin fluctuation. The electron/hole concentration x 0 that corresponds to the area of these Fermi surface pockets vanishes as the on-site Coulomb repulsion diverges. Sufficiently strong electron doping x > x 0 will then produce a rigid shift of such a renormalized band structure, with electron Fermi surface pockets that are characteristic of heavily electron-doped FeSe.
We have also shown that the extended two-orbital Hubbard model leads to a localmoment model in the limit of strong on-site Coulomb repulsion that harbors the same type of hidden magnetic order [25]. Recent calculations by one of the authors also find electron-type Fermi surface pockets at the corner of the two-iron Brillouin zone when electrons are added to the local moments [18]. Furtheremore, in the previous section, we have pointed out that the low-energy spin excitations predicted by the local-moment model in the hidden magnetic order phase resemble the "floating" ring of spin-excitations that has been observed recently in heavily electron-doped FeSe by inelastic neutron scattering [20,21,22]. The extended two-orbital Hubbard model therefore is promising phenomenologically. Recent exact calculations on finite clusters at the local-moment limit by one of the authors find evidence for Cooper pairs near a quantum critical point to "stripe" spin-density wave order [18]. It remains to be seen if the extended Hubbard model introduced here also harbors superconductivity, and if so, of what type. the limit lim t ⊥ 2 →0 sin 2δ(k) = 0. Next, observe by (41a) that the limit lim U 0 →∞,t ⊥ 2 →0 I 1 is equal to Fe k ([ε + (k)/U(π)] 2 + [ m 0,0 sin 2δ(k)] 2 ) −1/2 . (C.5) Figure 3 indicates that the hidden magnetic moment m 0,0 vanishes roughly as the hybridization between the 3d xz and 3d yz orbitals, |t ⊥ 2 |. We conclude that the limit lim U 0 →∞,t ⊥ 2 →0 I 1 diverges at least linearly with U 0 . Second, observe that the quotient in expression (41b) for I 2 is equal to 2 δ(∆ 0 [sin 2δ(k)]) in the limit U 0 → ∞. By (23) which coincides with the product of U(π) with the density of states of ∆ 0 [sin 2δ(k)] at zero energy. Now notice by (6b) that sin 2δ(k) disperses hyperbolically near (π/a, 0) and (0, π/a). This implies that lim U 0 →∞,t ⊥ 2 →0 I 2 diverges roughly as (|t ⊥ 1 |/|t ⊥ 2 |) 2 ln[U(π)/W bottom ]. Equating I 1 with I 2 then yields that t ⊥ 2 → 0 as U 0 → ∞, which is consistent with the original assumption.
Appendix D. Eliashberg equations at non-zero temperature Equation (56) in the text lists the three Eliashberg equations at non-zero temperature in terms of sums over Matsubara frequencies. The sums can be evaluated in closed form after a series of decompositions into partial fractions. That procedure yields Above, q = k − k ′ − Q AF . Also, n F (ε) and n B (ω) denote the Fermi-Dirac and the Bose-Einstein distribution functions.