Magnetic fluctuations and spin-spirals in single-layer FeSe

The magnetic properties of monolayer FeSe films are investigated via first-principles spin-spiral calculations. Although the (pi,pi) collinear antiferromagnetic (CL-AFM) mode is lowest in energy, the spin-wave energy E(q) - which exhibits intrinsic non-Heisenberg behavior - is found to be extremely very flat over a large region of the two-dimensional Brillouin zone centered at the checkerboard antiferromagnetic (CB-AFM) q=0 configuration, giving rise to a sharp peak in the spin density of states. Considering the paramagnetic state as an incoherent average over spin-spiral states, we find that resulting electronic band states around the Fermi level closely resemble the bands of the CB-AFM configuration - not the CL-AFM one - and thus providing a natural explanation of the angle-resolved photoemission observations. The presence of the SrTiO3(001) substrate, both with and without interfacial oxygen vacancies, is found to reduce the energy difference between the CB-AFM and CL-AFM states and hence enhancing the CB-AFM-like fluctuations.

Single-layer FeSe films grown on SrTiO 3 (001) (STO) have generated intense interest because of their reported high superconducting critical temperature T C ∼40-100 K. [1][2][3][4][5][6][7][8][9][10] Angle-resolved photoemission spectroscopy (ARPES) experiments 4-10 reveal a number of distinct features: As-prepared single-layer FeSe/STO films are not superconducting, but are insulating or semiconducting. After intensive vacuum annealing, however, the monolayer films are superconducting with a metallic Fermi surface. 6,7 Contrary to most iron-based superconductors, this Fermi surface is characterized by electron pockets centered at the Brillouin zone corner (M), with the zone center (Γ) states pushed below the Fermi level, posing a challenge for pairing theories relying on Fermi surface nesting between the Γ and M-centered pockets 11,12 through (π,π) collinear stripe antiferromagnetic (CL-AFM) spin fluctuations (c.f., Fig. 1(c) inset). Oxygen vacancies formed at the interface during the annealing process are believed to play an important role in the transition to the metallic phase and high T C superconductivity. Furthermore that bilayer FeSe films are not superconducting 7 suggests that the FeSe-STO interface may play a key role in the superconducting mechanism.
A number of first-principles density-functional theory (DFT) calculations have been reported for FeSe thin films. 10,[13][14][15][16][17][18][19][20][21][22][23] Detailed comparisons of the experimental ARPES data 10 for single-layer FeSe/STO to the DFT calculations corresponding to non-magnetic and various ordered antiferromagnetic configurations find that the calculated bands of only the checkerboard antiferromagnetic (CB-AFM) configuration ( Fig. 1(c) inset) are consistent with the experimental data: Both the Fermi surface and the band structure throughout the whole Brillouin zone (BZ) are reasonably reproduced, particularly when a small Hubbard U correction is included that pushes the Γ-centered hole-like band completely below the Fermi level. Other DFT studies have examined the effect of interfacial oxygen vacancies (O-vac): 13-17 oxygen vacancies are found to electron-dope the FeSe layer, and to modify the band structure of the CB-AFM state, including the hole-like band around Γ so that satisfactory agreement with the ARPES data is achieved without the addition of a phenomenological U . 17 This seemingly successful agreement between the ARPES data and the DFT calculated bands for the CB-AFM configuration of monolayer films of FeSe/STO, however is apparently inconsistent with the fact that the CB-AFM configuration is not the calculated ground state. Rather, DFT total energy calculations consistently find that the CL-AFM state is lower in energy than the CB-AFM state for both monolayer FeSe and FeSe/STO, but the calculated CL-AFM bands do not resemble the ARPES ones. Although the energy difference between the CB-and CL-AFM states is found to be reduced by electron doping due to the substrate, the CL-AFM state remains more stable at reasonable doping levels. 14,17,22 Moreover, there is no experimental evidence for longrange magnetic order in either bulk or monolayer FeSe -recent inelastic neutron scattering experiments 24 have found that CB-AFM and CL-AFM correlations coexist even at low excitation energies in the bulk -suggesting that a quantum paramagnetic state with strong fluctuations between the CB-AFM and CL-AFM states might also exist in the monolayer system. Furthermore, it has been suggested that a nematic paramagnetic state resembling the observed bulk FeSe nematic state is a result of a near degeneracy between the CB-AFM and CL-AFM states. 25 Even in the case of coexisting CB-AFM and CL-AFM correlations, the question of why ARPES experiments measure the CB-AFM-like band structure remains.
The apparent paramagnetic nature of FeSe implies (i ) there is no net magnetic moment or long-range antiferromagnetic order, and (ii ) the translational (and space group) symmetry is consistent with the crystallographic 4 atom unit cell. These conditions, however, do not imply or require that FeSe is non-magnetic. Equating the paramagnetic state with the non-magnetic state leads to in- (c) Spin-spiral total energies E(q) per Fe of the freestanding monolayer FeSe film along high-symmetry lines, relative to the collinear stripe energy. The horizontal dotted line is the energy of the non-magnetic state. Modes with relative atomic phases ϕ=0, π, π/2 are shown. Insets: Schematics of representative antiferromagnetic configurations at high-symmetry q points for ϕ=π (green background) and π/2 (yellow background) modes.
consistencies in the calculated properties such as requiring a large band renormalization to match the ARPES data. Instead, the Hund's rule tendency of Fe to form local moments and the existence of magnetic order in related Fe-based superconducting materials, strongly indicate that inclusion of magnetic interactions and (local) moments are essential for a proper description of the ground state properties.
While ordered magnetic configurations can be calculated straightforwardly using standard DFT methods, treating the paramagnetic state -or disordered/random states in general -is more difficult: Spin-dynamics DFT calculations 26,27 can, in principle, directly simulate the time evolution of the spins but require both large supercells and long times, as well as being most applicable for configurations around an ordered state. Coherent potential approximation-type approaches 28,29 effectively proceed by considering averaged interactions. Yet another approach is to model the paramagnetic state as an explicit configurational average of a (large) number of different "snap shots" to mimic the various short-range interactions and configurations that might occur.
We address the properties of magnetism in FeSe in the spirit of this last approach and show that spin fluctuations centered around the CB-AFM configuration are unexpectedly important due to the large near degeneracy of spin wave states with energies just above the CL-AFM state; our approach to the paramagnetic state provides the first natural resolution to the ARPES-DFT paradox. The different magnetic configurations are represented by spin-spiral states of spin-wave vector q, each of which by construction has no net magnetic moment but does have a local magnetic moment. The paramagnetic state is then represented as an incoherent sum of these states. These individual planar spin-spiral states are calculated by DFT non-collinear total energy calculations using the generalized Bloch theorem. 30 In Fig. 1(a), an example of a spin-spiral is shown for a spin configuration close to the CB-AFM configuration, but with the Fe moments rotating slowly with spin wave vector q: An Fe at R has its moment pointing in a direction given by the azimuth angle φ=q·R+ϕ α , where ϕ α is an atomic phase for the α atom in the (magnetic) unit cell. Representative magnetic states at high symmetry points for two-Fe (primitive chemical cell) spin-spiral configurations are shown in insets of Fig. 1(c). While standard DFT approaches using supercells can sample these "discrete" states and compare their energies, spin-spiral calculations for varying q provide information on the connections of one representative state to another in the context of the total energy E(q) and the electronic band structure nk (q).

