An effective spin model on the honeycomb lattice for the description of magnetic properties in two-dimensional Fe$_3$GeTe$_2$

Fe$_3$GeTe$_2$ attracts significant attention due to technological perspectives of realizing room temperature ferromagnetism in two-dimensional materials. Here we show that due to structural peculiarities of the Fe$_3$GeTe$_2$ monolayer, short distance between the neighboring iron atoms induces a strong exchange coupling. This strong coupling allows us to consider them as an effective cluster with a magnetic moment $\sim$5 $\mu_B$, giving rise to a simplified spin model on a bipartite honeycomb lattice with the reduced number of long-range interactions. The simplified model perfectly reproduces the results of the conventional spin model, but allows for a more tractable description of the magnetic properties of Fe$_3$GeTe$_2$, which is important, e.g., for large-scale simulations. Also, we discuss the role of biaxial strain in the stabilization of ferromagnetic ordering in Fe$_3$GeTe$_2$.


I. INTRODUCTION
The discovery of graphene in 2004 by Novoselov et al. [1] attracted enormous attention to the twodimensional (2D) materials due to their unique properties and possibility to control them by means of gate voltage, chemical doping, or strain [2][3][4].Since then the search and description of novel 2D materials is an actively developing research direction in material science.The discovery of 2D magnetic materials [5][6][7][8] has further enhanced the interest to the field of 2D materials, opening new ways to control magnetism in two dimensions, which is prospective for applications.One of the prominent example among 2D magnets is CrI 3 , which was successfully exfoliated from the bulk crystal [8] demonstrating the great advantage to control magnetic properties by electrostatic doping [9] or hydrostatic pressure [10].The measured Curie temperature of this material is about 45 K, slightly lower than in the bulk crystal.The other typical representative of the 2D magnet family is Fe 3 GeTe 2 (FGT) [11].The Curie temperature of FGT in the bulk phase is considerably higher T C ∼ 220 K [11][12][13][14], making this material more perspective from the practical point of view.Despite a decrease of the Curie temperature down to T C ∼ 130 K [15] in the monolayer limit, the corresponding value is still three times larger compared to CrI 3 .Moreover, the gate control of FGT allows one to enhance T C giving the opportunity to realize room-temperature ferromagnetism in 2D [11].Other studies report large anomalous Hall effect [16], anomalous Nernst effect [17] and magnetic stability [18] making FGT a promising candidate for spintronics, caloritronics and other applications [11,19].* a.rudenko@science.ru.nlAnother interesting aspect of FGT is the absence of inversion symmetry [20], which gives rise to nontrivial physics and complex spin textures such as magnetic skyrmions and spin-spirals [21][22][23], presumably emerging due to the Dzyaloshinskii-Moriya interaction (DMI) [24].On the other hand, some authors introduced four-spin interaction to explain the stabilization of complex magnetic patterns [25,26].Therefore, previously proposed conventional magnetic models need to be treated carefully, considering peculiarities of the crystal structure.In particular, there are two valence types of Fe atoms in FGT: Fe 2+  II and Fe 3+ I , where Fe I atoms with comparably short distance (∼2.47 Å) form a dimer-like structure (Fig. 1).The formation of such dimers is inherent for some insulating systems [27][28][29], which challenges the description of their magnetic properties.For instance, the short distance between the iron atoms in the dimer increases the overlap of wave functions, enhancing the hybridization and exchange interactions between the magnetic atoms.At the same time, strong interaction between the atoms leads to an ambiguity in the definition of local magnetic moments, challenging the description of magnetic properties within the spin models.
In this paper, we propose an alternative description for the magnetic properties of Fe 3 GeTe 2 monolayer, taking advantage of its structural peculiarities.We demonstrate that the iron dimer can be considered as an effective cluster, such that the spin model of FGT reduces to a bipartite honeycomb lattice (Fig. 1).The proposed model is further justified by means of the Monte Carlo and spin dynamics simulations, which perfectly reproduce the results of the conventional spin model in the relevant energy region at low computational cost.Additionally, we discuss the role of in-plane biaxial strain in the stabilization of magnetic order.The rest of this paper is organized as follows.Sec.II briefly describes numerical methods used in this paper.In Sec.III, we present our main results, including the alternative magnetic model, and thermodynamics of FGT.The effect of biaxial strain is also discussed.In Sec.IV, we summarize our results and conclude the paper.

