Spin wave excitations in low dimensional systems with large magnetic anisotropy

The low-energy excitation spectrum of a two-dimensional ferromagnetic material is dominated by single-magnon excitations that show a gapless parabolic dispersion relation with the spin wave vector. This occurs as long as magnetic anisotropy and anisotropic exchange are negligible compared to isotropic exchange. However, to maintain magnetic order at finite temperatures in extended systems, it is necessary to have sizable anisotropy to open a gap in the spin wave excitation spectrum. We consider four real two-dimensional systems for which ferromagnetic order at finite temperature has been observed or predicted. Density functional theory calculations of the total energy differences for different spin configurations permit us to extract the relevant parameters and connect them with a spin Hamiltonian. The corresponding values of the Curie temperature are estimated using a simple model and found to be mostly determined by the value of the isotropic exchange. The exchange and anisotropy parameters are used in a toy model of finite-size periodic chains to study the low-energy excitation spectrum, including single-magnon and two-magnon excitations. At low energies, we find that single-magnon excitations appear in the spectrum together with two-magnon excitations. These excitations present a gap that grows particularly for large values of the magnetic anisotropy or anisotropic exchange, relative to the isotropic exchange.


I. INTRODUCTION
The appearance of ferromagnetic order at finite temperatures in two-dimensional systems requires the existence of magnetic anisotropy, so that a gap appears in the magnon (spin wave) excitation spectrum.Indeed, the magnitude of this gap is determinant for the Curie temperature T C that marks the quenching of ferromagnetism.T C also depends on the exchange coupling between spins, the magnitude of the spins, and the number of nearest neighbours 1 .In the case of localized spins at the atomic sites, e.g., in insulating systems, one typically uses a Heisenberg-like spin Hamiltonian that includes, at least, nearest neighbours isotropic exchange interactions between sites, as well as the single-ion anisotropy at each site and, additionally, anisotropic exchange or even antisymmetric exchange 2 .Single-ion anisotropy is determined by the crystal field around the magnetic atoms and the strength of spin-orbit interaction, this latter being also relevant for the anisotropic exchange.However, explaining the particular type of exchange interactions for each and every system is far from trivial and different models have been proposed 3 .
From the experimental side, there are a bunch of exper-imental techniques that can be used to explore the magnetic order of thin films with magnetic anisotropy and, thus, the signatures of both single and multi-magnon processes.This includes X-ray magnetic circular dichroism (XMCD) 4 and its depth-resolved variants 5,6 , neutron scattering, which has the additional advantage of direct access to the dispersion relation on the whole wave vector space 7,8 , magnetometry using SQUID 9 , Raman spectroscopy [10][11][12][13] , ferromagnetic resonance (FMR) [14][15][16] or local techniques, such as spin-polarized scanning tunneling microscopy 17 and inelastic tunneling spectroscopy [18][19][20] .Interestingly, only recently the effect of magnons in atomically thin layers has been accessible experimentally through Raman scattering 21 , allowing for instance to determine the optical selection rules established by the interplay between crystal symmetry, layer number, and magnetic states in CrI 3 .Nowadays, the use of first-principles calculations permits obtaining accurate total energy values for different spin configurations of a given magnetic system.A proper choice of these configurations with different spin orientations allows fitting parameters defined in the spin Hamiltonian by mapping total energy differences, although with some limitations 22 .In this work, we have consid-ered four two-dimensional systems of different kinds: 1) single septuple layer of the magnetic topological insulator MnBi 2 Te 4 , 2) single triple layer of a Mott-Hund´s magnetic insulator CrI 3 , 3) 2D metal-organic coordination network Fe-DCA on the Au(111) substrate, and 4) monolayer of Co on the heavy metal substrate Pt(111).The four systems have been shown [23][24][25][26][27] to present ferromagnetic order at finite temperatures and, therefore, are interesting to consider.Our aim here is to understand the nature and origin of the anisotropy 28,29 , as well as the spin wave excitation spectrum and Curie temperatures.In order to achieve a qualitative understanding, we use an auxiliary toy model of finite ferromagnetic spin chains with exchange and anisotropy parameters extracted from the previous mapping of the spin Hamiltonian.

