Frustrated spin order and stripe fluctuations in FeSe

Abstract

The charge and spin dynamics of the structurally simplest iron-based superconductor, FeSe, may hold the key to understanding the physics of high temperature superconductors in general. Unlike the iron pnictides, FeSe lacks long range magnetic order in spite of a similar structural transition around 90 K. Here, we report results of Raman scattering experiments as a function of temperature and polarization and simulations based on exact diagonalization of a frustrated spin model. Both experiment and theory find a persistent low energy peak close to 500 cm⁻¹ in B_1g symmetry, which softens slightly around 100 K, that we assign to spin excitations. By comparing with results from neutron scattering, this study provides evidence for nearly frustrated stripe order in FeSe.

Introduction

Fe-based pnictides and chalcogenides, similar to cuprates, manganites or some heavy fermion compounds, are characterized by the proximity and competition of various phases including magnetism, charge order, and superconductivity. Specifically the magnetism of Fe-based systems has various puzzling aspects which do not straightforwardly follow from the Fe valence or changes in the Fermi surface topology^1,2,3,4. Some systems have a nearly ordered localized moment close to 2μ_B⁵, such as FeTe or rare-earth iron selenides, whereas the moments of AFe₂As₂-based compounds (A = Ba, Sr, Eu or Ca) are slightly below 1μ_B⁶ and display aspects of itinerant spin-density-wave (SDW) magnetism with a gap in the electronic excitation spectrum⁷. In contrast others do not order down to the lowest temperatures, such as FeSe⁸ or LaFePO⁹.

The material specific differences are a matter of intense discussion, and low- as well as high-energy electronic and structural properties determine the properties^{1,2,4,10,11,12,13}. At the Fermi energy E_F, the main fraction of the electronic density of states originates from t_2g Fe orbitals, but a substantial part of the Fe–Fe hopping occurs via the pnictogen or chalcogen atoms, hence via the xz, yz, and p_z orbitals. For geometrical reasons, the resulting exchange coupling energies between nearest (J₁) and next nearest neighbor (J₂) iron atoms have the same order of magnitude, and small changes in the pnictogen (chalcogen) height above the Fe plane influence the ratio J₂/J₁, such that various orders are energetically very close¹².

The reduced overlap of the in-plane xy orbitals decreases the hopping integral t and increases the influence of the Hund’s rule interactions and the correlation energy U, even though they are only in the range of 1–2 eV. Thus the electrons in the xy orbitals have a considerably higher effective mass m* and smaller quasiparticle weight Z than those of the xz and yz orbitals. This effect was coined orbital selective Mottness^14,15,16 and was observed by photoemission spectroscopy (ARPES) in Fe-based chalcogenides¹⁷. It is similar in spirit to what was found by Raman scattering in the cuprates as a function of momentum¹⁸. In either case some of the electron wave functions are more localized than others. This paradigm may explain why the description remains difficult and controversial in all cases.

Therefore we address the question as to whether systematic trends can be found across the families of the Fe-based superconductors, how the spin excitations are related to other highly correlated systems, and how they can be described appropriately.

As an experimental tool we use Raman scattering since the differences expected theoretically^1,3 and indicated experimentally in the electronic structure⁷ can be tracked in both the charge and the spin channel. Another advantage is the large energy range of approximately 1 meV to 1 eV (8–8000 cm⁻¹) accessible by light scattering¹⁹.

Early theoretical work on Fe-based systems considered the Heisenberg model the most appropriate approach²⁰, and the high-energy maxima observed by Raman scattering in BaFe₂As₂ were interpreted in terms of localized spins^21,22. On the other hand, the low-energy spectra are reminiscent of charge density wave (CDW) or SDW formation^22,23,24,25. In principle, both effects can coexist if the strength of the correlations varies for electrons from different orbitals, where itinerant electrons form a SDW, while those on localized orbitals give rise to a Heisenberg-like response.

