Temperature dependence of spin-model parameters in antiferromagnets

The temperature dependence of mesoscopic spin-model parameters is derived in two-sublattice antiferromagnetically aligned systems based on Green's function theory. It is found that transversal spin correlations decrease the anisotropy terms while increasing the Heisenberg and Dzyaloshinsky--Moriya exchange interactions and the latter's contribution to the anisotropy. The obtained temperature dependences show quantitative agreement with the results for ferromagnets, and they also agree well with numerical atomistic simulations which treat the spin correlations without approximations. Possible applications of the results in multiscale modelling are discussed.


I. INTRODUCTION
The most widely studied class of antiferromagnets contains two sublattices on which the magnetic moments point oppositely to each other. Materials where the magnitude of the moments on the sublattices is different are known as ferrimagnets. Both antiferromagnets [1][2][3][4] and ferrimagnets [5,6] have attracted much attention recently as a material platform for spintronics. Their dynamics is typically orders of magnitude faster than that of ferromagnets, including higher spin-wave frequencies [1]; relativistic domain wall motion [7,8]; enhanced magnetization switching rates induced by current [9], thermal excitations [10,11] or ultrashort laser pulses [12][13][14]; and increased demagnetization speeds [15]. Antiferromagnets may transform into a phase where the sublattice magnetizations have a weak parallel or ferromagnetic component. This transformation can be achieved by applying an external field [16], increasing the temperature [17], or even in the ground state in the presence of the Dzyaloshinsky-Moriya interaction [18,19]. Due to the different temperature dependence of the sublattice magnetizations and angular momenta, they can become compensated in certain ferrimagnets, influencing the velocity and the movement direction of domain walls and skyrmions driven by spin-polarized currents or thermal gradients [20][21][22]. At the angular momentum compensation point, ferrimagnets combine the advantages of both ferromagnets and antiferromagnets: easy control and detection of their net magnetization by an external field, antiferromagnetic-like ultrafast dynamics and the potential for high-density devices.
These phenomena can often be successfully modelled by theoretical approaches agnostic to the underlying atomic structure, such as finite-temperature macrospin models like the Landau-Lifshitz-Bloch equation [23] or continuum theories. In these models, the on-site anisotropy contributions and the pair-wise interactions, such as exchange or Dzyaloshinsky-Moriya terms, are intrinsically temperature dependent due to averaging over the fluctuations and correlations of atomic spins in a finite volume. Computer simulations using atomistic spin models [24] can natu- * levente.rozsa@uni-konstanz.de rally describe the equilibrium thermodynamics and nonequilibrium dynamics of antiferromagnets and ferrimagnets, but they are considerably more resource intensive. The price paid for the reduced computational cost of the mesoscopic methods is that the temperature-dependent effective parameters of these models are difficult to determine. While first-principles methods have proven successful in calculating atomistic [25][26][27] or zero-temperature continuum [28,29] spin-model parameters, their application to finite-temperature mesoscopic models remains limited. The temperature dependence of the magnetocrystalline anisotropy energy has been calculated based on a disordered local moment scheme [30][31][32]. However, this method treats spin fluctuations on a mean-field level, since modelling correlations in density functional theory is inherently challenging. Methods for calculating exchange interactions at finite temperature have been proposed in Refs. [33,34], but the degree of spin disorder was determined from atomistic simulations in these cases. Therefore, obtaining the effective parameter values for the models require fitting to experimental data obtained in a wide temperature range which are not always available, or to the results of atomistic spin-model simulations what counteracts the reduced computational cost of the mesoscopic models.
This difficulty can be circumvented by applying analytical methods which treat correlations accurately, and can approximate these temperature-dependent parameters at low computational costs. Such methods are often based on Green's functions, where the difficulty arises in choosing the decoupling scheme, i.e., the procedure for truncating the infinite series of correlation functions of increasing order. Non-linear spin-wave theory is based on conventional diagrammatic perturbation methods developed for bosonic and fermionic systems, and is widely used for describing quantum fluctuations [35,36] at low temperatures. These methods are often inaccurate for spin models at elevated temperatures; for example, they predict a finite jump in the magnetization at the critical temperature in simple ferromagnets [37,38]. For accurately modelling the temperature-induced phase transitions, semi-empirical decoupling schemes for spin Green's functions have been developed by Tyablikov [39], known as the random-phase approximation, and by Callen [40], which have proven to be especially accurate for low and high spin values, respectively. The method has been generalized to two-sublattice antiferromagnets by Anderson and Callen [41]. The applications of these methods to two-dimensional ferromagnetic and antiferromagnetic quantum spin systems is summarized in Ref. [42], and antiferromagnets with Dzyaloshinsky-Moriya interactions have been treated within the random-phase approximation in Refs. [43,44]. These works primarily focused on the calculation of the magnetizations, the correlation functions and the magnon frequencies at finite temperatures, which can be used to describe phase transitions, but not on the connection between the atomistic and mesoscopic models. The latter topic has been investigated for ferromagnets in the classical limit of infinite spin in Refs. [45][46][47], but the corresponding investigations for antiferromagnetically aligned systems seem to be lacking.
Here, we derive the temperature dependence of the effective interaction parameters in mesoscopic models of two-sublattice antiferromagnets and ferrimagnets. We extend the Green's function theory in Ref. [41] by including all terms preserving rotational symmetry around the axis of the magnetizations, namely Heisenberg and Dzyaloshinsky-Moriya exchange interactions as well as single-ion and two-ion anisotropy terms, and discuss both the classical and quantum cases. A comparison with atomistic Monte Carlo simulations demonstrates the accuracy of the method in treating spin correlations at a fraction of the computational cost of the numerical simulations.
The paper is organized as follows. In Sec. II A, we present the self-consistency equations of Green's function theory, which we apply to derive the correspondence between the atomistic and mesoscopic models in Sec. II B. We discuss the scaling exponents of the effective parameters in Sec. II C. We apply the method to a square lattice in Sec. III and compare the predictions with atomistic simulations.