Results
The monolayer FeSe/STO system is modeled by a slab consisting of Se-Fe 2 -Se/TiO 2 /SrO as shown in Fig. 1(b). Although this is a minimal representation of the substrate STO, our results reasonably reproduce previous calculations that use thicker STO substrates. 14,17,18,22 The planer lattice constant a is set to the STO parameter ∼3.9Å. The vertical heights of Fe and Se above the TiO 2 are relaxed in the CB-AFM state (as summarized in Ta- ble I of the Supplementary Information), and then held fixed during the spin-spiral calculations. To isolate the interfacial effects, we consider a free-standing FeSe monolayer using the same structural parameters. In addition, we also calculate spin-spirals for bulk FeSe at the experimental lattice parameters 31 in order to compare to the inelastic neutron scattering experiment 24 and to address the possible correspondence with magnetic fluctuations in single-layer FeSe/STO.
We first consider the freestanding FeSe monolayer film. Figure 1(c) shows the spin-spiral total energy E(q) with q along high-symmetry lines in the BZ. Since there are two Fe atoms per cell, several spin-spiral "modes" appear, characterized by the relative atomic phase ϕ=ϕ 2 −ϕ 1 . A mode with ϕ=π (solid line connecting the filled circles) stably exists throughout the BZ and forms a "band" of lowest energy. This band includes the CB-AFM state at Γ, the CL-AFM state at M, and the non-collinear orthogonal state 32 at X. Around the minimum at q=M, E(q) is mostly parabolic. In contrast, E(q) around Γ (CB-AFM state) is markedly flat, i.e., there is almost no energy cost to slightly rotate the Fe moments from the CB-AFM configuration and there are many spin-wave states with (almost) the same energy as that of the CB-AFM one.
The dashed line represents the ϕ=0 mode. This mode is degenerate with the π mode at M on the zone boundary, but splits away and corresponds to a ferromagnetic configuration at Γ. This mode, however, is highly unstable in the sense that the Fe moments collapse immediately as q goes away from the zone boundary (c.f., Fig. 2 of Supplementary Information) and a self-consistent ferromagnetic configuration does not exist. The horizontal dotted line represents the energy of the non-magnetic (NM) state. When E(q) approaches or exceeds this energy, the spiral calculation of the ϕ=0 mode either converges to the NM state or never converges to a selfconsistent solution. We also find π/2 and 3π/2 modes (represented by filled squares connected by the red line) that exist only at the zone boundary and are degenerate in energy. At M, these modes give the non-collinear vortex state, 33 and at X the bicollinear stripe state.
The energy splitting of the half-integer and integer π modes is very large, amounting to ∼64 meV/Fe at M, strongly indicating that non-Heisenberg interactions are important. 31 In particular, in a simple Heisenberg model, the integer and half-integer π modes would be degenerate at M. The inclusion of the previously considered fourthorder bi-quadratic term, 31 (S i · S j ) 2 , however, is not sufficient; from extensive fits of E(q), more general 4-spin terms 34,35 are essential to account for the DFT results. This non-Heisenberg behavior has implications for the stability of possible paramagnetic states, including the nematic paramagnetic phase suggested for bulk FeSe. 25 A description of the paramagnetic state in terms of the spin-spiral states requires knowledge of E(q) not only along the high symmetry directions, but throughout the BZ. The energy landscape of the π mode for the freestanding FeSe monolayer is shown in Figs. 2(a)-(c). The region around Γ is found to be exceedingly flat: for the Γ-centered circle of radius 0.2 (2π/a) shown in Fig. 2(b), the energies are in a window of −3 to 1 meV relative to the CB-AFM energy E CB =56 meV, and the flat (yellow) region covers a substantially wide area of the 2D BZ, resulting in a sharp peak in the density of states of the π mode, Fig. 2(c). At high temperatures, entropy considerations suggest that many spin-wave states around Γ will be excited, and may be frozen in as the temperature is lowered.
The question that then arises is how the electronic structure varies for these q states near Γ compared to the CB-AFM state. To obtain the electronic energy bands in a spin-wave q state, it is necessary to carry out a k-projection procedure 36,37 (see "Method" for details), leading to an approximate "spectral weight" A q (k, ε) rather than a sharp band dispersion ε nk . The spectral weight of a paramagnetic state approximated by a q-average of spin-spiral states is shown in Fig. 2(d) together with the pure CB-AFM band structure; the qresolved evolution of the bands is given in the Supplementary Information. The bands in the pure CB-AFM state (upper panel) are consistent with previous calculations: there is an electron pocket at M and a hole-like band around Γ, which touches the Fermi level. (This band around Γ is pushed below the Fermi level, in agreement with experiment, when oxygen vacancies in the STO substrate are included. 17 ) The q-averaged spectral weight, averaged over the Γ-centered circle given in (b), shows features very similar to the pure CB-AFM bands. Although the CB-AFM spin configuration itself is easily deformed by spin wave formation, CB-AFM-like band features are very robust. These results thus provide a natural explanation for the ARPES observation that the electronic bands look like those of the CB-AFM configuration: entropy effects lead to a paramagnetic state built up from spin-wave states centered around Γ (extending over a large part of q-space), and these approximately degenerate q states have similar electronic bands, such that the spectral weight looks similar to the CB-AFM.
The objection to this scenario is that the CL-AFM phase still is lower in energy, and so the predicted electronic bands should also show CL-AFM features, in contrast to the ARPES data. The band structure of the pure CL-AFM state (the top panel of Fig. 2(e)) has hole pockets at M and Γ, in addition to two bands crossing the Fermi level along Γ-M, in good agreement with previous calculations. Taking a q-average over a circle of radius 0.2 (2π/a) centered at the BZ corner, the spectral weight is drastically changed from the pure CL-AFM bands as shown in the bottom panel of Fig. 2(e). At k=M, the convex (hole-like) and concave (electron-like) bands repel each other, with the result that states around M are pushed far away from the Fermi level. Unlike the CB-AFM case, the electronic bands of spin-spiral states with q near M are very sensitive to small changes in magnetic configuration; thus, the superposition of spin-wave states to form the paramagnetic phase results in electronic bands that do not resemble the ideal CL-AFM ones. Moreover, even if there are contributions from the ϕ=π states throughout the BZ, the electronic bands around M close to the Fermi level would be dominated by the CB-AFM-like bands since those are in the correct energy range and have large spectral weight. (Around Γ near the Fermi level, the bands are somewhat similar. There are, however, more significant differences in the electronic bands that distinguish the different magnetic configurations at energies less than about −0.4 eV.) Previous DFT calculations have shown that oxygen vacancies at the STO interface not only provide electron doping to the FeSe layer, but also modify the electronic structure in the CB-AFM state significantly. Selfdeveloped electric fields cause spin splittings of the bands and open a gap for the 3d zx/zy electron-like bands at M. 17,20 These effects have been attributed to hybridization with Se 4p z state deformed by the electric field. 20 The 3d z 2 hole-like band around Γ, which touches the Fermi level if a Coulomb correction U is not added, 10 has its energy lowered relative to the 3d zx/zy band at M and its band width reduced, leading to a disappearance of the Γ-centered hole pocket even without relying on U . 17 It is also reported that oxygen vacancies modify the magnetic interaction and reduce the CB-AFM energy E CB relative to CL-AFM. 14,17 Because of the possible importance of oxygen vacancies to the properties of FeSe/STO, we have also considered their effect on E(q). Here we use the conventional virtual crystal approximation without reducing the crystal symmetry. A single oxygen vacancy on the surface or at interface provides excess charge of nominally ∼−2e. A vacancy concentration α (0≤ α ≤1) is simulated by replacing the atomic number of interface oxygen with Z=8+2α. In Fig. 3, E(q) for the π modes of (i ) free-standing FeSe, and STO-supported FeSe (ii ) with a perfect interface (Z=8) and (iii ) with oxygen vacancies (Z=8.1) are compared. Although the flatness of E(q) around Γ remains in all cases, the CB-AFM energy E CB is significantly reduced by the presence of the STO substrate from 66 to 40 meV, and further reduced to 24 meV by the introduction of oxygen vacancies. However, the overall energy landscape still has the same shape as in Fig. 2(a)-(c), and the variation of the electronic bands with q will be similar. Thus, presence of the STO substrate -and oxygen vacancies -will increase the propensity CB-AFM-like electronic structure. For bulk FeSe we obtain a similar E(q) landscape (Supplementary Information): again the CL-AFM state is lower in energy than the CB-AFM state, and E(q) is flat around q=Γ, and is consistent with the experimentally observed coexistence of CB-and CL-AFM correlations in the bulk. 24 The similarity of E(q) for bulk and films strongly support the idea that the same coexistence is inherent in FeSe/STO as well. In addition, the increased flatness of bulk E(q) around the CL-AFM configuration (q=M) implies that CL-AFM correlations are more important in the bulk than the films while CB-AFM correlations are stronger in monolayer FeSe/STO. These quantitative differences in the bulk and film spin-spiral dispersions may be relevant to the difference in superconducting T C and to the appearance of nematic order in the 3-D bulk.