II. METHODS
To study the magnetic properties of monolayer Fe 3 GeTe 2 , we consider the following spin model: Here, J ij stands for the isotropic exchange interactions between iron ions, and A is the single-ion anisotropy.
If A > 0, Eq. ( 1) describes an easy-axis magnet with the ground-state magnetization perpendicular to the 2D plane.For simplicity reason we map our model to unified spin S = 1.Magnetic exchange interactions were calculated using the local force theorem approach [30,31] where m, m ′ , n, n ′ are orbital quantum numbers, ∆ mm is the intra-orbital spin-splitting energy and Ĝ(ϵ) = 1/(ϵ− Ĥ) is single-particle Green's functions.In case of conventional model orbital indices run within orbitals of single site, while for conventional model they run within orbitals of two sites (cluster).In contrast to the mapping procedure of total energies of collinear spin configurations [32,33], this approach allows one to estimate long-range exchange couplings without using large unit cells.We do not consider DMI in the present study, as we are mostly interested in the thermodynamic properties of FGT as well as in stability of the FM order, where intersite anisotropic interactions play a minor role.Density-functional (DFT) band-structure calculations were performed within the Perdew-Burke-Ernzerhof (PBE) exchange-correlation functional [34] as implemented in Vienna ab initio simulation package (vasp) [35,36].Additionally, we also performed DFT+U calculations with simplified rotationally invariant scheme [37] using the effective parameter U = 4 eV [38], and further use these results for comparative analysis.In all these calculations, we set the energy cutoff of the plane-wave basis to 400 eV, the energy convergence criteria to 10 −8 eV and use a (18×18×1) Γ centered grid for the Brillouin zone integration.The experimental crystal structure of bulk Fe 3 GeTe 2 was used [39], where a vacuum space more than 16 Å between monolayer replicas in the vertical z direction was introduced.The positions of atoms were allowed to relax until all the residual force components on each atom were less than 10 −3 eV/ Å.From the calculated electronic structure, maximally localized Wannier functions [40] were constructed using the wannier90 package [41] projected onto the 3d and 5p states of iron and tellurium, respectively [42].This ensures that the constructed Wannier functions for both spin channels are well localized, which provides an atomic-like basis for the application of Eq. ( 2), and for the analysis of magnetic properties.
Classical Monte Carlo (MC) simulations of the constructed spin Hamiltonian [Eq.( 1)] as well as spin-wave dispersion calculations were carried out using the Uppsala Atomistic Spin Dynamics (UppASD) package [43,44].All the results for both conventional and alternative models are obtained using the exact lattice structures presented in Fig. 1.At each temperature, we perform 20 000 MC steps for thermalization and 150 000 for measurements.Magnetization and specific heat curves are calculated for a lattice containing 90 × 90 unit cells with periodic boundary conditions.Importantly, to obtain clear peaks in the specific heat, we averaged the graphs over ten different MC runs and additionally smoothed the resulting curve.Spin-spin correlators ⟨S 1 S 2 ⟩ for nearest neighbor Fe I and Fe II atoms were calculated for each considered temperatures averaging the results over 100 MC configurations.
Spin-wave dispersion is calculated within the framework of the classical spin dynamics approach by solving the Landau-Lifshitz-Gilbert (LLG) equation.
where γ is the gyromagnetic ratio, α is the damping parameter and b i (t) is a stochastic magnetic field with a Gaussian distribution arising from the thermal fluctuations.For this, we calculate the space-and timedisplaced correlation functions where ⟨. . .⟩ denotes the ensemble average and k is the Cartesian component.The corresponding Fourier transform gives the dynamical structure factor with N being the total number of magnetic atoms.This quantity can be probed in neutron scattering experiments of bulk systems [45].The energy dispersion for spin waves presented in the simulated system is described by the positions of the peaks in χ k (q, ω) [43,46].Thus, the resulting values at each q vector are convoluted with a Gaussian filter and then normalized to make peaks at any q positions visible regardless of their relative intensity [47].
All the results presented in this paper are obtained for a lattice of (450 × 450) unit cells with periodic boundary conditions.Adiabatic magnon frequencies ω qν were calculated via diagonalizing the spin-wave Hamiltonian, defined for magnetic sublatices µ and ν in the unit cell [48]: ) where J µν (q) are the Fourier transform of the exchange interactions.The magnon frequencies allow us to estimate the Curie temperature within the random phase approximation (RPA) as [48,49]: (5)