In contrast to the AFe₂As₂-based compounds, FeSe seems to be closer to localized order with a larger mass renormalization than in the iron pnictides¹. Apart from low lying charge excitations, the remaining, presumably spin, degrees of freedom in FeSe may be adequately described by a spin-1 J₁-J₂-J₃-K Heisenberg model¹² which provides also a consistent description of our results shown in this work and allows for the presence of different spin orders. Since various types of spin order are energetically in close proximity^12,26,27, frustration may quench long-range order down to the lowest temperatures⁸, even though neutron scattering experiments in FeSe find large values for the exchange energies^27,28.

Recent experiments on FeSe focused on low-energies and B_1g (x² − y²) symmetry, and the response was associated with particle-hole excitations and critical fluctuations²⁹. Here, we obtain similar experimental results below 1500 cm⁻¹. Those in the range 50–200 cm⁻¹ show similarities with the other Fe-based systems while those above 200 cm⁻¹ are distinctly different but display similarities with the cuprates^30,31. In addition to previous work, we analyze all symmetries at higher energies up to 3500 cm⁻¹, to uncover crucial information about the behavior of the spin degrees of freedom.

By comparing experimental and simulated Raman data we find a persistent low-energy peak at roughly 500 cm⁻¹ in B_1g symmetry, which softens slightly around 100 K. We assign the B_1g maximum and the related structures in A_1g and B_2g symmetry to spin excitations. The theoretical simulations also aim at establishing a link between light and neutron scattering data with respect to the spin degrees of freedom and to furnish evidence for nearly frustrated stripe order at low temperature. We arrive at the conclusion that frustrated order of localized spins dominates the physics in FeSe, while critical spin and/or charge fluctuations are not the main focus of the paper.

Results

Experiments

Symmetry-resolved Raman spectra of single-crystalline FeSe (see Methods) in the energy range up to 0.45 eV (3600 cm⁻¹) are shown in Fig. 1. The spectra are linear combinations of the polarization dependent raw data (see Methods and Supplementary Fig. 1 in Supplementary Note 1). For B_1g symmetry (Fig. 1a) we plot only two temperatures, 40 and 300 K, to highlight the persistence of the peak at ~500 cm⁻¹. The full temperature dependence will be shown below. For A_1g, A_2g, and B_2g symmetry we show spectra at 40, 90, and 300 K (Fig. 1b–d). Out of the four symmetries, the A_1g, B_1g, and B_2g spectra display Raman active phonons, magnons or electron-hole excitations, while the A_2g spectra are weak and vanish below 500–1000 cm⁻¹. As intensity in A_2g symmetry appears only under certain conditions not satisfied in the present study, we ignore it here.

Fig. 1

Symmetry-resolved Raman spectra of FeSe at various temperatures for large energy transfers. a B_1g spectra at temperatures as indicated. The spectrum at 90 K is omitted here for clarity but is displayed in a separate figure below The weak structure at T = 40 K in the range 20–25 cm⁻¹ is left over from the fluctuation peak which is most pronounced right above T_s as shown below. b A_1g, c A_2g, and d B_2g spectra at temperatures as indicated. In A_1g and B_2g symmetry particle-hole excitations dominate the response. In agreement with the simulations weak additional peaks from spin excitations appear at low temperature (blue shaded areas). B_2g shows a loss of spectral weight (shaded red). The narrow lines close to 200 cm⁻¹ are the A_1g and B_1g phonons. In the 1 Fe unit cell used here the B_1g phonon appears in B_2g symmetry since the axes are rotated by 45° with respect to the crystallographic (2 Fe) cell. The A_2g intensity vanishes below 500 cm⁻¹ and the cross section is completely temperature independent