II. DENSITY FUNCTIONAL THEORY CALCULATIONS AND MAPPING TO SPIN HAMILTONIAN
In this section we briefly describe the way density functional theory (DFT) calculations have been performed for the four magnetic systems of interest, as well as the extraction of the exchange and anisotropy parameters that are later on used to model the spin wave excitations.
DFT calculations were performed using the projector augmented-wave method 30 , implemented in the VASP code [31][32][33] .The generalized gradient approximation (GGA) was employed to describe the exchangecorrelation energy 34 .To take into account the strongly localized character of the 3d-states (except those of Co in the Co/Pt system, which is metallic), we resorted to the GGA + U approach 35 .The U eff = U − J 36 parameter was chosen to be equal to 2, 6, and 5.34 eV for Cr, Fe, and Mn 3d-states, respectively, following the literature 24,26,28 .For each system of interest, the atomic coordinates were relaxed using a force tolerance criterion for convergence of 0.01 eV/ Å.
The calculations were performed in a model of repeating films separated by a vacuum gap of a minimum of 10 Å.In the Co/Pt and MnBi 2 Te 4 systems [CrI 3 and Fe-DCA/Au(111) systems], in which the 3d-atom layer has a hexagonal [honeycomb] lattice, the (1 × √ 3) [(1 × 1)] magnetic periodicity was used, so that the unit cell contains two magnetic atoms in each case (see Fig. 1).
The values of the nearest-neighbor Heisenberg exchange coupling constants J, anisotropic exchange λ [(J z − J) in Table I], and single-ion anisotropy D parameters [see equation (1) below where the spin Hamiltonian is explicitly given] were obtained via accurate static totalenergy calculations, performed for the optimized crystal structures.While J is obtained via a scalar relativistic DFT calculation of ferromagnetic and antiferromagnetic configurations, determining λ and D requires (i) inclusion of spin-orbit coupling and (ii) consideration of the latter two magnetic configurations for both in-plane and outof-plane moment directions.In all of these static total- energy calculations we used a strict convergence criterion of 10 −8 eV, the Γ-centered k-points sampling of the 2D Brillouin zone, and the tetrahedron integration method with Blöchl corrections.The specific Brillouin zone samplings are given in the Supplementary Information.
For the systems with a honeycomb lattice (Cr and Fecontaining systems), there are two magnetic atoms in the surface unit cell, and each atom is coordinated with three nearest neighbors (Fig. 1a).These three nearest neighbours can have ferromagnetic −3JS 2 or antiferromagnetic 3JS 2 coupling between them, which translates into a total energy difference of 6JS 2 .For Co/Pt(111) and MnBi 2 Te 4 , the coordination of Co and Mn atoms is six and we have to use a (1 × √ 3) unit cell that contains two magnetic atoms.For the cell shown in Fig. 1b, the ferromagnetic coupling includes the six equal −JS 2 contributions, i.e., a total of −6JS 2 , while the antiferromagnetic coupling includes two spins coupled ferromagnetically and four spins coupled antiferromagnetically, i.e., a total of 2JS 2 , which translates into a total energy difference of 8JS 2 .This number of aligned or antialigned spin pairs counting is the same for the anisotropic exchange contribution, while the single-ion anisotropy term includes only two magnetic atoms per cell.
The results obtained for the four systems under study are given in Table I.In all cases, the largest energy scale corresponds to the isotropic exchange strength, which can be as large as tens of meV in the Co/Pt case but of the order of meV, or even less, in the three other cases.As shown in the next section, the isotropic exchange is crucial to determine the Curie temperature.The anisotropic exchange λ is, at least, two orders of magnitude smaller than the isotropic exchange J.The single-ion anisotropy D is larger than the anisotropic exchange, except for CrI 3 .The largest value of D/J = 0.26  I. 2D systems under study with the spin S of each magnetic atom (within brackets), the corresponding values of exchange couplings J, axial anisotropy D and anisotropic exchange λ = Jz − J as extracted from DFT according to the mapping with Hamiltonian (1).TC was calculated according to the expressions provided in Ref. 37 .
corresponds to the case of MnBi 2 Te 4 .These variations in the anisotropic exchange and single-ion anisotropy, relative to the isotropic exchange, have consequences in the spin wave excitation spectrum that are discussed in section IV C: they will induce an energy gap in the spin waves excitation spectra.