A. Magnetic models
The resulting isotropic exchange couplings are represented in Supplemental Material [42].The nearestneighbor FM interaction J 1 = −113 meV is three times larger than J 2 and an order of magnitude stronger than AFM J 3 .The values of these exchange couplings were additionally checked by using the total energies of collinear structures [32,33], which yields J 1 = −113.9meV and J I 3 = 10.9 meV.The obtained parameters are also in reasonable agreement with recently reported values in Ref. 50, where J 1 = −146.3meV, J 2 = −36.2meV and J I 3 /J II 3 = 10.9/4.5 meV.Short-range interactions J 2 and J 3 originate from the ligand-mediated superexchange mechanism.Due to a metallic character [51], the exchange couplings demonstrate an oscillating RKKY-like behavior at increasing distances between the iron atoms (Fig. 2).These long-range interactions are essential for the simulation of magnetic systems by means of the spin models [52].The resulting magnetic moments in these calculations are m I = 2.43 µ B and m II = 1.48 µ B .
Strong coupling between nearest neighbour iron atoms allows us to propose the new model, where two Fe I atoms in the dimer are to be considered as an effective cluster with ten 3d orbitals.Within this alternative model the exchange couplings can be reformulated in terms of the interactions between these effective clusters and Fe II atoms (Fig. 1).The proposed alternative model, which has the form of a bipartite honeycomb lattice, is more tractable compared to the original 3-site model as it allows to reduce the number of interactions needed for simulations, and helps to eliminate the strong coupling J 1 within the Fe I -Fe I cluster.Further we will use notation of exchange interactions for alternative model with tilde Ji .Below, we perform a systematic comparison between the magnetic properties simulated using the two different spin models, and discuss the performance of the alternative model in more detail.

B. Model simulations
Monte Carlo results for both spin models are represented in Fig. 3.In these simulations, we use the single-ion anisotropy parameter A = 0.35 meV/Fe ion, which was calculated from the total energy difference between the configurations with the magnetic moments oriented in-plane and out-of-plane directions taking spinorbit coupling into account [38,50].Without loss of accuracy, we consider the exchange interactions around the iron atoms within the radius r < 13 Å (Fig. 2).From Fig. 3, one can see that the calculated magnetization curves for both spin models perfectly match and drop down at T ∼ 340 K.At the same time, the estimated critical temperature via the adiabatic magnon spectra gives T C = 336 K and 328 K for the conventional and alternative models, respectively.
Figure 5 shows the adiabatic magnon dispersion relation calculated for the two models.The resulting low energy (acoustic and optical) magnon branches are in good agreement between the models.The spin-wave stiffness constants are also close to each other and equal to D ≈ 445 meV• Å2 and D ≈ 408 meV• Å2 for the conventional and alternative models, respectively.The magnon gap for both models is 1.36 meV since it is independent of the exchange interactions.The most prominent feature in the magnon spectrum of the conventional model is the presence of a nearly flat branch at high energies, which is mainly originated from the nearest-neighbor J 1 interaction between iron atoms.In the alternative model, the flat branch is absent (Fig. 5) because the two Fe I atoms are replaced by a single cluster.The energies of the corresponding flat branch are extremely high for magnons.This suggests that these states are strongly damped by other spin (or even electron) excitations in the system.
To examine whether one can expect to observe this flat branch in the experiments, we calculate the magnon spectral function for the excited spin waves using spin dynamics simulations, which is a theoretical analog of the inelastic neutron scattering.Our calculations do not take coupling with the Stoner excitations into account, but allows us to elucidate on the magnon damping resulting from the localized spin excitations.As can be seen from Fig. 5, both models behave in a similar way, reproducing the adiabatic magnon dispersion with decaying intensity toward high energies.Even for the conventional model, the high energy flat branch is strongly damped.The obtained behavior is not unusual.For example, the authors of Ref. 53 demonstrate that the intensities of the optical modes near the center of the Brillouin zone are suppressed because of the specific dynamical properties of the Heisenberg model, leading to strong damping of the standing modes observed in the adiabatic spin-wave calculations.We suppose that the absence of the flat branch in the spin dynamics simulation of Fe 3 GeTe 2 is in line with this scenario.Another reason why high-energy optical mode is absent may be associated with the strong short range magnetic order between Fe I atoms, which is demonstrated by the corresponding spin-spin correlation functions in Fig. 4. Strong interaction within the Fe I -Fe I pair suppresses excitations of individual spins.
The obtained Curie temperature overestimates the experimental T C ∼ 200 K [11][12][13][14], which is likely due to overestimation of the exchange interactions.For 3d metals it is instructive to apply a Hubbard U correction within DFT+U method.However, our results show that this method results in even further overestimation of exchange interactions [42] and T C .This discrepancy could be attributed to a limited applicability of the localized spin models to FGT, which neglect coupling to the electronic subsystem, especially important for metals.Partially itinerant character of FGT is also suggested by the non-integer magnetic moments on Fe atoms.Previous works demonstrated a significant renormalization of the electron spectrum [16] as well as site-dependent correlations [54], as follows from dynamical mean-field studies.An explicit consideration of the electron subsystem as well as electron correlation effects indeed contribute to the reduction of the Curie temperature in Fe 3 GeTe 2 , as it has been recently demonstrated [55].At the same time, the description of magnetism within the localized spin models still remain possible with renormalized exchange interactions.This, therefore, is not expected to limit the applicability of the alternative model we propose in our study.The constructed models allow us to study the magnetic properties of FGT under in-plane biaxial strain.For this purpose, we vary the lattice constant a from −5% (compressive) to 5% (tensile) and recalculated the isotropic exchange couplings.The resulting magnetic moments of Fe II remain nearly the same (1.50 µ B ), while the moments of Fe I increase up to 2.6 µ B in case of the tensile strain, and reduce down to 1.53 µ B under compression.This is in agreement with the results reported previously [56], also showing that tensile strain significantly enhances the FM stability and increases T C [57].
For the conventional model, the nearest-neighbor FM interaction J 1 increases its absolute value under pressure (Fig. 6), which, however, does not affect the Curie temperature keeping in mind the adiabatic magnon spectra.The next-nearest-neighbor FM interaction J 2 demon- strates only a moderate change under strain.More importantly, the third-nearest-neighbor interaction between Fe I atoms, i.e.J I 3 changes its sign from AFM to FM at ∆a ≈ 3%.A similar tendency can be seen for JI 2 within the alternative model in Fig. 6(b).Such a behavior reduces magnetic frustration, making the FM order more preferable under tensile strain.This is further demonstrated by the Monte Carlo simulations as well as by calculations within RPA, both leading to an increase of the Curie temperature at ∆a > 0 (see Fig. 2 and Table I).
At the same time, compressive strain suggests a destabilization of the FM order due to the presence of AFM exchange interactions.Indeed, from the adiabatic magnon spectra (Fig. 7) one can see that the acoustic mode becomes imaginary at finite wave vectors, indicating instability of the FM ground state in compressed Fe 3 GeTe 2 .