In the high-energy limit the intensities are smaller in all symmetries than those in other Fe-based systems such as BaFe₂As₂ (see Supplementary Fig. 2 in Supplementary Note 2). However, in the energy range up to ~3000 cm⁻¹ there is a huge additional contribution to the B_1g cross section in FeSe (Fig. 1a). The response is strongly temperature dependent and peaks at 530 cm⁻¹ in the low-temperature limit. Between 90 and 40 K the A_1g and B_2g spectra increase slightly in the range around 700 and 3000 cm⁻¹, respectively (indicated as blue shaded areas in Fig. 1b, d). The overall intensity gain in the A_1g and B_2g spectra in the shaded range is a fraction of ~5% of that in B_1g symmetry. The B_2g spectra exhibit a reduction in spectral weight in the range from 600 to 1900 cm⁻¹ (shaded red) which is already fully developed at the structural transition at T_s = 89.1 K in agreement with earlier work²⁹. In contrast to A_1g and B_2g symmetry, the temperature dependence of the B_1g intensity is strong, whereas the peak energy changes only weakly, displaying some similarity with the cuprates³². This similarity, along with the considerations of Glasbrenner et al.¹², motivated us to explore a spin-only, Heisenberg-like model for describing the temperature evolution of the Raman scattering data.

Simulations at zero temperature

We performed numerical simulations at zero temperature for a frustrated spin-1 system on the basis of a J₁-J₂-J₃-K Heisenberg model¹² on a 16-points cluster as shown in Fig. 2a and described in the Methods section. Figure 2b shows the resulting phase diagram as a function of J₂ and J₃. K was set at 0.1 (repulsive) in order to suppress ordering tendencies on the small cluster. The parameter set for the simulations of the Raman and neutron data at finite temperature is indicated as a black dot.

Fig. 2

Model and resulting phase diagram. a A 4 × 4 cluster was used for the simulations. The red spheres represent the Fe atoms, each of which carries localized spin S_i, with S = 1. The nearest, next-nearest, and next-next-nearest neighbor interactions J₁, J₂, and J₃, respectively, are indicated. K is the coefficient of the biquadratic term proportional to (S_i · S_j)². b J₂ − J₃ phase diagram as obtained from our simulations at T = 0 and for K = 0.1. The black dot shows the parameters at which temperature-dependent simulations have been performed

In Fig. 3 we show the low-temperature data (Fig. 3a) along with the simulations (Fig. 3b). The energy scale for the simulations is given in units of J₁ which has been derived¹² to be 123 meV or 990 cm⁻¹, allowing a semi-quantitative comparison with the experiment. As already mentioned, the experimental A_1g and B_2g spectra are not dominated by spin excitations and we do not attempt to further analyze the continua extending to energies in excess of 1 eV, considering them a background. The opposite is true for B_1g symmetry, also borne out in the simulations. For the selected values of J₁ = 123 meV, J₂ = 0.528J₁, J₃ = 0, and K = 0.1J₁, the positions of the spin excitations in the three symmetries and the relative intensities are qualitatively reproduced. The choice of parameters is motivated by the previous use of the J₁-J₂ Heisenberg model, with J₁ = J₂ to describe the stripe phase of iron pnictides²⁰. Here we use a value of J₂ smaller than J₁ to enhance competition between Néel and stripe orders when describing FeSe. This approach and choice of parameters is strongly supported in a recent neutron scattering study²⁷.

Fig. 3

Symmetry-resolved Raman spectra of FeSe for large energy shifts at low temperature. a Experimental results for symmetries as indicated at 40 K. The B_1g peak at 500 cm⁻¹ dominates the spectrum. In A_1g and B_2g symmetry the electron-hole continua dominate the response, and the magnetic excitations yield only small additional contributions at approximately 700 and 3000 cm⁻¹, respectively. b Simulated Raman spectra at T = 0 K including only magnetic contributions. The A_1g and B_1g symmetries have peaks solely at low energies whereas the B_2g contributions are at high energies only. The B_1g response is multiplied by a factor of 0.02

The comparison of the different scattering symmetries, the temperature dependence, and our simulations indicate that the excitation at 500 cm⁻¹ is an additional scattering channel superimposed on the particle-hole continuum and fluctuation response, as shown in Supplementary Note 3 with Supplementary Figs. 3 and 4. Here we focus on the peak centered at ~500 cm⁻¹ which, in agreement with the simulations, originates from two-magnon excitations in a highly frustrated spin system, although the features below 500 cm⁻¹ also are interesting and were interpreted in terms of quadrupolar orbital fluctuations²⁹.

Temperature dependence