III. CRITICAL TEMPERATURES
The magnetic susceptibility of a ferromagnet obeys Curie-Weiss law where, above a critical temperature T C , it behaves essentially as a paramagnet with an enhanced susceptibility 38 , while below it, the material can keep a finite magnetization in the absence of applied field.Thus, it is important to estimate reasonable T C values for the different systems.Torelli and Olsen 37 proposed a simple parametric dependence of T C on the model parameters for different two-dimensional materials based on fitting to classical Monte Carlo simulations.The proposed T C reproduces both the critical temperature in the large magnetic anisotropy limit, and in the low anisotropy limit where the commonly used random phase approximation (RPA) 38 fails.Table I shows the critical temperatures extracted for the four materials studied in Sec.II.Notice that since D/J, λ/J ≪ 1 in all cases, the value of the critical temperature is essentially determined by the value of the isotropic exchange J, the number of nearest neighbours and the magnitude of the spin magnetic moment S 37 .However, it requires finite values of the anisotropy to differ from T C = 0. Incidentally, these calculated T C are in reasonable agreement with available data or estimated values by other methods 21,23,[25][26][27] .

IV. SPIN WAVE EXCITATIONS IN FINITE FERROMAGNETIC SPIN ARRAYS
We consider a system of N interacting spins described by the following spin Hamiltonian: Values of the spin S, relative local magnetic anisotropy D/J and anisotropic exchange (J z − J)/J for four different spin rings whose parameters resemble those of Co/Pt, CrI3, Fe-DCA/Au and MnBi2Te4, respectively.For simplicity, we have denoted J ⊥ = J.
where J z ij ≥ 0 represents the exchange coupling of the longitudinal component z between spins i and j, J ⊥ ij is the exchange coupling on the perpendicular plane, and thus, λ ij = J z ij − J ⊥ ij corresponds to the anisotropic exchange between the pair of spins (i, j).In addition, D is the uniaxial magnetic anisotropy (D > 0 conforms with an easy axis magnetic anisotropy).The operators S z j and S ± j correspond to the z-component of the j-spin operator and the corresponding ladder operators S ± j = S x j ± iS y j .Here the sum over indices ⟨i, j⟩ denotes the sum over the first neighbors.This Hamiltonian can represent either a one-dimensional or a multidimensional system.For onedimensional rings, we impose periodic boundary conditions where S 1 = S N .
For the description of the spin waves, it is convenient to introduce the Fourier transforms of the exchange interaction where we have assumed that the local S j -spin is at the R j position, and the dimensionless parameter Notice that the periodicity of the lattice implies that J z/⊥ (q) does not depend on the k-lattice site.For simplicity, we assume that all spins are of equal magnitude S ≥ 1/2 and a uniform exchange interaction between first-neighbors, i.e., J z/⊥ kj = J z/⊥ if j and k are firstneighbors and zero otherwise.
The ground state |G⟩ of Hamiltonian (1) for J z , J ⊥ > 0 corresponds to the Weiss state that can be written as |G⟩ = N i=1 | − S⟩, being S the total number of spins in the lattice.?For finite-size systems, one can in principle solve Hamiltonian (1) and find the discrete energy levels ϵ n and the corresponding eigenvectors |ϵ n ⟩.These states can be written in terms of the spin configurations |α⟩ ≡ |m 1 , m 2 , . . ., m N ⟩, where m i is the spin quantum number associated with S z i .Then, we can write the eigenvectors of H as |ϵ n ⟩ = α ⟨α|ϵ n ⟩|α⟩.