A. Green's function theory
We consider a two-sublattice magnet described by the atomistic spin Hamiltonian Here r, s ∈ {A, B} denote the two sublattices, J rs ij is the Heisenberg exchange interaction between atoms at sites i and j, D rs ij is the Dzyaloshinsky-Moriya vector, ∆J rs ij is the two-ion anisotropy, K r is the single-ion magnetocrystalline anisotropy, µ r is the magnetic moment, and B z is the external magnetic field. We note that this model describes an antiferromagnet when µ A = µ B and a ferrimagnet when µ A ̸ = µ B . S ir stands for the spin vectors; for most considerations they will be treated as classical unit vectors |S ir | = 1, since atomistic spin-model simulations used for comparison are easier to perform in the classical limit. However, at certain points the quantum-mechanical case with spin operators will be discussed. The number of unit cells in the lattice will be denoted by N c , corresponding to the number of atoms in both the A and the B sublattices. Note that fully analytical results in the limit N c → ∞ are only available when the types of interactions are more restricted; therefore, in most cases we will rely on semi-analytical techniques where lattice sums over a finite number of lattice sites must be performed. It will be assumed that in the classical ground state, spins on sublattice A are oriented along the +z direction, while spins on sublattice B point along the −z direction. To stabilize the antiferromagnetic alignment, it will be assumed that the antiferromagnetic intersublattice coupling J AB ij < 0 is dominant compared to the intrasublattice coupling and the Dzyaloshinsky-Moriya interaction, while the anisotropy prefers spin alignment along the z axis.
The equation of motion generated by the Hamiltonian H reads (2) in the classical limit, where γ is the gyromagnetic ratio. Equation (2) describes the precession of the spins around the effective magnetic field B eff ir = −µ −1 r ∂H/∂S ir . We introduce a local coordinate system where all spins are oriented along the +z direction in the classical ground state, iB denotes the ladder operators. In linear spin-wave theory, the dynamical equation is linearized around the classical ground state in the quantities S ± ir . We introduce the shorter notationsS (1) ir ∈ S + iA ,S − iB and S (2) ir ∈ S − iA ,S + iB since this linearized equation couples the + and − indices between the A and B sublattices. After performing spatial and temporal Fourier transformation ∂ t → −iω, we will use the notations Here, R i − R j denotes the vector connecting the lattice sites i and j, where R i = (x i , y i , z i ) stands for the position of the spin i in the lattice. Note that only the z component of the Dzyaloshinsky-Moriya vectors appears in these expressions. In these variables, the linearized equation of motion reads where single and double underlines denote vectors and matrices in sublattice indices r, s. Here µ = diag µ and the spin-wave Hamiltonian is with σ = [1, −1] T and σ z = diag (σ) a Pauli matrix, which appear due to the antiparallel alignment of the sublattices.
Note that the (1) and (2) components of the spins are decoupled in the linearized equation of motion, because the Hamiltonian in Eq. (1) contains all possible single-spin and two-spin terms that are rotationally invariant around the z direction. This simplification may be justified in systems with at least a threefold rotational symmetry around the Néel vector, such as the trigonal antiferromagnet α-Fe 2 O 3 [17], hexagonal ferrimagnets like GdCo 5 [48] and cubic ferrimagnets including Mn 2 Ru x Ga [49]. Even if the rotational symmetry is lower, such as in atomically thin antiferromagnetic Mn layers on Nb(110) [50], this model should provide a useful approximation if the anisotropy terms are considerably weaker than the exchange interactions.
The eigenvalues of Eq. (6) correspond to the magnon frequencies. The thermal occupation of the magnon modes enables the calculation of the sublattice magnetizations and the two-spin correlation functions at low temperatures [51]. At elevated temperatures, a higher number of magnons becomes excited, leading to a temperature-dependent renormalization of the frequencies. A self-consistent procedure for treating this renormalization based on Green's functions was introduced in Refs. [40,41], which we apply to the present system here. Details of the derivation are given in Appendix A. In this method, the single-particle excitation spectrum is given by the eigenvalues of the matrix γµ −1Γ q , where the matrixΓ q replacing the spin-wave Hamiltoniañ where S z = diag S z and • denotes element-wise multiplication. The coefficient α 0 is a phenomenological constant required for the decoupling of the Green's functions (see Appendix A); here it will be set to α 0 = 1/2 in the classical limit of the decoupling scheme proposed by Callen and Anderson [40,41]. The Φ q matrix is related to the transversal spin correlation functions via The magnon frequencies are expressed as Note that ω + q ≥ 0 and ω − q ≤ 0 if the antiferromagnetic alignment of the sublattices is stable. For example, in the antiferromagnetic limit with identical sublattices and no external field, one obtains γµ −1 The opposite signs of the frequencies describe that the modes have opposite circular polarizations, which has also been demonstrated experimentally in a ferromagnet recently [52].
Equation (8) is made self-consistent by determining the correlation functions from the spectral theorem [42] This integral is evaluated in Appendix A. The function n (ω) = k B T /ω corresponds to the occupation number of the magnon mode with frequency ω in the classical limit in units of action. Substituting this function simplifies Eq. (12) to The sublattice magnetizations are given by the Langevin function where Φ r = q Φ rr q may be interpreted as the total spin carried by the magnons on sublattice r. At zero temperature, Φ r = 0 holds as can be seen from Eq. (13), and the sublattice magnetizations are saturated S z r = 1. It is worth noting that S z r = 1 does not hold in the quantum case even for T = 0, as is already known from linear spin-wave theory [51]. Using spin operators in the quantum case, Eqs. (8) and (12) remain valid, but the function n quantum (ω) = ℏ e ℏω/(k B T ) − 1 −1 now gives the Bose-Einstein occupation number for ω > 0. The magnetic moments µ r have to be replaced by gµ B , where g is the spin gyromagnetic factor of the electron and µ B is the Bohr magneton, leading to γ/µ r = ℏ −1 , since the magnitude of the moments is described by the spin quantum number S in this case. Due to this different normalization, the decoupling coefficient reads α 0 = 1/ 2S 2 for quantum spins. On the left-hand side of Eq. (9), the product of the spin components S (1)