It is enlightening to look at the B_1g spectra across the whole temperature range as plotted in Fig. 4. The well-defined two-magnon peak centered at ~500 cm⁻¹ in the low temperature limit loses intensity, and becomes less well-defined with increasing temperature up to the structural transition T_s = 89.1 K. Above the structural transition, the spectral weight continues to decrease and the width of the two-magnon feature grows, while the peak again becomes well-defined and the energy increases slightly approaching the high temperature limit of the study. What may appear as a gap opening at low temperature is presumably just the reduction of spectral weight in a low-energy feature at ~22 cm⁻¹. The intensity of this lower energy response increases with temperature, leading to a well-formed peak at an energy around 50 cm⁻¹ near the structural transition. Above the structural transition this feature rapidly loses spectral weight, hardens, and becomes indistinguishable from the two-magnon response in the high temperature limit. This low-energy feature develops in a fashion very similar to that found in Ba(Fe_1−xCo_x)₂As₂ for x > 0^33,34,35.

Fig. 4

Raman spectra in B_1g symmetry at temperatures a below and b above T_s = 89.1 K. The spectrum at 91 K appears in both panels for better comparison. The black dashed line in a and the gray shaded area in b indicate the approximate positions of the low-energy peak resulting from critical fluctuations. The peak centered close to 500 cm⁻¹ results from excitations of neighboring spins which are studied here. A tentative decomposition is shown in Supplementary Fig. 4

Now we compare the measurements with numerical simulations for the temperature dependence of the Raman B_1g susceptibility in Fig. 5a, b, respectively. For the simulations (Fig. 5b) we use the same parameters as at T = 0 (black dot in Fig. 3). At zero temperature the simulations show a single low energy B_1g peak around 0.3J₁. As temperature increases, a weak shoulder forms on the low-energy side of the peak, and the whole peak softens slightly and broadens over the simulated temperature range. Except for the additional intensity at low energies, Ω < 200 cm⁻¹ (Fig. 5a), there is good qualitative agreement between theory and experiment. As shown in Supplementary Fig. 5 in Supplementary Note 4, a similar agreement between experiment and simulations is obtained for the temperature dependence in A_1g and B_2g symmetries, indicating that both the gain in intensity (blue shaded areas in Fig. 1) as well as the reduction in spectral weight in B_2g from 600 to 1900 cm⁻¹ (shaded red in Fig. 1d) can be attributed to the frustrated localized magnetism.

Fig. 5

Temperature dependence of the B_1g response. a Experimental spectra at selected temperatures as indicated. The spectra include several excitations the decomposition of which is shown in Supplementary Fig. 4. b Simulated Raman response at temperatures as indicated. Only magnetic excitations are included. The coupling constant was derived as J₁ = 123 meV in ref. ¹²

Connection to the spin structure factor

To support our explanation of the Raman data, we simulated the dynamical spin structure factor S(q, ω) and compared the findings to results of neutron scattering experiments²⁷. While clearly not observing long-range order, above the structural transition neutron scattering finds similar intensity at finite energy for several wave vectors along the line (π, 0) − (π, π). Upon cooling, the spectral weight at these wave vectors shifts away from (π, π) to directions along (π, 0), although the respective peaks remain relatively broad. In Fig. 6a, b we show the results of the simulations for two characteristic temperatures. As the temperature decreases, spectral weight shifts from (π, π) toward (π, 0) in agreement with the experiment²⁷. In Fig. 6c we show the evolution of the spectral weights around (π, π) and (π, 0) in an energy window of (0.4 ± 0.1)J₁ as a function of temperature, similar to the results shown in ref. ²⁷. In the experiment, the temperature where the integrated dynamical spin structure factor changes most dramatically is close to the structural transition. From our simulations, the temperature where similar changes occur in comparison to neutron scattering corresponds to the temperature at which the simulated B_1g response (Fig. 5) shows the most pronounced shoulder, and the overall intensity begins to decrease. Not surprisingly, the low-energy peak in the Raman scattering experiment is also strongest near the structural transition.

Fig. 6

