Magnetic fluctuations in single-layer FeSe

The electronic structure of single-layer FeSe films on SrTiO3 presents a quandary: experimentally there is no long-range magnetic order, but the observed bands are reasonably well described by density functional calculations assuming the checkerboard antiferromagnetic (CB-AFM) ordering despite this configuration not being the calculated ground state. Here we investigate the paramagnetic nature of this system via first-principles spin-spiral calculations. Fits of the spin-spiral dispersion to spin models place this S = 1 spin system in a region of parameter space where CB-AFM quantum fluctuations lead to a magnetically disordered paramagnetic state. Modeling the paramagnetic state as an incoherent superposition of spin-spiral states arising from thermal and/or quantum fluctuations, the resulting electronic bands around the Fermi level are found to closely resemble those of the ordered CB-AFM configuration, thus providing a consistent explanation of the angle-resolved photoemission observations. These results suggest that CB-AFM fluctuations play a more important role than previously thought. The paper reports a theoretical investigation of the paramagnetic states of monolayer FeSe. The authors use first-principles spin-spiral calculations to show that spin fluctuations around the checkerboard states explain the paradox between angle-resolved photoemission experiments and density-functional theory.

S ingle-layer FeSe films grown on SrTiO 3 (001) (STO) have generated intense interest because of their reported highsuperconducting critical temperature T C~4 0-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: single-layer FeSe/ STO films as prepared 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 (BZ) 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 Mcentered pockets 11,12 through (π,π) collinear stripe antiferromagnetic (CL-AFM) spin fluctuations. Oxygen vacancies (Ovac) 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, the observation that bilayer FeSe films are not superconducting 7 suggests that the FeSe-STO interface (Fig. 1a) 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 are consistent with the experimental data: both the Fermi surface and the band structure throughout the whole 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 O-vac 13-17. O-vac are found to electron-dope the FeSe layer, and to modify the band structure of the CB-AFM state, including shifting the hole-like band around Γ down 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 "long-range" 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 of "paramagnetic" FeSe appear to agree with the "ordered" CB-AFM 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 (2 Fe) unit cell. These conditions, however, do not imply or require that FeSe is nonmagnetic. Equating the paramagnetic state with the non-magnetic state leads to inconsistencies 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 Schematic of a planer spin-spiral with wave vector q = (1/36,0) r.l.u. and relative atomic phase φ = π. The shaded square represents the crystallographic primitive cell. c Spin-spiral total energies E(q) per Fe of the free-standing 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 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 electronic properties related to the magnetism of 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 spinspiral states with energies just above the CL-AFM state and to intrinsic quantum fluctuations arising from magnetic frustration; our approach to the paramagnetic state provides a natural resolution to the ARPES-DFT paradox.