−qrS
(2) qs has to be replaced by the an- The expectation values of the sublattice magnetizations can be calculated from the Brillouin function as with X r = 2 arcoth (2Φ r ) using the definition of Φ r given above. Although the notations are different, it can be shown that the system of equations presented here is equivalent to Ref. [41] when the Dzyaloshinsky-Moriya and twoion anisotropy terms are set to zero. For T = 0, one obtains Φ r > 0 and S z r < S, indicating that the classical ground state is not the correct quantum ground state. Note that although the Brillouin function and in the classical limit the Langevin function also define the magnetization in meanfield theory, the argument of the functions differs from the mean-field model in Green's function theory; see Ref. [53] for a detailed discussion.
The main result of the Green's function formalism is the calculation of the frequencies of the two magnon modes ω + q and −ω − q in Eq. (10), and of the sublattice magnetizations in Eq. (14). These expressions allows us to calculate the temperature-dependent mesoscopic parameters via direct comparison to the excitation frequencies of the continuum model. Our theory is numerically validated in Sec. III, where the proposed analytical expressions are compared to numerical Monte Carlo simulations for a specific spin model, where the excitation frequencies can be given in a simpler form.

B. Effective temperature-dependent parameters
In the long-wavelength limit, the Hamiltonian in Eq. (1) can be approximated by the free-energy functional in d spatial dimensions. Here, the magnetization fields m r are required to be of unit length, being defined as M r m z r = µ r S z r /V c , where V c is the volume of the unit cell and M r is the saturation magnetization of the sublattice. The m subscript denotes effective mesoscopic parameters, while α, β and γ are Cartesian indices. The first line of Eq. (16) describes energy contributions from a spatially inhomogeneous magnetization. In the Dzyaloshinsky-Moriya term, the Lifshitz invariant L rs,αβ = 1/2 γ,δ ε αγδ m γ r ∂ β m δ s − m δ r ∂ β m γ s depends on the considered symmetry class [54]. The second line of Eq. (16) remains finite for homogeneous sublattice magnetizations, describing the energy contribution depending on the global orientation of the magnetization vectors with respect to the easy axis, to the external field, and to each other in the two sublattices.
The equation of motion in the continuum model reads which transforms into a form analogous to Eq. (6) in the local coordinate system and in Fourier space. Requiring that the excitation frequencies of the continuum model coincide with Eq. (10) of the atomistic model in the long-wavelength limit, the following temperature dependence is obtained for the parameters: where R ij = R i − R j . The comparison based on the magnon spectrum only provides information on the z component of the Dzyaloshinsky-Moriya interaction; the other components may be obtained by comparing the atomistic and continuum models for ground states oriented along different directions. The parameters Eq. (18)- (21) show similar trends to what has been calculated in single-sublattice ferromagnetic systems in Refs. [45][46][47]. The main contribution to the temperature dependence of all the parameters comes from S z r S z s , which corresponds to the mean-field or random-phase approximations. By considering a decoupling scheme different from the random-phase approximation, i.e., α 0 ̸ = 0, the temperature dependence of the micromagnetic parameters is corrected by taking spin correlation effects into account. The correlation corrections proportional to α 0 have a positive sign for the isotropic exchange and the Dzyaloshinsky-Moriya terms, which makes these terms decrease slower in magnitude with the temperature. Note that for only nearest-neighbor interactions, the relative correction to the isotropic exchange and Dzyaloshinsky-Moriya interactions turn out to be precisely the same, similar to what is observed in a single-sublattice ferromagnet in Ref. [46]. In contrast, the correlation corrections are negative for the anisotropy terms, indicating a faster decrease. In ferromagnets, this is known to correspond to the Callen-Callen law K ∼ ⟨S z ⟩ 3 for the temperature dependence of the uniaxial anisotropy [55] based on the first term in Eq. (20), and to a scaling exponent K ∼ ⟨S z ⟩ 2+ε slightly larger than 2 for the twoion anisotropy [47] in the second term. The last term in Eq. (20) gives a positive contribution to the anisotropy, meaning that the Dzyaloshinsky-Moriya interaction stabilizes collinear order in the presence of thermal fluctuations [46]. Equation (18) defines different mesoscopic exchange stiffness parameters for the intrasublattice coupling J AA m , J BB m and for the intersublattice coupling J AB m , J BA m . Together with the anisotropy, the value of these exchange stiffness parameters are relevant for the estimation of domain wall width δ w (T ) ∝ K(T )/J (T ). Equation (19) describes the temperature dependence of the mesoscopic intrasublattice and intersublattice Dzyaloshinsky-Moriya parameters, which are necessary for the estimation of the skyrmion radius [56,57]. The competition between the different contributions to the anisotropy term in Eq. (20) gives rise to fluctuation-driven spin reorientation transitions induced by the Dzyaloshinsky-Moriya interaction [58] and unusual exponents in the temperature dependence of the anisotropy parameter [47], similarly to what has been observed before in ferromagnetic systems.
As mentioned in Sec. II A, in the quantum caseS ir + = 1/4), which is consistent with the fact that the single-ion anisotropy just acts as a constant energy term. Indeed, this condition was one of the motivations behind choosing the value of the decoupling coefficient α 0 in Ref. [41].