Simulations of the dynamical structure factor S(q, ω) of localized spin excitations integrated over an energy window of (0.4 ± 0.1)J₁. a, b display cuts through the first Brillouin zone at T = 0.25 and 0J₁, respectively. At high temperature there is intensity at (π, π) indicating a tendency toward Néel order. At low temperature the intensity at (π, π) is reduced and the stripe-like antiferromagnetism with (π, 0) ordering wave vector becomes stronger. c S(q, ω) integrated over an energy window (0.4 ± 0.1)J₁ for fixed momenta (π, π) and (π, 0) intensities as a function of temperature

Discussion

The agreement of experiment with theory in both neutron and Raman scattering suggests that a dominant contribution to the FeSe spectra comes from frustrated magnetism of essentially local spins. The differences between the classes of ferro-pnictides and -chalcogenides, in particular the different degrees of itineracy, may then originate in a subtle orbital differentiation across the families¹.

If FeSe were frustrated, near such a phase boundary between magnetic states, then its behavior would be consistent with the observed sensitivity to intercalation^36,37, layer thickness³⁸, and pressure³⁹, which could affect the exchange interactions through the hopping. Relative to the theoretical results below 200 cm⁻¹, critical fluctuations of any origin, which are characterized by a diverging correlation length close to the transition, can neither be described nor distiguished in such a small cluster calculation. Here, only experimental arguments can be applied similar to those in ref. ³⁵, but will not be further discussed, since they are not the primary focus of the analysis. A brief summary may be found in the Supplementary Note 3.

It is remarkable how clearly the Raman spectra of an SDW state originating from a Fermi surface instability and a magnet with local moments can be distingiuished. For comparison, Fig. 7 shows B_1g Raman spectra for La₂CuO₄ and BaFe₂As₂ at characteristic temperatures. La₂CuO₄ (Fig. 7a) is an example of a material with local moments on the Cu sites^30,31 having a Néel temperature of T_N = 325 K. The well-defined peak at ~2.84J₁^40,41 possesses a weak and continuous temperature dependence across T_N³². The origin of the scattering in La₂CuO₄ and other insulating cuprates⁴² can thus be traced back to Heisenberg-type physics of local moments⁴³, which, for simplicity, need only include the nearest-neighbor exchange interaction J₁.

Fig. 7

Examples of localized and itinerant magnets. a B_1g Raman spectra of La₂CuO₄. From ref. ³¹. b B_1g spectra of BaFe₂As₂ at four characteristic temperatures as indicated

In contrast, most iron-based superconductors are metallic antiferromagnets in the parent state exhibiting rather different Raman signatures. In BaFe₂As₂ (Fig. 7b) abrupt changes are observed in B_1g symmetry upon entering the SDW state: the fluctuation peak below 100 cm⁻¹ vanishes, a gap develops below some 500–600 cm⁻¹, and intensity piles up in the range 600–1500 cm⁻¹ (ref. ^22,44), the typical behavior of an SDW or CDW²⁴ in weak-coupling, resulting from Fermi surface nesting. Yet, even for itinerant systems such as these, longer range exchange interactions can become relevant and lead to magnetic frustration⁴⁵.

In summary, the Raman response of FeSe was measured in all symmetries and compared to simulations of a frustrated spin-1 system. The experimental data were decomposed in order to determine which parts of the spectra originate from particle-hole excitations, fluctuations of local spins, and low energy critical fluctuations. Comparison of the decomposed experimental data with the simulations gives evidence that the dominant contribution of the Raman spectra comes from magnetic competition between (π, 0) and (π, π) ordering vectors. These features of the Raman spectra, which agree qualitatively with a spin-only model, consist of a dominant peak in B_1g symmetry around 500 cm⁻¹ along with a peak at similar energy but lower intensity in A_1g and at higher energy in B_2g symmetry. These results will likely help to unravel the mechanism behind the superconducting phase found in FeSe.

Methods

Experiment