Conclusions
The calculation of spin-spiral states provides a number of insights into the magnetic and electronic properties of the monolayer FeSe, as well as providing an approach for consistently including (local) magnetic effects in the treatment of the paramagnetic state. Although the CB-AFM is not the calculated lowest energy ordered magnetic configuration, E(q) is extremely flat around the Γ point, resulting in a high density of states and entropy. These almost degenerate spin-wave states have electronic bands similar to the pure CB-AFM state (q=0). The electronic bands for states around low-energy CL-AFM configurations, on the other hand, are very sensitive to spin fluctuations, with the electronic states near M are pushed away from the Fermi level. Thus for the paramagnetic phase described in terms of spin-wave states, the resulting electronic structure around the Fermi level is expected to be dominated by CB-AFM-like features, as observed in ARPES experiments. The STO substrate, including oxygen vacancies, acts to reduce the energy of the CB-AFM (q=0) vs. CL-AFM (q=(π, π)), thus making both CB-AFM-like magnetic correlations and electronic structure more favorable. Although we do not propose any mechanism for the superconductivity, the magnetic fluctuations intrinsic to the paramagnetic state, enhanced by proximity to the STO substrate, and the resulting effect on the description of the electronic structure in the normal state are essential aspects that will need to be addressed in any comprehensive theory of the superconductivity in the FeSe system; for example, a recent effective k·p-based theory predicts 38 that the CB-AFM type fluctuations produces a fully-gapped nodeless d-wave superconducting state on the M-centered electron pockets, naturally accounting for the gap anisotropy observed in single-layer FeSe films. 39