C. Scaling exponents
To obtain a simpler formula for the temperature dependence of the parameters in the continuum model, we only keep the Heisenberg interactions which are typically the largest in magnitude, and calculate the expressions which are connected to Eqs. (18) and (21). The correlation functions are real in the absence of the Dzyaloshinsky-Moriya interaction, and we substitute them into Eq. (22) from Eq. (9) using Eq. (13), and approximateΓ q with the spin-wave HamiltonianH SW,q from Eq. (7) in the lowtemperature limit. This results in As mentioned above, the temperature dependence of the parameters of the continuum model is often expressed in terms of a power law of the magnetization. This is common practice partially because it is easier to implement numerically in a micromagnetic framework and partially because in non-equilibrium situations the value of the magnetization represents better the thermodynamic state of the system than the temperature of the heat bath. Both antiferromagnetic and ferrimagnetic systems may be characterized by the sublattice magnetizations S z A , S z B . As it is shown in Eqs. (18)-(23), the sublattice magnetizations are a more natural choice for expressing the effective parameters than the combinations µ A S z A ±µ B S z B resulting in the total and staggered magnetizations, respectively. The temperature T in Eq. (23) may also be expressed by either sublattice magnetization using the low-temperature expansion of Eq. (14), Equation (24) connects the two sublattice magnetizations to each other as well, meaning that either one can be used to express the effective parameters. To simplify the expressions further, we go to the classical limit, assume that all intrasublattice interaction terms are the same and there is no external magnetic field. In this case, it can be shown based on the definition Eq. (14) that S z A = S z B , meaning that the temperature dependence of all magnetizations is precisely the same if they are normalized to their zerotemperature value. This is the case for antiferromagnets with identical sublattices, but is also a good approximation for ferrimagnets with µ A ̸ = µ B if the intrasublattice interactions are negligible compared to the intersublattice ones. This follows from the fact that the self-consistency Eqs. (8), (13) and (14) do not depend explicitly on the magnetic moments. In the following we restrict our attention to this limit, and leave the case of different sublattice magnetizations observable in, e.g., ferrimagnets with a compensation point or in the quantum limit, to later studies.
As a specific example, we consider nearest-neighbor antiferromagnetic exchange J AB where S z is the magnetization on either sublattice and the correction to the intersublattice (AB, BA) and intrasublattice (AA, BB) exponents read with the geometrical factors  26) and (27). The inter-and intersublattice Heisenberg interactions are scaled from only antiferromagnetic coupling between the sublattices λ = 0 to two decoupled ferromagnetic sublattices λ = 1.
The ε values from Eqs. (26) and (27) are shown in Fig. 1 for the rock-salt structure, where the antiferromagnetically coupled sublattices together form a simple cubic lattice and each sublattice is an fcc lattice. For λ = 1 (ferromagnet), we numerically obtain an exponent close to the analytical value ε intra = 0.255, which was reported for the ferromagnetic fcc lattice in Ref. [59]. The intersublattice correction to the exponent vanishes in this limit as can be seen from Eq. (26), meaning that for weak intersublattice coupling parameters, the corresponding effective parameter closely follows a mean-field scaling. In the λ = 0 case (antiferromagnet), the intersublattice exponent also converges to the analytical value for a simple cubic lattice ε sc = 0.341, although the coupling between the sublattices is antiferromagnetic instead of ferromagnetic as in Ref. [59]. The correction to the intrasublattice exponent remains finite in this case, meaning that the influence of spin correlations on the temperature dependence of intrasublattice coupling can be observed even in this limit. Increasing λ leads to a decrease in both the intersublattice and intrasublattice ε values at first, because distributing the exchange interactions between nearest and next-nearest neighbors brings the system closer to a mean-field behavior. It was similarly found in Ref. [59] that the ε value is lower for ferromagnetic FePt where interactions with several neighbors were taken into account than for any nearest-neighbor cubic lattice. However, the scaling performed here demonstrates that this decreasing trend is reversed for λ values close to 1 in the intrasublattice term, while the intersublattice ε value vanishes as discussed above.
The agreement for the exponent corrections between the ferromagnetic and the antiferromagnetic case mentioned for λ = 0 in the simple cubic lattice is a general property of the model. If it is assumed that the antiferromagnetic alignment of the spins is realized in a system where all atoms together form a Bravais lattice, then one obtains γ AB q+Q = −γ AB q with Q the wave vector of the antiferro-magnetic ordering, making it possible to rewrite Eq. (26) for λ = 0 as where the summation now runs over the atomic or ferromagnetic Brillouin zone which is twice the size of the antiferromagnetic one. For infinite lattices where the summations can be replaced by integrals, the correction to the exponent is ε sc inter = 0.341 for the simple cubic and ε bcc inter = 0.282 for the body-centered cubic lattice [59], both of which can accommodate a two-sublattice ordering. Even for systems where all atoms together do not form a Bravais lattice (e.g., the honeycomb lattice), it can be derived that the expectation values and the correlation functions in the antiferromagnetically aligned model precisely coincide with those of the ferromagnetic model where the sign of all intersublattice coupling terms is reversed. Consequently, Eqs. (26) and (27) may also be used for ferromagnets containing both nearest-neighbor and next-nearest-neighbor interactions. The agreement between the ferromagnetic and antiferromagnetic cases essentially relies on the fact that in the classical limit, the self-consistency Eqs. (8), (13) and (14) do not depend on the magnon frequencies in Eq. (10), which are different between the ferromagnetic and the antiferromagnetic alignment. This is different in the quantum case, where Eq. (12) does depend on the frequencies as shown in Appendix A.
For weak Dzyaloshinsky-Moriya interaction, the same correction ε to the scaling exponent can be used. For twosite anisotropy between the same pairs of atoms as the exchange, the exponent is 2 + ε owing to the opposite sign of the correlation correction between Eqs. (18)- (19) and the second term in Eq. (20), respectively. For the on-site anisotropy term, the exponent is close to 3 as in the ferromagnetic case [55] since the on-site correlations are stronger than the two-site terms.
In two-dimensional systems, the sums in Eqs. (26) and (27) diverge for infinite lattice sizes, as is known from, e.g., the proof of the Mermin-Wagner theorem [60]. This implies that the exponents may only be calculated if a finite anisotropy is taken into account, in which case they have to be evaluated numerically. Since the correlation corrections are expected to be enhanced in low-dimensional systems, this procedure is carried out and compared to numerical simulations in Sec. III.

