Thermal evolution of spin excitations in honeycomb Ising antiferromagnetic FePSe3

We use elastic and inelastic neutron scattering (INS) to study the antiferromagnetic (AF) phase transitions and spin excitations in the two-dimensional (2D) zig-zag antiferromagnet FePSe$_3$. By determining the magnetic order parameter across the AF phase transition, we conclude that the AF phase transition in FePSe$_3$ is first-order in nature. In addition, our INS measurements reveal that the spin waves in the AF ordered state have a large easy-axis magnetic anisotropy gap, consistent with an Ising Hamiltonian, and possible biquadratic magnetic exchange interactions. On warming across $T_N$, we find that dispersive spin excitations associated with three-fold rotational symmetric AF fluctuations change into FM spin fluctuations above $T_N$. These results suggest that the first-order AF phase transition in FePSe$_3$ may arise from the competition between $C_3$ symmetric AF and $C_1$ symmetric FM spin fluctuations around $T_N$, in place of a conventional second-order AF phase transition.

Two-dimensional (2D) spin models are one of the most well-studied magnetic systems in condensed matter physics.The pioneering solution by Onsager in 1944 [1] revealed that the 2D ferromagnetic (FM) Ising model, characterized by nearest neighbor coupling J, undergoes a magnetic phase transition at k B T C = 2.269J, exhibiting an order parameter critical exponent β = 1/8.The 2D XY model, on the other hand, lacks longrange order, but it exhibits the ability to undergo a Berezinskii-Kosterlitz-Thouless transition, during which the correlation function transitions from an exponential decay to a power-law behavior as a function of distance [2,3].In real magnetic materials, spins exhibit free 3D rotation rather than adhering to binary orientations or 2D rotations, and the Hamiltonian for actual magnetic systems will take the form of a Heisenberg model plus magnetic anisotropy [4].According to the Mermin-Wagner theorem, in the case of an isotropic Heisenberg Hamiltonian with short-range magnetic exchange couplings, the thermal fluctuations are strong enough to prevent long-range magnetic order in the 2D limit at any finite temperature [5], with the correlation length diverging only at absolute zero.Conversely, in the presence of easy-axis/easy-plane anisotropy, the renormalization group flow will favor the selection of the anisotropic component of the Hamiltonian as the relevant parameter, resulting in critical behavior that aligns with predictions from the Ising/XY models [6].
The discovery of long-range magnetic order at nonzero temperatures in the 2D monolayer of several honeycomb lattice van der Waals (vdW) ferromagnets and antiferromagnets suggests a suppression of the thermal fluctuations, most likely due to the formation of an Ising-type magnetic anisotropy gap in these materials [7][8][9][10][11].To understand the microscopic origin of the long-range magnetic order in the 2D limit for different classes of vdW materials, it is therefore important to determine the magnetic properties of their 3D bulk compounds.While the FM Ising Hamiltonian on a 2D square lattice has been solved [1], the situation for Ising spin systems in 2D honeycomb lattices is more complicated.As a non-Bravais lattice, the honeycomb structure can host a few different collinear magnetic structures, including simple ferromagnets, Néel antiferromagnet with c-axis aligned moments, stripy antiferromagnet, and zig-zag antiferromagnet [Fig.1(e)], depending on the relative strengths of nearest J 1 (NN), next nearest J 2 (NNN), and next-next nearest neighbor J 3 (NNNN) magnetic exchange interactions [Fig.1(a)] [12].The first two magnetic structures have the same in-plane ordering wavevectors as the lattice vectors, and respect the three-fold rotational (C 3 ) symmetry of the honeycomb lattice.The stripy and zigzag antiferromagnetic (AF) structures, however, break the C 3 rotational symmetry of the honeycomb lattice to C 1 and fold the first Brillouin zone (FBZ) into a smaller rectangular magnetic FBZ.In these structures, there will be three magnetic domains separated by 60 • degrees, meaning that the overall magnetism still obeys the underlying C 3 lattice symmetry.The magnetic order parameter is not a direct measure of the ordered or staggered moment because it has multiple components.For example, the order parameter for the zig-zag AF order is composed of three components (Ψ 1 , Ψ 2 , Ψ 3 ), each reflecting the staggered moment along one of the three C 3 axes [13].By probing the static magnetic order and spin excitations across the magnetic phase transition, one can obtain information concerning the behaviors of the phase transition and compare the results with the expectations of the conventional Ising or XY models.
For the FM honeycomb system CrSiTe 3 , the critical behavior of the magnetic phase transition exhibits 2D Ising characteristics [14,15], with a dimensional crossover from 2D to 3D near T C [16].For other honeycomb lattice ferromagnets such as CrI 3 , VI 3 , and CrGeTe 3 , the critical behavior is more reminiscent of the tricritical mean field model [17][18][19][20][21].For AF ordered MnPS 3 , NiPS 3 and CoPS 3 , the critical behavior near T N shows 3D XY or Ising characteristics [22][23][24].However, none of these AF materials follow the genuine 2D Ising model primarily due to either a small magnetic anisotropy or an XY-type anisotropy, along with interlayer interactions that impact the magnetic phase transition [25,26].
FePS 3 and FePSe 3 constitute a potential material family for exploring an Ising-type antiferromagnetic phase transition accompanied by C 3 symmetry breaking.This study specifically centers on FePSe 3 [Figs.1(a) and 1(c)].FePSe 3 belongs to the rhombohedral R 3 space group, with hexagonal layers stacked through weak vdW interactions along the c-axis [Fig.1(b)] [27].Compared with FePS 3 in the monoclinic C/2m symmetry group, FePSe 3 preserves the C 3 symmetry and therefore is ideal to study any possible magnetic order induced in-plane C 3 symmetry breaking without significant interlayer coupling effects [28].Since FePSe 3 has an in-plane zig-zag AF structure below T N =110K, we expect an FM exchange interaction J 1 between NN, and AF exchange interactions in J 2 and J 3 between NNN and NNNN, respectively [Fig.1(a)] [27].In the 3D bulk limit, the AF ordering wavevector is Q AF = (1/2, 0, 1/2) [Fig.1(b) and 1(f)] due to a weak AF interlayer coupling J c [Fig.1(b)] [27].The magnetic Fe 2+ ion has a 3d 6 electronic orbital with four unpaired spins, giving S = 2.A strong spin-orbit coupling (SOC) induces a large magnetic anisotropy with the c-axis as the easy axis [27,29].Raman scattering experiments on bulk and monolayer FePSe 3 have observed a ∼15meV spin gap, suggesting that the system is Ising-like [11].
In this work, we use elastic and inelastic neutron scattering (INS) to study the magnetic order and spin dynamics in bulk single crystal FePSe 3 .The temperature dependence of the magnetic order parameter from elastic neutron scattering experiments suggests that the AF phase transition is first-order in nature [Fig.1(d)].In the AF ordered state, our INS experiments reveal that spin waves are gapped below ∼15meV, consistent with Raman scatteirng results [11], and are highly 2D with weak dispersion along the c-axis [Fig.1(g)].By fitting the spin-wave dispersion spectra using linear spin wave theory (LSWT) [30], we determine magnetic exchange couplings, magnetic anisotropy, and find evidence for a biquadratic term in the spin Hamiltonian.On warming above T N , we find dispersive spin excitations that can be explained by AF and FM excitations from hon-eycomb lattice clusters with the C 3 symmetry.These results are different from the expectation of a conventional AF second-order phase transition, suggesting that the long-range AF order in FePSe 3 is replaced by C 3 magnetic honeycomb lattice clusters as the temperature is raised above T N .The uncorrelated paramagnetic scattering from individual Fe ions is only established at temperatures well above T N .
Single crystals of FePSe 3 are synthesized using the chemical vapor transport method described in ref. [31].Elastic and INS experiments were respectively performed at the CORELLI spectrometer [32] on one single crystal sample and the ARCS spectrometer [33] on ∼0.5g of coaligned crystals at the Spallation Neutron Source, Oak Ridge National Laboratory.The momentum transfer Q is referenced in reciprocal lattice units (rlu) with respect to the rhombohedral unit cell of FePSe 3 with a = b =6.26 Å, c =19.71 Å.The magnetic ordering wavevector observed at Q AF = (1/2, 0, 1/2) confirms the zig-zag magnetic structure, and the temperature dependence of the order parameter is shown in Fig. 1(d).A power law fit of the order parameter with I = I 0 (1 − (T /T N ) 2β ) yields a critical exponent β = 0.063, which is much smaller than the 2D Ising model prediction β = 1/8, suggesting a first-order nature of the phase transition.This is further supported by measurements of correlation lengths which jump abruptly at T N [29].
Figures 1(g) and 2(b) show the measured spin wave dispersions of FePSe 3 at 5 K along the c axis and within the 2D honeycomb lattice plane, respectively.The overall dispersion has a band top of ∼40 meV, with an anisotropy gap of ∼15 meV at the Γ and M point.The large anisotropy gap value is also observed in the sister compound FePS 3 , which is due to the combined effects of SOC in the 3d 6 orbital of Fe 2+ and the distortion of the FeSe 6 octahedron [Fig.1(c)].For comparison, the isostructural compounds MnPSe 3 and MnPS 3 do not exhibit such large gaps because their 3d 5 electrons have quenched orbital moments [25,34].The spin waves propagating along the c-axis exhibit significantly less dispersion in contrast to the in-plane spin waves, featuring a bandwidth of approximately 1 meV for the former as opposed to a 20 meV bandwidth for the latter, indicating the presence of a very weak interlayer coupling J c [Fig.1(g)].This is consistent with the fact that the T N of FePSe 3 changes little as a function of layer numbers, suggesting that the interlayer exchange interactions have minimal effect on the magnetic phase transition in FePSe 3 [9,11].
To describe the in-plane spin wave dispersion, we first use LSWT to calculate the spin waves with the spin Hamiltonian where the bilinear (Heisenberg) exchange interaction J ij is summed over the 1st, 2nd and 3rd NN, and D z is the single-ion anisotropy with z-axis as its easy axis.
It is worth noting that within the actual system, magnetic anisotropy may stem from either single-ion effects or Ising-type exchanges.Therefore, the Hamiltonian can encompass Ising exchange terms J z S iz S jz in addition to the single-ion anisotropy.Nevertheless, within the framework of linear spin wave theory, the Ising term's impact on dispersion will be identical to that of singleion anisotropy, provided that the respective parameters are appropriately adjusted (See the supplementary information).For the sake of simplicity, we have chosen to solely include the single-ion anisotropy term here.Using a least-square-error fitting method with S = 2, we extract the magnetic exchange coupling parameters as shown in TABLE I.However, the best fitting using this model with uniform J 1 does not precisely reproduce the dispersion, especially the low-energy part perpendicular to the zig-zag direction.In fact, the exact same scenario is observed in the sister compound FePS 3 [36] where a simple Heisenberg model plus anisotropy cannot account for the spin excitations accurately.In the case of FePS 3 , two approaches have been utilized to resolve this problem, and here we apply them to FePSe 3 as well.The first is to introduce bond-dependent J 1 with J 1a bonding the parallel spins and J 1b for anti-parallel spins; The second is to introduce a biquadratic interaction where the biquadratic exchange K ij is summed over the 1st NN only.Both methods stabilize the zig-zag AF order, and with proper fitting parameters in TABLE I, they yield nearly identical dispersions that can accurately reproduce the experimental data.For the J 1a -J 1b model, the difference between J 1a and J 1b is large, indicating that a lattice distortion should accompany the magnetic phase transition as in the case of FePS 3 .However, no obvious lattice parameter change is observed across T N , indicating that the J 1a -J 1b model is unlikely to be correct.On the other hand, The existence of biquadratic interaction has been theoretically proposed in many such 2D magnetic systems with edge-shared octahedron structures [39,40], and therefore may be the more suitable model for describing the spin waves in FePSe 3 .
The situation is similar in the case of iron-based superconductors, where a biquadratic interaction has been used to account for the observed in-plane spin wave dispersions [37,38].A noteworthy fact is that the lowenergy magnons at the Γ point are linearly coupled to phonons, which introduces a magnon-polaronic gap that lifts the 2-fold degeneracy of the AF magnons [11].This suggests the necessity of introducing a magnon-phonon  coupling term into the spin Hamiltonian.However, the impact of this term is minimal, exhibiting a gap of approximately 0.6 meV, a value below the resolution of the instruments.Hence, it is omitted from equation (2).Furthermore, the required strength of the Kitaev interaction (0.03 meV) to induce the magnon-phonon gap is significantly less than that of the Heisenberg exchanges, and thus, these are not taken into account in the fitting parameters.
The first-order nature of the AF phase transition also displays itself in the spin fluctuations in the neighborhood of T N .Fig. 3 shows a summary of the INS spec- trum near T N = 110 K.At T = 101 K, in-plane spin excitations are well-defined and similar to spin waves at 5 K [Fig.3(a)].On warming to T = 107 K, spin excitations remain well defined, but become broader with an anisotropy spin gap above 10 meV [Fig.3(b)].Upon further warming to 110 K, the spin gap drastically decreases to zero at the M point, while keeping a non-zero value at the Γ point [Fig.3(c)].A cut along the [H, 0] direction at E = 5 ± 1 meV shows clearly that the M point spin fluctuations are enhanced from 110 K to 125 K, different from the expectation of critical magnetic scattering associated with a second-order phase transition.Although the spin gap closing at the M point may arise from the vanishing static magnetic moment across T N , the differences in temperature dependence of the spin gap at Γ and M points cannot be explained by a second-order AF phase transition since the Γ and M points are equivalent within the spin wave theory which requires long-range AF order and folding of the Brillouin zone below T N .This is consistent with the observed first-order transition in the AF order parameter(fig.1d).In the paramagnetic state, the rectangular magnetic Brillouin zone unfolds to the hexagonal lattice Brillouin zone.However, we find similar spin excitations as in the AF ordered state [Figs.Assuming that the spins form clusters of hexagons with zig-zag order [Fig.4(a), we calculate the neutron intensity I(Q) with where f 2 (Q) is the magnetic form factor, m, n ∈ 1, ..., 6 is the index for spins in Figs.The spin fluctuations observed at 110K encode substantial information regarding the spin correlations in FePSe 3 .The pillar-like feature of the spin excitations at the M point in Figure 3(c) indicates that the Qdependence of low-energy spin excitations remains unchanged below ∼8meV.This implies that an interval over the low-energy excitations can reasonably approximate the instantaneous spin correlations.Consequently, despite our zig-zag cluster calculation being used to reproduce the spin fluctuation in the energy range E= [3,7] meV, it is reasonable to infer that this zig-zag arrangement accurately reflects the nature of these instantaneous spin correlations.
A noteworthy observation is that even in the AF state at 4K, the dynamical spin-spin correlation centers around the FM wavevector located at Γ, rather than the AF wavevector at M [Figure 4(g) and 4(h)].The reason is that as the magnetic anisotropy grows, the AF spin wave excitations are gradually replaced by local spin-flip excitations.In an isotropic AF Heisenberg model, the magnetic excitations are spin waves that give zero structure factor at the FM wavevector, while in a pure AF Ising model, the spin-flip excitations will introduce a non-zero  I. intensity centered at the Γ point.In the case of FePSe 3 , the robust magnetic anisotropy shifts the primary spin excitations from AF spin waves to spin-flip excitations.As the temperature rises to T N , the local spin-flip excitation remains gapped, indicating that the dynamic susceptibility associated with spin-flip scattering remains nonzero.Consequently, it cannot be the primary driving force behind the first-order magnetic phase transition.Instead, the excitations of zig-zag clusters at T N become gapless, implying that the dominant component of spin configurations near T N is the zig-zag short-range order.It can be speculated that the magnetic phase transition is driven by a discontinuous condensation of these zig-zag hexagonal clusters, corresponding to the C 3 -C 1 symmetry breaking.This distinguishes it from the fluctuations observed in conventional Ising magnets that break the Z 2 time-reversal symmetry through spin-flip scatterings.
In summary, we utilized neutron scattering to study the spin excitations and spin fluctuations in FePSe 3 across T N .By analyzing the spin wave dispersion, we find evidence for biquadratic magnetic exchange interactions in the effective spin Hamiltonian terms.The microscopic origin of the biquadratic term remains to be determined.In addition, our magnetic order parameter measurements indicate that the AF phase transition is first-order in nature.From measurements of the temperature-dependent spin excitations across T N , we infer that the AF phase transition may be driven by the C 3 to C 1 symmetry breaking from low-energy spin clusters associated with zig-zag AF honeycomb lattice, respectively.Our results provide an enriched perspective on the intricate interplay between magnetic interactions and structural symmetries in the Ising magnetic systems with additional C 3 symmetry breaking.

FIG. 1 .
FIG. 1. Real/reciprocal space of FePSe3 and spin waves along in-plane and c-axis directions.(a) Schematics of the intralayer AF magnetic structure and exchange interactions.(b) Interlayer magnetic structure and coupling Jc.(c) Local Fe 2+ environment in FePSe3, with the 3d 6 spin configuration in an octahedral environment.(d) Neutron intensity of magnetic Bragg peak (1/2 0 1/2) as a function of temperature.(e) The in-plane zig-zag magnetic structure.The bold black rhombus shows the lattice unit cell and the dashed red rectangle shows the magnetic unit cell.(f) The reciprocal space with black hexagons showing the lattice FBZ and the dashed red rectangle shows the magnetic FBZ.The blue dots are the magnetic Bragg peaks from the magnetic structure shown in (e).(g) Out-of-plane spin wave dispersion along [1/2 0 L] from the experiment.
3(d), 3(e), and 3(f)].Although the AF phase transition in FePSe 3 is first-order, the short-range zig-zag order may still persist near the phase transition.A natural way to reconcile the existence of the zig-zag order and the C 3 symmetry is to combine all the zig-zag domains with equal weight [Fig.4(a)], while the size of the fluctuation domains determines the correlation lengths of the short-range order.To understand the nature of the dispersive paramagnetic spin excitations in FePSe 3 , we plot constant-E slices in the [H K] plane at E = 5 ± 2 meV [Fig.4(c)] and 11 ± 2 meV [Fig.4(d)] at 110 K, which shows the spin excitations at the M points and their connections, respectively.The low-energy spin excitations show diffusive hexagonal patterns which are hollow at the Γ point, and the high-energy excitations show plate-like hexagonal patterns.Inspired by the spin cluster methods used in refs.[45][46][47],we apply a similar analysis on FePSe 3 .
4(a) and 4(b).The calculated I(Q) is then averaged over six zig-zag spin configurations.Figures 4(d) and 4(f) show the calculated intensity from the zig-zag clusters [Fig.4(a)], which matches well with the experimental result at low energy.For the higher energy part at 11meV, due to the enhanced Γ point spin fluctuations, an additional FM cluster is required [Fig.4(b)] in order to reconstruct the plate-like excitation pattern.

FIG. 4 .
FIG. 4. Simulation on spin fluctuations at Neel temperature in FePSe3.(a,b) Zig-zag AF and FM clusters used to simulate the spin fluctuation, respectively.(c) Experimental data at T =110K and E=[3,7]meV.(d) Calculated neutron intensity of the AF cluster in (a).(e) Experimental data at T =110K and E=[9,13]meV.(f) Calculated neutron intensity of the AF cluster (a) plus the FM cluster (b).The AF clusters and FM clusters are weighted with a ratio of 2:1 in (f).(g) experimentalw3 data at T =4K and E=[14,18]meV.(h) Calculated neutron intensity at the same energy range as (g) by LSWT with parameters from TABLE I.

TABLE I .
Magnetic exchange interaction strength in different models in FePSe3.All units in meV.In the J1a-J 1b model, the J1a indicates FM interactions in the zig-zag chain, and J 1b indicates interactions between chains.J1, J2, J3 indicate the first, second, and third nearest neighbor Heisenberg exchange, respectively.K1 refers to the nearest-neighbor biquadratic interactions, and Dz stands for single-ion anisotropy.