The FeSe crystals were prepared by the vapor transport technique. Details of the crystal growth and characterization are described elsewhere⁴⁶. Before the experiment the samples were cleaved in air and the exposure time was minimized. The surfaces obtained in this way have several atomically flat regions allowing us to measure spectra down to 5 cm⁻¹. At the tetragonal-to-orthorhombic transition T_s twin boundaries appear and become clearly visible in the observation optics. As described in detail by Kretzschmar et al.³⁵ the appearance of stripes can be used to determine the laser heating ΔT_L and T_s to be (0.5 ± 0.1) Km W⁻¹ and (89.1 ± 0.2) K, respectively.

Calibrated Raman scattering equipment was used for the experiment. The samples were attached to the cold finger of a He-flow cryostat having a vacuum of ~5 × 10⁻⁵ Pa (5 × 10⁻⁷ mbar). For excitation we used a diode-pumped solid state laser emitting at 575 nm (Coherent GENESIS MX-SLM 577-500) and various lines of an Ar ion laser (Coherent Innova 304). The angle of incidence was close to 66° for reducing the elastic stray light entering the collection optics. Polarization and power of the incoming light were adjusted in a way that the light inside the sample had the proper polarization state and, respectively, a power of typically P_a = 4 mW independent of polarization. For the symmetry assignment we use the 1 Fe unit cell (axes x and y parallel to the Fe–Fe bonds) which has the same orientation as the magnetic unit cell in the cases of Néel or single-stripe order (4 Fe cell). The orthorhombic distortion is along these axes whereas the crystallographic cell assumes a diamond shape with the length of the tetragonal axes preserved. Because of the rotated axes in the 1 Fe unit cell the Fe B_1g phonon appears in the B_2g spectra. Spectra at low to medium energies were measured with a resolution σ ≈ 5 cm⁻¹ in steps of ΔΩ = 2.5 or 5 cm⁻¹ below 250 cm⁻¹ and steps of 10 cm⁻¹ above where no sharp peaks need to be resolved. Spectra covering the energy range up to 0.5–1 eV were measured with a resolution σ ≈ 20 cm⁻¹ in steps of ΔΩ = 50 cm⁻¹.

Simulations

We use exact diagonalization to study a Heisenberg-like model on a 16-site square lattice, which contains the necessary momentum points and is small enough that exact diagonalization can reach high enough temperatures to find agreement with the temperature dependence in the experiment. This was solved using the parallel Arnoldi method⁴⁷. The Hamiltonian is given by

$$\begin{array}{*{20}{l}} {\cal H} \hfill & = \hfill & {\mathop {\sum}\limits_{{\mathrm{nn}}} \left[ {J_1{\bf{S}}_i \cdot {\bf{S}}_j + K\left( {{\bf{S}}_i \cdot {\bf{S}}_j} \right)^2} \right]} \hfill \\ {} \hfill & {} \hfill & { + \mathop {\sum}\limits_{{\mathrm{2nn}}} {\kern 1pt} J_2{\bf{S}}_i \cdot {\bf{S}}_j + \mathop {\sum}\limits_{{\mathrm{3nn}}} {\kern 1pt} J_3{\bf{S}}_i \cdot {\bf{S}}_j} \hfill \end{array}$$

(1)

where S_i is a spin-1 operator reflecting the observation that the local moments of iron chalcogenides are close to 2μ_B⁴⁸. The sum over nn is over nearest neighbors, the sum over 2nn is over next nearest neighbors, and the sum over 3nn is over next next nearest neighbors.

We determine the dominant order according to the largest static spin structure factor, given by

$$S({\bf{q}}) = \frac{1}{N}\mathop {\sum}\limits_l {\kern 1pt} {\mathrm{e}}^{i{\bf{q}} \cdot {\bf{R}}_l}\mathop {\sum}\limits_i \left\langle {{\bf{S}}_{{\bf{R}}_i + {\bf{R}}_l} \cdot {\bf{S}}_{{\bf{R}}_i}} \right\rangle .$$

(2)

Due to the possible spontaneous symmetry breaking we adjust the structure factor by the degeneracy of the momentum. To characterize the relative strength of the dominant fluctuations we project the relative intensity of the dominant static structure factor onto the range [0, 1] using the following