III. SIMULATIONS
To probe the accuracy of the analytical method described in Sec. II, its predictions will be compared to the numerical simulations based on the Hamiltonian Eq. (1). While the magnetization and the static correlation functions may be directly determined from averaging over spin configurations from the different simulation steps, the frequencies required for determining the temperature dependence of the parameters in the continuum model are more difficult to access. Equations (9)-(13) establish the relations between the expectation values and the frequencies. They may be reformulated as The product of the frequencies ω prod = ω + q ω − q is given by Eq. (30), which also yields the sum ω sum = ω + q + ω − q from Eq. (31). The individual frequencies may be calculated as Note that in Eqs. (30) and (31), the correlation functions are given in the global coordinate system for easier implementation in the simulations. Since these equations establish a connection between the eigenfrequencies, the correlation functions and the temperature, they may be considered as a form of the equipartition theorem. Although Eqs. (30) and (31) were determined from the Green's function formalism, they do not depend on the explicit form of the decoupling α 0 , only on the assumption that the spectral density is concentrated in single-particle excitations. Therefore, substituting the expectation values obtained from the simulations into Eq. (30) and (31) enables the calculation of the frequencies of the simulated system. Furthermore, this method allows for determining the frequencies based on Monte Carlo simulations, which accurately describe thermal equilibrium properties but do not provide direct access to the real-time dynamics of the system. The simulated model system is illustrated in Fig. 2. It consists of a square lattice with equivalent sublattices µ A = µ B = µ S , only considering nearest-neighbor intersublattice Heisenberg exchange −J AB = J > 0 and Dzyaloshinsky-Moriya interactions of magnitude D AB = D, with the Dzyaloshinsky-Moriya vectors being perpendicular to the lattice vectors connecting the neighbours following a C 4v symmetry. The easy axis K A = K B = K was assumed to lie along one of the nearest-neighbour directions, which enables the investigation of the Dzyaloshinsky-Moriya vectors parallel to the z direction on the spin-wave spectrum. The external magnetic field was set to zero. We performed Monte Carlo simulations on a 64 × 64 lattice using the single-spin Metropolis algorithm where the trial spin direction is chosen uniformly on the surface of the unit sphere. The lattice was equilibrated for 2 · 10 5 Monte Carlo steps at each temperature, then the expectation values were calculated from data obtained over 10 8 Monte Carlo steps. To further improve the accuracy, 50 independent simulations were averaged in the end.
Due to the symmetry of the sublattices, we obtain S z A = S z B = ⟨S z ⟩, which also coincides with the dimensionless staggered magnetization n. The simulated and calculated values of n are compared in Fig. 3, demonstrating good agreement. Including the Dzyaloshinsky-Moriya interaction decreases the staggered magnetization at a fixed temperature. The critical temperature of the system is around k B T c ≈ 0.84J from Green's function theory. Note that ⟨S z ⟩ is not possible to calculate accurately at temperatures close to T c , since due to the relatively small system size and the long simulation length the system starts to switch between the +z and −z directions.
The spin-wave spectrum at finite temperature was calculated based on Eq. (30), since the symmetry of the sublattices implies ω + q = −ω − q , and both sides of Eq. (31) vanish. For the considered system, the two branches of the spinwave dispersion relation are given by The spectrum is illustrated in Fig. 4. The Dzyaloshinsky-Moriya interaction lifts the degeneracy of the two branches and shifts the minimum of the spectrum away from q x = 0. The anisotropy opens a gap in the spectrum, which is exchange enhanced compared to the ferromagnetic case:   the parameters J , D = |D ij |, and K are defined as  Eqs. (19) and (20) appear due to the sign change in J and the antiferromagnetic alignment of the sublattices, respectively. Figure 4 supports the high accuracy of the Green's function formalism up to intermediate temperature values of k B T = 0.40J. The scaling of the parameters in the magnon spectrum with the staggered magnetization n is shown in Fig. 5, displaying a power-law behavior as discussed in Secs. II B and II C. The explicit temperature dependence is illustrated in Fig. 7 in Appendix A for comparison. These results confirm the reliability of the Green's function method in predicting the simulation results. The Dzyaloshinsky-Moriya interaction D decreases slower in temperature than the anisotropy K, as discussed in Sec. II B for the general case. The temperature dependence of the Heisenberg term J is identical to that of the Dzyaloshinsky-Moriya interaction in Green's function theory and agrees with it in the simulations within error bars; therefore, it is omitted from the figure. Based on a fit to the simulation data, the scaling exponent is 1.58 for the Dzyaloshinsky-Moriya interaction, decreased by the correction ε inter = 0.42 compared to the uncorrelated value. The scaling exponent agrees with the value of 1.54 − 1.57 obtained for the ferromagnetic case in Ref. [46]. For the anisotropy, an exponent of 3.03 is obtained without Dzyaloshinsky-Moriya interaction, rather close to the well-known power law predicting an exponent of 3 [55]. In the presence of the Dzyaloshinsky-Moriya interaction, the exponent is slightly reduced to 2.92, i.e., there is an additional positive contribution to the temperature dependence of the uniaxial anisotropy due to the Dzyaloshinsky-Moriya interaction.
The accuracy of the decoupling scheme may be better visualized after subtracting n 2 from the normalized parameters, leaving only the correlation corrections shown in Fig. 6. Note that in the random-phase approximation obtained for α 0 = 0, the curves would be zero as indicated by the dashed line in the figure. Comparing Figs. 5 and 6, it is clear that the correlation corrections are not negligible, contributing around 10% of the total value of the Dzyaloshinsky-Moriya interaction and around 50% of the total anisotropy at the highest simulated temperatures. As mentioned earlier, for D = 0 the correction to the anisotropy will be n 3 − n 2 , i.e., it results in the Callen-Callen power law [55]. The corrections are positive for the exchange and negative for the anisotropy terms as mentioned above, leading to increased and decreased effective exponents, respectively. While in this plot the deviations between Green's function theory and the simulations become apparent, even for the anisotropy terms the analytical description reproduces about 2/3 of the corrections observed in the simulations. The accuracy appears to be higher for the Dzyaloshinsky-Moriya interaction itself and its correction to the anisotropy (difference of the orange and yellow lines).