METHOD
All DFT calculations were carried out using the fullpotential linearized augmented plane wave method 40 as implemented in the HiLAPW code, including noncollinear magnetism 41 and generalized-Bloch-theorem spin-spiral magnetism. Orientation of the spin density is fully relaxed throughout the entire space. 42 The muffintin sphere radius was set to 0.8Å for oxygen and 1.1Å for other atoms, and the wave function and density and potential cutoffs were 16 and 200 Ry, respectively. The Perdew, Burke, and Ernzerhof form of the Generalized Gradient Approximation was used for exchange correlation. The Brillouin zone was sampled with 20×20×1 and 20×20×15 k-point meshes for the film and bulk calculations, respectively. The density of two-dimensional spin-wave states (Fig. 2(c)) was calculated by the triangle method. 43 In the spin-spiral calculations, the cell-periodic partρ of the spin off-diagonal density matrix ρ ↑↓ for wave vector q has the formρ = e iq·r ρ ↑↓ in the interstitial region, and in the sphere regionρ = e iq·τ ρ ↑↓ with atomic position τ .
To obtain the band dispersion for a spin-wave state q, a k-projection in the same spirit as the unfolding technique for supercell calculations is needed. 36,37,44 The wave function of a band ε nk (q) is a two component spinor, where each spin component has a different Bloch phase, ψ q nk = e i(k−q/2)·r u q ↑ nk e i(k+q/2)·r u q ↓ nk . (1) Thus this band must be projected onto two different k vectors, to obtain a band structure consistent with standard (supercell) calculations and with the momentum-resolved measurement in ARPES experiments. The projection weights for k ± are w + = u ↑ k |u ↑ k and w − = u ↓ k |u ↓ k . They can be found from the expectation value of σ z by w ± = 1 ± ψ k |σ z |ψ k 2 (3) since ψ k |σ z |ψ k = w + −w − and the orthonormality condition gives w + + w − = 1. (The band index and q have been omitted for simplicity.) The integrals implicit in the brackets are performed over the primitive cell. Figure 2(f) gives an example of how to obtain the dispersion along the Γ-M line: The band energies are calculated along two k-rods shifted from the Γ-M segment by ±q/2, and then projected onto Γ-M. To get the q-averaged spectral weights shown in Fig. 2(d)(e), the two paths shown in (g) are averaged since they become inequivalent for arbitrary q.