IV. CONCLUSION
In this work, we present a systematic description of the magnetic properties of monolayer Fe 3 GeTe 2 using a combination of spin Hamiltonians and ab initio calculations.A strong coupling between the nearest iron atoms motivates us to consider them as an effective cluster with a magnetic moment of ∼5 µ B , giving rise to a simplified spin model on a bipartite honeycomb lattice.This lattice is in many respects a more tractable and fundamentally important model in physics, which permits analytical calculations and might be important for the development or ∆a = 0% ∆a = 3% ∆a = 5% ∆a = -3% ∆a = -5% FIG. 7. Evolution of the adiabatic magnon spectra calculated for the conventional [Fig.1(a)] and alternative [Fig.1(b)] spin models of Fe3GeTe2 under strain.
extensions of many-body theories related to collective excitations in 2D honeycomb magnets [58].The alternative spin model perfectly reproduces the results of the conventional three-site model in the moderate-energy region, which is demonstrated by simulating the magnon spectra as well as thermodynamical properties.The alternative model allows us to reduce the number of long-range interactions needed for modeling of metallic systems, which is important, e.g., for large-scale simulations.We also find that the stability of the FM ordering in monolayer Fe 3 GeTe 2 can be enhanced under tensile strain.

VI. ISOTROPIC EXCHANGE INTERACTIONS
In Table I we present the isotropic exchange interactions between iron atoms (or clusters) in Fe 3 GeTe 2 , calculated at ambient pressure.
FIG. 1. (Top) Crystal structure of Fe3GeTe2 monolayer: Side and top views.FeI and FeII denote two inequivalent iron atoms.(Bottom) Panels (a) and (b) demonstrate the conventional and alternative magnetic models.

FIG. 2 .
FIG.2.Evolution of exchange couplings as a function of the distance between FeI atoms shown for the two different spin models.

FIG. 5 .
FIG. 5. Comparison of the magnon spectral function intensity map for conventional model (a) (left) and alternative model (b) (right), obtained via spin dynamics calculations at T = 5 K and damping constant α = 10 −3 .Red lines correspond to the adiabatic magnon spectra.

3 FIG. 6 .
FIG. 6. Evolution of the nearest-neighbor and next-nearestneighbor exchange interactions in the conventional [Fig.1(a)] and alternative [Fig.1(b)] spin models for various values of strain.For details see the crystal structure, corresponding exchange interactions are given in the same color scheme.Note that short distance interaction between FeI atoms (J1 shown with blue solid line) in Model (a) is absent in Model (b).

TABLE I .
Comparison of the Curie temperature TC (in K) calculated for Fe3GeTe2 via MC simulations and RPA approach for the conventional and alternative magnetic models of Fe3GeTe2 under strain.Values in brackets correspond to results of simulations using DFT+U parameters with U = 4 eV.

TABLE I .
Main exchange coupling (in meV) of Fe3GeTe2 within conventional and alternative magnetic model calculated using DFT and DFT+U (U = 4 eV) electronic structure.