A. Single-spin excitation
Single magnons are delocalized spin waves whose total spin differs in one unit of angular momentum with respect to the ground state of the magnet.In our finite-size system described by Eq. (1), we define the normalized local excitation of the ℓ spin of the Weiss state as Here we have included the normalization factors similarly to what is done in usual spin wave theory (see supporting information for further details).A finite-size spin wavelike state containing a single magnon has the form with R ℓ the position of the ℓ-spin.Notice that, for a closed one-dimensional chain, the momentum q is quantized due to the periodic boundary conditions (see supplementary material).Thus, we can define the overlap of the eigenvectors of our spin Hamiltonian with the |q⟩ spin wave as where one finds that The excitation energy of a single spin wave can be found from semiclassical arguments and it reads as 1 : where ξ = J ⊥ (q)/J z (q).Hence, the spinwave spectrum presents an energy gap of magnitude S J ⊥ (0)/ξ − J ⊥ (0) +D (2S − 1) at zero magnetic field.In finite-size systems, the momentum q is quantized into a set of discrete values {q i } and, thus, for each energy level ϵ n there will be a given number N n of states with overlaps ψ n qi above a given threshold ε.Moreover, when finite-size effects are negligible, not only does the lowest energy ϵ n reproduce the spin-wave dispersion E SW (q) perfectly 39 but, as we shall see below, it leads to a single state with non-negligible overlap for E SW (q i ) = ϵ n .

B. Two-spin excitation
Two-magnon excitations of a ferromagnet correspond to spin waves differing by two units of angular momentum from the totally aligned states, a set of states that comprises an invariant subspace under the dynamical motion defined by the Heisenberg Hamiltonian 40 .These twomagnon excitations depend on the momentum of the pair of magnons k 1 and k 2 or, alternatively, on the total momentum Q and the difference q, where k 1 = Q/2+q and k 2 = Q/2 − q.It embraces a continuum of states corresponding to two non-interacting magnons.In addition, there may be up to two branches below the continuum that correspond to two-magnon bound states resulting from the magnon-magnon interaction that confines both magnons together 40 .
Two-magnon excitations are essential to understand the low-temperature Raman scattering (Q = 0) where they manifest as resonances close to the bound states energies 41 , the enhanced relaxation of uniform modes observed in FMR 15 , or the strongly temperature-dependent peaks in neutron scattering [42][43][44] .Moreover, exchange Hamiltonians like (1) display two-magnon bound states bellow a two-magnon continuum 45,46 .These multiplemagnon bound states can be detected in the resonant spectra 47,48 or in electron spin resonance (ESR) experiments 49 , with intensities that decrease exponentially upon lowering the temperature or increasing the multiplicity of the magnon mode.However, due to the dominant statistical weight of the two-magnon continuum in thermodynamic measurements, the common forbidden transition character of the Q = 0 point explored in most spectroscopic techniques 49 , the inherently weak intensity of the two-magnon cross-section makes the observation of these bound states quite challenging.Nevertheless, they can play an important role in the entanglement and non-equilibrium dynamics as explicitly in trapped ions quantum simulators 50,51 .
The problem of the two magnons has been solved exactly by Tonegawa 52 in one and two dimensions (the main results can be found in the supplementary material).He found that there are two branches of the energy of the two-magnon bound state E 2M (Q).He observed that below a certain threshold momentum Q th , there is a pair of complex conjugate solutions and one real solution.This indicates that the corresponding bound-state branch is no longer stable and the bound-state dispersion crosses the continuum.Moreover, for certain regions of the parameter space, negative values of the (excitation) energy of the two-magnon-bound state can be obtained.This reflects that one-dimensional ferromagnets may show unstable fully aligned states along the quantization axis z even for parameters for which the single magnon excitation energy is positive.In addition, it is possible to find single-magnon energies E SW (q) above the two-magnon bound states E 2M (Q).
In particular, he found that the amplitude ψ Q (ℓ, ℓ ′ ) for finding two spin-deviations at the ℓ-th and ℓ ′ -th sites decays quite fast with |ℓ − ℓ ′ |, with ψ Q (ℓ, ℓ ′ ) = δ ℓ,ℓ ′ at the Brillouin zone boundary for one of the two energy branches, the so-called Ising type, while ψ Q (ℓ, ℓ ′ ) = δ ℓ ′ ,ℓ+1 for the second one, called as Bethe type.For an arbitrary value of Q, one of the two magnon bound states has a larger amplitude on nearest neighbors while the 0 .0 0 . 2 0 .4 q a / π C ( S = 2 ) q a / π q a / π q a / π A ( S = 1 )  II.Each non-zero momentum q = 2πn/(N a) is associated with a single chain of size N .The symbol's sizes are proportional to the weights ψ n q for the single magnon, and only states with projections ψ n q above a threshold value of ϵ = 0.1 have been depicted.
second one has a larger amplitude on the same site, but they can not be classified rigorously into either type 52 .In the same way, the amplitudes ψ Q (ℓ, ℓ ′ ) decay much slower with |ℓ − ℓ ′ | as we move away from the first Brillouin zone boundaries.
Let us formulate the problem using the ideas of spinwave theory.The normalized two-spin excitation corresponding to sites ℓ and ℓ ′ is defined as where χ S ℓℓ ′ = 1/2S if ℓ ̸ = ℓ ′ and χ S ℓℓ = 1/ 4S(2S − 1).Similarly to the single spin excitations, we can define a two-magnon wave function as where the normalization condition imposes that We now define the overlaps of the eigenvectors |ϵ n ⟩ with the two-magnon waves where n.n., ℓ indicates that the sum is realized over the nearest neighbors of spin ℓ.Markedly, while for the Isingtype the two-magnon is fully defined by a total momentum Q, the Bethe-type two-magnon state also depends on the momentum difference q and so does also the overlap ψ Bethe Q,q (n).These overlaps contain, in addition to the two-magnon bound state contribution, a portion of two unbound magnons that, in average, occupy either the same or neighboring sites for each type, respectively 51 .As it happens for the single magnon, for each energy level ϵ n , there will be a finite number of states with overlaps ψ n Q,q above a given threshold.