IV. CONCLUSION
We applied Green's function theory to calculate the magnon frequencies in two-sublattice antiferromagnetically aligned systems, and to determine the temperature dependence of the interaction parameters in the magnon spectrum. We found that transversal spin correlations stabilize the Heisenberg and Dzyaloshinsky-Moriya exchange interactions against thermal fluctuations, but induce a faster decay of the anisotropy terms with the temperature. The Dzyaloshinsky-Moriya interaction also contributes to the uniaxial anisotropy term via the spin correlations, increasing its value at finite temperature in contrast to the typical decrease. We obtained good agreement between the predictions of the theory and Monte Carlo simulations performed on a square lattice, where the correlations play a pronounced role due to the reduced dimensionality.
Remarkably, these observations do not simply qualitatively agree with previous calculations for ferromagnets [45][46][47], but a mathematical correspondence can also be established. The self-consistency equations may be exactly transformed into each other in the classical limit when reversing the magnetization direction on one sublattice simultaneously with the sign of all intersublattice coupling terms. If the intrasublattice interactions are identical, the sublattice magnetizations and consequently the total and staggered magnetizations show precisely the same temperature dependence, even if the magnetic moments on the sublattices are different. Therefore, the scaling relations of the inter-and intrasublattice coupling terms discussed here can also be applied to ferromagnets with nearest-neighbor and next-nearest-neighbor interactions.
The calculated temperature dependence of the parameters are fundamental for the development of multiscale models connecting first-principles spin-model parameters to finite-temperature mesoscopic computational approaches, such as micromagnetism or the Landau-Lifshitz-Bloch equation. Most of the multiscale approaches proposed so far rely on an intermediate step based on classical spinmodel simulations, which could be replaced by the considerably more efficient semi-analytical expression presented here. Multiscale methods would be able to access the dynamics of and the phase transitions in antiferromagnetically aligned systems, for example for a realistic and computationally efficient description of all-optical ultrafast switching processes in ferrimagnets [61] or of magnetic domain wall motion in antiferromagnets [62].
Deviations in the equilibrium parameters from singlesublattice ferromagnets are expected to be observed in systems where the intrasublattice terms are not equivalent, such as ferrimagnets with a compensation point, or particularly when quantum effects are taken into account. Validating the predictions of Green's function theory in the quantum limit would require comparisons with classical spinmodel simulations augmented by a semi-quantum thermostat [63] or with renormalized heat-bath temperatures [64], or to quantum spin-model simulations based on quantum Monte Carlo [65] or tensor-product states [66]. The multiscale quantum approach would be completed by using the calculated temperature-dependent parameters in the quantum version of the Landau-Lishitz-Bloch equation [67]. by MCIN/AEI/10.13039/501100011033 and by "ERDF A way of making Europe" and "ESF Investing in your future".