Results
Spin-spirals and ordered magnetic states. The relevant independent degrees of freedom when considering magnetic fluctuations are the individual Fe spins. The resulting different magnetic configurations are represented (via the standard Fourier transformation) by normalized spin-spiral states of wave vector q. Each spin-spiral state by construction has no net magnetic moment but does have a local magnetic moment. The paramagnetic state is then approximated 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. 1b, an example of a spinspiral is shown for a spin configuration close to the CB-AFM configuration, but with the Fe moments rotating slowly with wave vector q: an Fe at position 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 as insets in Fig. 1c. 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).
The monolayer FeSe/STO system is modeled by a slab consisting of Se-Fe 2 -Se/TiO 2 /SrO as shown in Fig. 1a. Although this is a minimal representation of the substrate STO, our results    Figure 1). 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 singlelayer FeSe/STO. We first consider the free-standing FeSe monolayer film. Figure 1c 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-spiral 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 (Supplementary Figure 2) 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 self-consistent 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 per Fe at M, strongly indicating that non-Heisenberg interactions are important 31 . In particular, in a simple Heisenberg model, the integer and halfinteger π modes would be degenerate at M. The inclusion of the previously considered fourth-order bi-quadratic term 31 , (S i ⋅ S j ) 2 , however, is not sufficient; from extensive fits of E(q), more general four-spin terms 34,35 are essential to account for the DFT results.
The fits including (not including) fourth-order terms give 35 Heisenberg parameters J 2 /J 1 of 0.59 (0.67) for the free FeSe monolayer, which reduce to 0.55 (0.63) for ideal FeSe/STO, and farther to 0.53 (0.57) for FeSe/STO with O-vac, in large part because of the reduction in the CB-CL energy difference (Supplementary Table 2). For the classical spin case (S → ∞), the CL-AFM is stable for J 2 =J 1 > 1 2 , consistent with the calculated DFT total energies for the two ordered phases. However, this system is magnetically frustrated, with the result that Néel/CB-AFM quantum fluctuations in this S = 1 system 36-39 are important for J 2 /J 1 up to around 0.7, leading to a disordered paramagnetic phase. Inclusion of bi-quadratic terms in the spin models 40 similarly favors Néel fluctuations beyond the classical stability limit for the CB-AFM state.
The connection between the DFT spin-spiral calculationswhich do not include quantum fluctuations-and the (extended) S = 1 Heisenberg model has several implications: (i) The result that the DFT total energy of the CL-AFM configuration is lower than the CB-AFM (J 2 =J 1 > 1 2 ) is necessary for the observed paramagnetism; if the CB-AFM phase were lower or degenerate in energy with that of the CL-AFM, an ordered Néel phase would be expected. On the other hand, if the energy differences between the CL-AFM and CB-AFM configurations were much larger, the increased J 2 /J 1 values then would be in a part of phase space where the CL-AFM phase would be stable and the Néel fluctuations would no longer be important. (ii) CB-AFM-like quantum fluctuations corresponding to q states around Γ will be important at all temperatures, including T = 0. In addition, thermal effects may also introduce CB-AFM contributions. (iii) The fourth-order non-Heisenberg contributions have implications for the types 35 and 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 spinspiral states requires knowledge of E(q) not only along the highsymmetry directions, but throughout the BZ. The energy landscape of the π mode for the free-standing FeSe monolayer is shown in Fig. 2a-c. The region around Γ is found to be exceedingly flat: for the Γ-centered circle of radius 0.2 (2π/a) shown in Fig. 2b, the energies fall 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 twodimensional BZ, resulting in a sharp peak in the density of states of the π mode, Fig. 2c. At high temperatures, entropy considerations suggest that many spin-spiral states around Γ will be excited, and may be frozen in as the temperature is lowered; this contribution is in addition to the quantum fluctuations discussed above.
Electronic spectra of paramagnetic phase. Given that our DFT calculations (and experiment) imply the existence of Néel-like quantum fluctuations, the question arises of how the electronic structure varies for these q states near Γ compared to that of the CB-AFM state (Fig. 2d). To obtain the electronic energy bands in a spin-spiral q state, it is necessary to carry out a k-projection procedure 41-43 (see Methods section for details), leading to an approximate "spectral weight" A q (k,ε) rather than a sharp band dispersion ε nk . To investigate the effect of magnetic fluctuations on the electronic bands, the spectral weight corresponding to a qaverage of spin-spiral states around q = 0 is shown in Fig. 2e; the q-resolved evolution of the bands is given in the Supplementary Figure 3. The bands in the pure CB-AFM state (Fig. 2d) 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 O-vac in the STO substrate are included 17 . The q-averaged spectral weight (Fig. 2e), averaged over the Γ-centered circle given in Fig. 2b, shows features very similar to the pure CB-AFM bands. Although the CB-AFM spin configuration itself is easily deformed by spin-spiral formation, CB-AFM-like electronic band features are very robust and extend over a large portion of q-space. Thus, because of the quantum (and possibly thermal) fluctuations, the electronic spectral weight should exhibit CB-AFM-like features.
In addition to the magnetic fluctuations around q = 0, the fluctuations around the CL-AFM state need to be considered since the CL-AFM phase is the calculated ground state and the nearby spin-spiral states will contribute to the paramagnetic state. Even in the presence of CB-AFM quantum fluctuations, the electronic bands would be expected to show CL-AFM-like features, but these are not seen in the ARPES data. The band structure of the pure CL-AFM state (Fig. 2f) has hole pockets at M and Γ, in addition to two bands crossing the Fermi level along Γ-M, in good agreement with previous calculations. Accounting for fluctuations around the CL-AFM state by 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 Fig. 2g. 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 contributions from spin-spiral states around M 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. To demonstrate that this is in fact the case, Fig. 3 and Supplementary Figure 4 present the electronic bands arising from a simple Boltzmann thermal average, 〈A q (k,ε)exp(−E(q)/k B T)〉, at different temperatures for the free monolayer assuming that the (normalized) spin-spiral states are the independent degrees of freedom; the CB-AFM-like features around M are clearly visible for T = 0.3 (Fig. 3a) even neglecting the quantum fluctuations that will enhance these features at all temperatures, even as T → 0, and become more pronounced as the temperature increases (Fig. 3b). For low temperatures the calculated thermally averaged spectra alone clearly will not agree with experiment since they neglect the Néel quantum fluctuations contributions (Fig. 2e). These thermal-only spectra, however, do demonstrate that the CB-like electronic structure already becomes dominant at an energy scale a fraction of the CB-CL energy difference.
These results provide a natural explanation for the ARPES observation that the electronic bands look like those of the CB-AFM configuration: because of both entropic (thermal) effects and quantum fluctuations, the paramagnetic phase will have contributions from q states centered around Γ (but extending over a large part of q-space), states that all have similar electronic bands; the spectral weight coming from spin-spiral states around q = M, on the other hand, does not have the characteristic CL-AFM features nor is it in the spectral range typically probed in the ARPES experiments.
Previous DFT calculations have shown that O-vac 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. Self-developed 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 orbitals deformed by the electric field 20 . Spin-orbit effects can also open gaps at M even without O-vac 17 (Supplementary  Tables 3 and 4). 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 the disappearance of the Γ-centered hole pocket even without relying on U 17 . It is also reported that O-vac modify the magnetic interaction and reduce the CB-AFM energy E CB relative to CL-AFM 14,17 . Because of the possible importance of O-vac 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 approximately −2e. A vacancy concentration α (0 ≤ α ≤ 1) is simulated by replacing the atomic number of interface oxygen with Z = 8 + 2α. In Fig. 4, E(q) for the π and π/2 modes of (i) free-standing FeSe, and STOsupported FeSe (ii) with a perfect interface (Z = 8), and (iii) with O-vac (Z = 8.1) are compared; the corresponding CB-AFM bands are given in Supplementary Figure 5. 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 56 to 40 meV, and further reduced to 24 meV by the introduction of O-vac. However, the overall energy landscape still has the same shape as in Fig. 2a-c, and the variation of the electronic bands with q will be similar. Thus, the presence of the STO substrateand oxygen vacancies O-vac-will increase the propensity CB-AFM-like electronic structure. For bulk FeSe we obtain a similar E(q) landscape (Supplementary Figure 6): 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 CBand 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 three-dimensional bulk.