C. Results for a one-dimensional toy model
To further extend our understanding of the spin waves and their connection with the magnetic anisotropy, we will analyze the energy spectrum of a closed chain of anisotropic spins.We stress that, from the practical point of view, the main effect of the change in dimensionality is to extend the sum over the first neighbors.For the ring, the momentum q of the spin waves is quantized, i.e., q n = 2πn/(aN ), n = 0, . . ., N − 1, where a is the distance between spins in the chain.Here, we have considered four scenarios that mimic the properties of the systems studied in Sec.II.The model parameters of each of these systems are summarized in Table II.The four cases A-D correspond to situations of increasing anisotropy, starting from the almost isotropic Heisenberg-like S = 1 chain in A and finishing with the largest anisotropy-induced gap for case D.
Figure 2 summarizes the main results of the spin waves in close rings with the parameters of Table II.The single spin wave solution E SW (q) has been plotted as a reference for the four cases.Case A represents the closest case to the isotropic Heisenberg spin chain.As observed, it does not show any apparent gap in the spectrum of spin waves.The presence of a local anisotropy D ̸ = 0 or an anisotropic exchange induces a finite energy gap at k = 0 of magnitude E gap SW = 2SJ ⊥ (1/ξ − 1) + D (2S − 1).This (n) for the Ising-type (red squares) two-magnon spin waves (only weights above 0.1 have been represented).The symbol's sizes are proportional to the weights.The single magnon solution ESW (Q) (dashed line) has been added as a reference.Each non-zero momentum q = 2πn/(N a) is associated with a single chain size N .