Appendix A: Derivation of the self-consistency equations
The dynamics of the spin system is generated by the Poisson brackets in the classical limit, where i, j are lattice indices, r, s are sublattice indices and α, β, γ are Cartesian indices as in the main text. In the quantum case, the Poisson brackets {S α ir , S β js } have to be replaced by the commutators as −iℏ −1 [S α ir , S β js ], and γ/µ r has to be replaced by ℏ −1 in Eq. (A1) and the following expressions. Note that in the quantum limit we introduced a sign change compared to the conventional commutation relations, since the S α ir operators represent the dimensionless magnetic moments of electrons which are antiparallel to the angular momenta.
The time-dependent Green's function is defined as where θ (t) is the Heaviside function, ⟨⟩ denotes averaging in thermal equilibrium, and u is a real parameter. Here, α and β will primarily denote the ladder operator indices + and −, but unless they are explicitly specified the expressions are also valid for the Cartesian indices. The Green's function satisfies the equation of motion The second term on the right-hand side of Eq. (A3) introduces higher-order Green's functions, which will be decoupled using Here we took advantage of the rotational symmetry of the system: only such expectation values are considered in the decoupling which are rotationally invariant around the z axis, namely ⟨S z ir ⟩ ≡ ⟨S z r ⟩ which is the same at all sites in the sublattice due to the translational invariance of the ground state, and S + js S − ir that is replaced by half of the anticommutator in the quantum limit. The decoupling coefficients are chosen as α s+r− = ⟨S z r ⟩ α sr 0 , with α sr 0 = 1/2 in the classical and α sr 0 = 1/ (2S r S s ) in the quantum case. This choice of the decoupling parameters will be motivated later.
We introduce the transformed coordinate system with the sublattice spins pointing along the local z direction as discussed in the main text, and perform temporal and spatial Fourier transformation via ∂ t → −iω and Following the decoupling and the Fourier transformation, Eq. (A3) reads and theΓ rs q coefficients are introduced in Eq. (8) in the main text, with the components given bỹ We again used the short-hand notationS (2) qr ∈ S − qA ,S + qB . As in linear spin-wave theory in Eq. (6), the corresponding equations for theS (1) qr ∈ S + qA ,S − qB spin components decouple from Eq. (A7). The equations for S (1) qr yield the other branch of the spin-wave dispersion relation shown in Fig. 4, but they need not be solved separately since they are connected to the ω ± q frequencies by particle-hole symmetry. Equation (A7) is solved as where the inverse matrix has poles at the real frequencies ω ± q given in Eq. (10). These poles may be used to evaluate the correlation functions via the spectral theorem (cf. for u = 0. Here Eq. (9) was used to introduce the Φ rs q quantities, andΘ A−;A+ (u = 0) = 2 ⟨S z A ⟩ andΘ B+;B− (u = 0) = −2 ⟨S z B ⟩ were substituted based on the Poisson brackets. In the classical limit, these equations simplify to with detΓ q =Γ AA qΓ BB q −Γ AB qΓ BA q . These equations are summarized in Eq. (13). Note that Φ rs q is symmetrized by the anticommutator in the quantum case in the convention used here, and the components are given by Eqs. (A15)-(A18). The sublattice-diagonal part of Eq. (A14) also yields a differential equation in u, which is of the same form as the one described in Ref. [40] in the quantum case and in Ref. [46] in the classical limit; the solution of this equation with the appropriate boundary conditions gives the final equation (14) or (15) [41]. This constraint is satisfied by the choice used in the main text and described after Eq. (A4). However, it also allows for using different decoupling schemes for the intrasublattice and the intersublattice terms. Figure 7. Dependence of the effective interaction parameters on the temperature. All quantities are normalized to their zero-temperature value. Results of the numerical simulations (symbols) are compared to Green's function theory calculations (lines). Simulation data were obtained by fitting the functions in Eqs. (33) and (34) to the simulated frequencies; error bars denote the error of this fit. The atomistic model parameters are K = 0.1J and D = 0.2J for the blue and orange curves, D = 0.0J for the yellow curves.