$${\mathrm{Intensity}} = 1 - \frac{{d_{{\bf{q}}_{{\mathrm{sub}}}}S\left( {{\bf{q}}_{{\mathrm{sub}}}} \right)}}{{d_{{\bf{q}}_{{\mathrm{max}}}}S\left( {{\bf{q}}_{{\mathrm{max}}}} \right)}}$$

(3)

where d_q is the degeneracy of momentum q, q_max is the momentum with the largest d_qS_q, and q_sub is the momentum with the second largest (subdominant) d_qS_q.

The Raman susceptibilities for B_1g, B_2g, and A_1g symmetries for non-zero temperatures were calculated using the Fleury-Loudon scattering operator²⁰ given by

$${\cal O} = \mathop {\sum}\limits_{i,j} {\kern 1pt} J_{ij}\left( {\widehat {\bf{e}}_{{\mathrm{in}}} \cdot \widehat {\bf{d}}_{ij}} \right)\left( {\widehat {\bf{e}}_{{\mathrm{out}}} \cdot \widehat {\bf{d}}_{ij}} \right){\bf{S}}_i \cdot {\bf{S}}_j$$

(4)

where J_ij are the exchange interaction values used in the Hamiltonian, $\widehat {\bf{d}}_{ij}$ is a unit vector connecting sites i and j and $\widehat {\bf{e}}_{{\mathrm{in/out}}}$ are the polarization vectors. For the symmetries calculated we use the polarization vectors

$$\begin{array}{*{20}{c}} {\widehat {\bf{e}}_{{\mathrm{in}}} = \frac{1}{{\sqrt 2 }}\left( {\widehat {\bf{x}} + \widehat {\bf{y}}} \right),\hskip 4pt \widehat {\bf{e}}_{{\mathrm{out}}} = \frac{1}{{\sqrt 2 }}\left( {\widehat {\bf{x}} + \widehat {\bf{y}}} \right)\,{\mathrm{for}}\,A_{1g} \oplus B_{2g},} \\ {\widehat {\bf{e}}_{{\mathrm{in}}} = \widehat {\bf{x}},\hskip 4pt \widehat {\bf{e}}_{{\mathrm{out}}} = \widehat {\bf{y}}\,{\mathrm{for}}\,B_{2g},} \\ {\widehat {\bf{e}}_{{\mathrm{in}}} = \frac{1}{{\sqrt 2 }}\left( {\widehat {\bf{x}} + \widehat {\bf{y}}} \right),\hskip 4pt \widehat {\bf{e}}_{{\mathrm{out}}} = \frac{1}{{\sqrt 2 }}\left( {\widehat {\bf{x}} - \widehat {\bf{y}}} \right)\,{\mathrm{for}}\,B_{1g},} \end{array}$$

(5)

(where $\widehat {\bf{x}}$ and $\widehat {\bf{y}}$ point along the Fe–Fe directions). We use this operator to calculate the Raman response R(ω) using the continued fraction expansion⁴⁹, where R(ω) is given by

$$R(\omega ) = - \frac{1}{{\pi Z}}\mathop {\sum}\limits_n {\mathrm{e}}^{ - \beta E_n}{\mathrm{Im}}\left( {\left\langle {{\mathrm{\Psi }}_n} \right|{\cal O}^\dagger \frac{1}{{\omega + E_n + i\epsilon - {\cal H}}}{\cal O}\left| {{\mathrm{\Psi }}_n} \right\rangle } \right)$$

(6)

with Z the partition function. The sum traverses over all eigenstates Ψ_n of the Hamiltonian ${\cal H}$ having eigenenergies E_n < E₀ + 2J₁ where E₀ is the ground state energy. The Raman susceptibility is given by $\chi \prime\prime (\omega ) = \frac{1}{2}\left[ {R(\omega ) - R( - \omega )} \right]$. The dynamical spin structure factor was calculated using the same method with ${\cal O}$ replaced with $S_{\bf{q}}^z = \frac{1}{{\sqrt N }}\mathop {\sum}\nolimits_l {\mathrm{e}}^{i{\bf{q}} \cdot {\bf{R}}_l}S_l^z$. Note added in proof: More details about the numerical study of the model can be found in ref. ⁵⁰.