is clearly observed in cases C and D (notice that the gap has been scaled by JS).
Let us now analyze what happens in finite-size rings and, in particular, to the overlaps with the singlemagnon.In all cases, the single-magnon excitations overlap exactly with the lattice solution due to the preserved translational symmetry of the ring.Hence, the only effect of the finite size is the quantization of the momentum q.Interestingly, only the energy level with ϵ n −ϵ 0 = E ring SW (q) shows a significant overlap with the single magnon spin wave, with ϵ 0 the ground state energy.This is depicted in Fig. 2 by the size of the corresponding symbols.
As explained in Sec.IV B, two-magnon excitations can be relevant in the dynamics of these spin systems, in particular, in the low-temperature limit.With this idea in mind, we have also explored the two-magnon excitations of the above systems, see Fig. 3.We have represented the two-magnon continuum by the shaded area together with the energy of the two-magnon bound states (see Supplementary Material for detailed expressions of continuum E 2SW (Q, q) and bound states energies E 2M (Q)).For the parameters in A-D, only one of the solutions corresponding to the bound states is positive and real.Notice that the two-magnon bound state energies constitute a lower bound to the energy of the continuum, and though the bound state energies approach the bottom of the continuum as Q → 0, it may or may not merge with it 52,53 .
Interestingly, an energy gap can also be opened on the two-magnon-bound states' excitations, and it can lie either above (as in B) or below (D) the single-magnon excitation, a situation typically found for large enough anisotropy D and where a large impact of two-magnon excitations is expected 52 .Indeed, Tonegawa 52 argued that this situation indicates that the states in which all the spins aligned along the z-axis are unstable.
Let us focus now on the projection of the energy eigenstates of (1) on the two-magnon states for finite rings.To facilitate the interpretation, only the N = 10 energies ϵ n have been depicted for q = 0, the only case where the different chain sizes correspond to the same momentum.Importantly, the long wavelength limit maintains the single magnon as the lowest energy excitation.As advanced in Sec.IV B, the overlaps defined in Eq. ( 10) contain a contribution of single-magnon type clearly reflected in Fig. 3.For almost isotropic rings (case A), there is a quite large overlap with the two-magnon-bound state for arbitrary momentum, which shows both Ising and Bethelike character in the low-k region.
The situation changes quite drastically when anisotropy is included.First, the quantized two-magnon excitation deviates appreciably from the semiclassical curve of the infinite lattice.In particular, the single magnon-like contribution lies above the single-magnon semiclassical solution, well within the two-magnon continuum and having a prominent Ising-type character, see cases B-D.Second, the two-magnon bound states, which for the isotropic case almost overlap with the semiclassical solution (see case A), now also move upwards in energy (cases B-D) when anisotropy increases.Indeed, for large anisotropy the long wavelength lowest energy excitation in the semiclassical limit corresponds to two-magnon excitations and, therefore, it represents a qualitative difference with respect to the isotropic case.Finally, contrary to the single magnon case where there is a single eigenstate of H that does overlap with the spin wave and, hence, has a well-defined character, now a large number of eigenstates with both Ising-and Bethe-type characters display finite but small overlap with the two-magnon sates (9).

V. SUMMARY AND CONCLUSIONS
Here, we have estimated the Curie temperatures of different two-dimensional materials based on the simple phenomenological expressions provided by Torelli and Olsen 37 .In so doing, the isotropic and anisotropic exchange coupling constants, as well as the single-ion anisotropy, are obtained by mapping the total energy differences obtained from DFT calculations for different spin configurations with a spin Hamiltonian model.
Later on, in order to explore the low-energy magnetic excitations in systems with magnetic anisotropy, we use a simple finite-size periodic chain model of spin Hamiltonian with the previously calculated exchange coupling constants and anisotropy.The corresponding eingenstates are found by exact diagonalization.Their projection onto single-magnon and two-magnon states reveals important changes in the spin wave excitation spectrum for large values of the magnetic anisotropy, both singleion and anisotropic exchange.We not only reproduce the well-known opening of a gap in the single-magnon excitations, but also find the importance of two-magnon excitations in the low-energy spin wave excitation spectrum.In addition, two-magnon excitations of different kinds, including two-magnon bound states 54 , are found.These results suggest that some of the two-dimensional materials considered in this work, particularly those with large magnetic anisotropy, may present a lower two-magnon energy gap compared to the usual single-magnon excitation.Therefore, we speculate that, in systems with large magnetic anisotropy, multi-magnon processes can play an important role in determining the low-energy magnetic excitations, although the intensity of the corresponding signal is expected to be weaker 55 and, thus, are relevant in the interpretation of the low-temperature properties of these two-dimensional ferromagnets.Its confirmation, of course, requires the use of precise and sensitive techniques, like Raman scattering 21 or ferromagnetic resonance [14][15][16] .

SUPPORTING MATERIAL SUPERCELL STRUCTURE AND BRILLOUIN ZONES IN DFT CALCULATIONS
Co/Pt.
For Co monolayer on Pt(111), the substrate contained five atomic layers.The two lowermost layers of Pt were fixed upon structure optimization, while all other layers in the cell were allowed to relax along the out-of-plane zdirection.The optimized Pt bulk lattice parameter a 0 = 2.813 Å was used to define the in-plane (1 × √ 3) cell of rectangular shape.The Co monolayer was expanded to match the Pt lattice constant.The static total-energy calculations were performed using the 43×25×1 k-mesh.