Discussion
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, the magnetic system is in a region of phase space where Néel (CB-AFM) quantum fluctuations are important and lead to a disordered paramagnetic phase as observed experimentally. The calculated E(q) is extremely flat around the Γ point, resulting in a high density of spin-spiral states and entropy. These almost degenerate spin-spiral 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 pushed away from the Fermi level. Thus for the paramagnetic phase described in terms of spinspiral 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 O-vac, 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 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 44 that the CB-AFM type fluctuations produce 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 45 .

Methods
All DFT calculations were carried out using the full-potential linearized augmented plane wave method 46 as implemented in the HiLAPW code, including noncollinear magnetism 47 and generalized Bloch-theorem spin-spiral magnetism. Orientation of the spin density is fully relaxed throughout the entire space 48 . The muffin-tin 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 49 of the generalized gradient approximation was used for exchange correlation. The BZ 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-spiral states (Fig. 2c) was calculated by the triangle method 50 .
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 electronic band dispersion for a spin-spiral state q, a k-projection in the same spirit as the unfolding technique for supercell calculations is needed [41][42][43] . The wave function of a band ε nk (q) is a two component spinor, where each spin component has a different Bloch phase, 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 þ ¼ hu " k ju " k i and w À ¼ hu # k ju # k i, and can be found from the expectation value of σ z by since 〈ψ k |σ z |ψ k 〉 = w + − w − and the orthonormality condition gives w + + w − = 1. The band index n and q have been omitted for simplicity. The integrals implicit in the brackets are performed over the primitive cell. Figure 5a 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. The electronic spectral weights along Γ-X-M-Γ line are calculated along the two k paths shown in Fig. 5b and added since these paths are inequivalent for a general q state. The triangle method is employed in generating the q-averaged spectral weights in Figs. 2d-g, 3. The self-consistent results on the coarse q-mesh points (20 × 20)total energy, band energy, Fermi energy, spectral weight, and so on-are linearly interpolated, and the q summation is done on a dense mesh (1000 × 1000).
Data availability. The data supporting the findings of this study are available within the paper, its Supplementary Information files, and from the corresponding author upon reasonable request. The two Γ-X-M-Γ paths of k used to calculate the electronic spectral weights