CrI3.
The optimized bulk in-plane lattice constant a 0 = 7.1 Å was used for the free-standing CrI 3 trilayer.This corresponds the Cr-Cr separation of 4.1 Å.The Cr-I interlayer distances were relaxed.The static total energy calculations were performed using the 25 × 25 × 1 k-grid. Fe-DCA/Au(111).
It was simulated using the experimentally found (4 26 , with the Fe atoms residing in the substrate's fcc hollow sites.The optimized Au bulk lattice parameter a 0 = 2.915 Å was used to construct the in-plane (4 √ 3 × 4 √ 3)R30 • cell (note that this is the periodicity with respect to the substrate, while for the Fe-DCA itself it is (1 × 1) honeycomb).This corresponds to the Fe-Fe distance of 11.66 Å.Three Au layers were used to simulate the Au(111) substrate, so that the Fe-DCA/Au(111) cell contained 224 atoms.The atoms of the lowermost Au layer were fixed during the structural relaxations, while the other two as well as the Fe-DCA layer were allowed to relax.The complex non-co-planar geometry that Fe-DCA acquires when placed on top of Au(111) is described in detail in Ref. 26 .The static totalenergy calculations were performed using the 7 × 7 × 1 k-mesh.

MnBi2Te4.
The in-plane lattice constant of the MnBi 2 Te 4 septuple-layer-thick film 23 was fixed to that of the bulk material (a 0 = 4.336 Å23,24 ) to construct the rectangular (1 × √ 3) cell.The interlayer distances were optimized.The static total-energy calculations were done using a k-point grid of 33 × 21 × 1.
Extracting the model parameters J, D, and λ.
Based on an energetic analysis of the solutions, one can obtain the magnetic anisotropy energy, as well as the isotropic and anisotropic exchange coupling between neighbour spins localized at the 3d magnetic atoms (Co, Cr, Fe and Mn) with different spins (S = 1, 3/2, 2, and 5/2, respectively).For the four systems under study, we find that they present out-of-plane anisotropy and favorable ferromagnetic coupling between the 3d magnetic atoms' spins.In first approximation, we can obtain the nearest neighbor (isotropic and anisotropic) exchange coupling strengths and single-ion anisotropy from total energy differences of magnetic configurations with parallel and anti-parallel spins along directions perpendicular and parallel to the plane that contains the magnetic atoms 28 .

SPIN WAVE THEORY
Let us revise the spin wave theory applied to the following spin Hamiltonian describing a ferromagnet in the presence of a homogeneous applied magnetic field.We will derive the energy of different non-interacting spin waves.We first notice that the doubly-degenerate ground state of the system corresponds to a state with all spins aligned.

Single magnon spin waves
If we restrict ourselves to low-energy excitations, the deviations of the different spins from the z-axis will be small, and one can treat them using spin-wave theory 38 .A spin wave with minimal spin deviation with respect to the Weiss state can be written as a linear combination of states where the n-state is flipped, Eqs.(3) and (4) in main text.Thus, we are looking for states of the form |SW 1 ⟩ = n A n S + n |G⟩, with the normalization condition n 2S|A n | 2 = 1.We assume that the exchange integrals J z ij and J ⊥ ij depend only on the difference between the spins positions R j − R i , so that for the spins located on an infinite lattice (in any dimensions), the Hamiltonian has translational symmetry.Thus, each eigenstate of this Hamiltonian must be a basis for an irreducible representation of the translational group of the lattice 1 , so that it satisfies Bloch theorem.Hence, the |q⟩ state (4) is an eigenstate of H.It is worth mentioning that for finite periodic systems, i.e., a ring in one dimension or a torus in two dimensions, the situation is fully analogous to the infinite lattice limit with the exception of the quantization of the momentum.?Taking into account the following relations for i ̸ = j one can evaluate the action of all the terms in our Hamiltonian on the |q⟩ state.In fact, one can demonstrate that 1,56 where the energy of the Weiss state can be written as and the excitation energy E SW (q) with respect to the ground state is given by Eq. ( 7) in the main text.
In a honeycomb lattice there are only three first neighbors for each lattice site, and the single-magnon spin wave has two possible solutions to E 0 + E SW (q) 28 and where ε 0 = 3SJ ⊥ /ξ + D (2S − 1) + gµ B B and f (q) = 1 + e iq•a1 + e iq•a2 is the form factor of the honeycomb lattice, with a 1 and a 2 the lattice vectors of the triangular lattice.
In the case of the 2d-triangular lattice, there are six first neighbors lying on the corners of a regular hexagon around each of the lattice sites, and only one element in the crystallographic basis.Introducing the lattice vectors a 1 = a î and a 2 = a/2 î + a √ 3/2 ĵ, the single magnon excitation energy is then given by E trian SW (q) = ε t 0 + J ⊥ S 6/ξ − f t (q) , (S6) with ε t 0 = D (2S − 1) + gµ B B and f t (q) = 2 [cos(q • a 1 ) + cos(q • a 1 ) + cos (q • (a 1 − a 2 ))] . (S7) Periodic chain: ring.
A periodic one-dimensional chain is in essence a closed ring.The only difference with respect to the infinite chain is the momentum quantization as a result of the boundary conditions.This leads to q n = 2πn/(aN ), n = 0, . . ., N − 1, where a is the distance between spins in the chain.The single-magnon energy ϵ q in Eq. ( 7) for a ring with first-neighbor exchange J z,⊥ then reads as Here we just recapitulate the main results of Tonegawa 52 (adapting the notation).Let us introduce the following dimensionless parameters: δ = 2D(2S − 1)ξ 4J ⊥ S (S9) First, the two non-interacting magnons give place to an energy continuum where the sum is over the first neighbors of the lattice and Q is the total momentum of the pair of magnons and whose wavevectors are k 1 = Q/2 + q and k 2 = Q/2 − q.
Thus, for a one-dimensional ring one finds an expression that equals twice the energy of a single magnon, i.e., In addition, there are two solutions below this energy continuum that correspond to bound states with localized wavefunctions of the two spin excitations whose wavefunctions can be written as with N the total number of spins.Thus, the amplitude ψ(l, l ′ ) for finding the two-spin deviations at the l and l ′ sites will be given by 4Sa(l, l ′ )/ √ N when l ̸ = l ′ and 2 S(2S − 1)a(l, l ′ )/ √ N when l = l ′52 .When cos (Q • r j ) = 0, one of these solutions corresponds to the single-ion type two-magnon bound state, with ψ(l, l) ̸ = 0 and ψ(l, l ′ ) = 0 for l ̸ = l ′ .This solution is known as "Ising type" bound state (this bound state only appears for S > 1/2).The other solution corresponds to a two magnon wavefunction where ψ(l, l ′ ) ̸ = 0 for l ′ a first neighbor of l, and ψ(l, l ′ ) = 0 it is known as "Bethe type" bound state.For a general Q, this simple separation is no longer possible.

FIG. 1 .
FIG. 1. Top view of the two-dimensional structures of the systems considered (only magnetic atoms are drawn) showing the unit cells used in the calculations for the (1×1) honeycomb (a) and (1 × √ 3) hexagonal (b) layer cases, both containing two atoms in the unit cell.The blue and red colors represent opposite spins in an antiferromagnetic configuration to illustrate the number of first neighbours with aligned and anti-aligned spins for each lattice: three anti-aligned spins in the honeycomb (a), while four anti-aligned and two aligned spins in the hexagonal (b).Notice that in a ferromagnetic configuration of spins the number of first neighbours with aligned spins coincides with the coordination, i.e, three in the honeycomb and six in the hexagonal lattice.

FIG. 2 .
FIG.2.Single magnon excitations of finite-size rings.Single-magnon energies (dashed lines) as calculated from ESW (q) and from the diagonalization of the finite-size Hamiltonian (1) for a one-dimensional ring (black circles) corresponding to A, B, C, and D systems in TableII.Each non-zero momentum q = 2πn/(N a) is associated with a single chain of size N .The symbol's sizes are proportional to the weights ψ n q for the single magnon, and only states with projections ψ n q above a threshold value of ϵ = 0.1 have been depicted.