From electronic structure to magnetism and skyrmions

Solid state theory, density functional theory and its generalizations for correlated systems together with numerical simulations on supercomputers allow nowadays to model magnetic systems realistically and in detail and can be even used to predict new materials, paving the way for more rapid material development for applications in energy storage and conversion, information technologies, sensors, actuators etc. Modeling magnets on different length scales (between a few Å ngström and several micrometers) requires, however, approaches with very different mathematical formulations. Parameters defining the material in each formulation can be determined either by fitting experimental data or from theoretical calculations and there exists a well-established approach for obtaining model parameters for each length scale using the information from the smaller length scale. In this review, this approach will be explained step-by-step in textbook style with examples of successful scale-bridging modeling of different classes of magnetic materials from the research literature as well as based on results newly obtained for this review.


Introduction
Materials science in the second half of the twentieth century has been marked by the development of theoretical methods and application of supercomputers for simulating and studying properties of real materials.One of the major breakthroughs of theory was the density functional theory (DFT) suggested by Hohenberg and Kohn [1] and the subsequent creation of the special ansatz by Kohn and Sham [2].DFT allows to describe, to a good approximation, the complex correlated behavior of electrons which determine most of the material properties.The Kohn-Sham ansatz suggests a practical way of doing so when a sufficiently accurate form of electronic correlation energy is arXiv:2310.08628v1[cond-mat.mtrl-sci]12 Oct 2023 constructed.These developments created a flurry of activities in theoretical materials science aimed at finding better approximations for electronic correlations and more efficient numerical methods for solving the DFT equations.During the decades after the seminal Hohenberg-Kohn and Kohn-Sham papers, DFT has proved to be an efficient and reliable theory for describing realistic material models and there are many success stories where DFT calculations played a major role in understanding a certain physical phenomena or designing new functional materials (for examples see [3][4][5]).
At this level of theory, all the details of the crystal structure are fully taken into account and the relevant quantities obtained in calculations are the electronic wavefunction Ψ({r i }), density n(r) and energy spectrum ε n (k) as a function of the wavevector k.The computational complexity of DFT calculations scales, at least, as N 3 where N is the number of electrons/atoms in the unit cell and DFT simulations are usually restricted to length scales of up to several nanometers.Numerical methods optimized for calculating electronic structure in real space [6] have been developed to tackle the DFT equations for a few thousand [7] or even several hundred thousand atoms [8], but this does not allow to go much beyond the aforementioned length scale.
Information about the nanometer-scale behavior is, however, not enough for understanding the magnetic phenomena in real materials, since new behavior can emerge on larger length scales.To circumvent the difficulty of DFT and quantum-mechanical methods, in general, to describe large groups of atoms, it is customary to construct effective models, which are supposed to describe a subset of materials properties.In the context of magnetism, a spin model with Heisenberg J ij and Dzyaloshinskii-Moriya (DM) ⃗ D ij interactions and on-site anisotropy (for this example, uniaxial anisotropy of magnitude K U ) is the most famous example: Instead of describing the whole electronic energy spectrum and all the electronic degrees of freedom, this generalized Heisenberg model and models similar to it focus instead on the low-energy part of the spectrum, close to the ground state (Fig. 1a).The electronic spin density is assumed to be well localized around each atom (i), so that its norm is constant and only the direction of spins (indicated by unit vectors ⃗ S i ) changes as a result of electron dynamics.In that case, expression (1) describes the energy of excited states where spins on different atoms deviate from the magnetic ground state, which includes spin waves and single-atom excitations.The advantage of the Heisenberg model is that parameters J ij contain a large part of complexity of electronic properties and the model itself is easier to solve than the equations of the density functional theory.This allows to perform simulations for up to 10 8 spins, which covers the length scales up to several hundred nanometers relevant for many interesting magnetic states, such as spin spirals, domain walls, skyrmions, vortices etc. [9,10] (more details in Sec. 5).
Further increase of the length scale of accessible magnetic phenomena is possible by means of micromagnetic simulations where the discrete atomistic character of material is neglected and instead the continuous description of magnetic properties is applied.The Figure 1.a) Main idea of the effective spin model.Low-energy part of the full electronic spectrum is used to construct a model that describes excitations from the magnetic ground state.b) Illustration of the micromagnetic approach where the material is described by a continuous magnetization ⃗ m(⃗ r), instead of discrete lattice of spins ⃗ S i sitting at positions ⃗ r i .The system is subdivided into small regions ∼ 1 nm with the local magnetization (shown by arrows) representing the average of atomic spins in that volume.
main idea of the micromagnetic approach is illustrated in Fig. 1b, where a continuous function ⃗ m(⃗ r) describes the local magnetization in different regions of the system and replaces the information about spin vectors ⃗ S i on different atomic sites at positions ⃗ r i , while ideally ⃗ m(⃗ r i ) = ⃗ S i should be satisfied.This approach works well when the magnetization varies on the length scale (l) much larger than the characteristic atomic distance (a), so that the nearest-neighbor spins have similar orientations.Magnetic ground state can be determined by minimizing the micromagnetic functional E = ε(⃗ r) dV with the energy density ε ≡ ε(⃗ r): The two terms in this functional directly correspond to the Heisenberg and DM interactions in the spin model (1) and the parameters, spin stiffness A and Dzyaloshinskii-Moriya matrix D, can be determined from the interatomic J ij and ⃗ D ij interactions, while the uniaxial anisotropy has the same form but different units in the atomistic and micromagnetic cases.In the following (Sec.4), we will discuss in detail how the micromagnetic functional (2) is derived.
Both in atomistic and micromagnetic approaches described by models (1) and ( 2), the magnetization dynamics can be studied by calculating the effective field ⃗ B eff acting on each spin or micromagnetic element and integrating the Landau-Lifshitz-Gilbert equations (further details in Sec. 5).The micromagnetic method allows to approach the length scales of up to a few micrometers, which is necessary, for example, for studying wide domain structures and relatively large groups of topological magnetic objects, which would be hard to model using the atomistic spin dynamics.However, the atomistic spin dynamics is a more general and accurate way of modelling magnetic properties in and out of equilibrium and can be applied to any type of magnetic system .Schematic illustration of the multiscale approach for modelling magnetic systems.First, the electronic properties are calculated from the first principles of quantum mechanics, leading to electron model (details in Sec. 2).From this, an effective spin model can be obtained, characterized by Heisenberg J ij and Dzyaloshinskii-Moriya D ij parameters that describe the interactions between spins ⃗ S i and ⃗ S j on different atomic sites (details in Sec. 3).Finally, micromagnetic description of magnetic texture (example of a typical simulation result is shown) is derived from the atomistic spin model of the previous step (detals in Sec. 4).The orders of magnitude of the relevant length scales and some of the material properties are listed for each model, together with the basic equations in the simplified form.Model parameters are highlighted by the green color.
From the introduction above, it is clear that there is a possibility of multiscale modelling of magnetic systems following three basic steps illustrated in Fig. 2: • Electronic structure: the electronic energy spectrum and related properties (hopping parameters t ij , orbital character of the wavefunction) are calculated from first principles for a given crystal structure and chemical composition; • Atomistic spin model: the energy of the system is parameterized using a spin model and the magnetic interactions are calculated based on the electronic properties; • Micromagnetic description: discrete atomic structure of the system is neglected and a continuous magnetic model is derived from the atomistic model under the assumption of slowly varying magnetization.
In the rest of the review, different steps of this multiscale approach will be explained in general and a few representative examples from the literature will be discussed.Most of the examples are concerned with systems having a significant Dzyaloshinskii-Moriya interaction which stabilizes non-collinear or even topological textures.Importantly, this review does not aim at providing a comprehensive picture of different magnetic phenomena that can emerge in electronic system but rather focuses on the fundamental aspects of magnetism and general theoretical methodology for predictive modelling of any given magnetic system fully from first principles.

Electronic structure
Soon after the foundations of quantum mechanics were laid in 1920's, it became clear that many properties of materials (conductivity, mechanical strength, magnetism etc.) are determined by the behavior of electrons in the crystal lattice.It was also clear that understanding this behavior is a formidable task for theory due to the long-range Coulomb interactions between all the electrons in the system, which decay rather slowly as 1/r (r -distance between electrons) and can lead to complex correlation effects.Even nowadays, the full picture of electronic correlations in solids is missing, although there has been a lot of progress in theoretical description of correlations using different approaches.One such approach is based on density functional theory (DFT), which was formulated by Pierre Hohenberg and Walter Kohn in 1964 [1] and has become a widely used and very successful tool for modelling realistic materials.In this section, general overview of the DFT methodology is given with a focus on how it can be used to understand magnetic systems.

Main idea of DFT
In principle, the quantum-mechanical model of a solid crystal is well defined in a sense that the exact form of the Hamiltonian is known: This Hamiltonian contains i) the kinetic energy of N electrons each with a mass m and charge e, ii) the Coulomb repulsion between them, iii) their electrostatic attraction to the nuclei with atomic numbers Z k , iv) the kinetic energy of the N 0 nuclei with masses M k , and v) their electrostatic repulsion energy.In many cases, the Born-Oppenheimer approximation [11] is valid, which assumes that the electron dynamics is much faster than the nuclei dynamics.One can write then a purely electronic Hamiltonian for given positions of the nuclei: However, even this simplified quantum-mechanical problem is still virtually impossible to solve, since the number of electrons in a solid can approach astronomical numbers.Already a relatively small sample with a size of (10 × 10 × 10) nm can contain around ∼10 6 electrons.Even storing the information about the full electronic wavefunction would be beyond the capabilities of modern computers, considering that the amount of required memory scales exponentially with the number of electrons.
Alternative approach, density functional theory (DFT), was suggested in the seminal paper by Hohenberg and Kohn [1], who proposed to use the electronic density n(⃗ r) instead of the full correlated wavefunction Ψ to describe electronic systems.In the paper, it is proved that the energy of the system can be expressed as a functional of the density only and the ground state electron density minimizes this energy functional.This greatly simplifies the task of solving the electronic problem, because i) density is a function of 3 arguments (components of ⃗ r), while the original wavefunction depends on the positions {⃗ r i } of all electrons, ii) the ground state can be found by variational principle.The later aspect has been further elaborated on in the paper by Kohn and Sham [2] who introduced a special ansatz to reformulate the DFT equations using an auxiliary independent-particle system.Each auxiliary particle has a wavefunction ψ i and experiences an effective potential which contains contributions from different interactions appearing in Eqn. 3 (ionic potential V , Hartree potential related to density n(⃗ r) and the so-called exchange-correlation potential V xc ): These independent particles have nothing to do with original electrons, neither relate their wavefunctions ψ i to the original wavefunction Ψ of the full electron system described by Hamiltonian (3).However, the actual electron density is postulated in the Kohn-Sham theory to correspond to the total density of auxiliary particles: This density also enters each of the N Kohn-Sham equations (5), meaning that these have to be solved self-consistently with Eqn. 6.The solution is usually found iteratively by updating the electronic density and DFT potential after each step until both quantities change negligibly based on the convergence criterion.

DFT approximations
In the density functional theory, all the correlations between electrons are encoded in the exchange-correlation potential V xc (⃗ r), the exact expression for which is unknown but has been a subject of intense research over many decades.In actual calculations for real materials one has to use approximations for this functional, of which there are several hundreds.Historically, the first one was the local-density approximation (LDA) [2,12] where the local correlation energy density ε xc (n) is postulated to coincide with that of a homogeneous electron gas of the same density.The later has been calculated using stochastic simulations of the quantum-mechanical electronic problem [13].Accordingly, the exchange-correlation energy functional E[n] and potential V xc (⃗ r) in LDA read: In this form, the local-density approximation does not contain any spin dependence and n(⃗ r) is the total electron density.Importantly, the LDA approach can be extended to magnetic systems, which was done by Barth and Hedin [14].The DFT potential now becomes spin-dependent which allows solutions of DFT equations with non-equal spinup and spin-down densities n ↑ and n ↓ , such that the total density is n(⃗ r) = n ↑ (⃗ r)+n ↓ (⃗ r) and the magnetization is m(⃗ r) = n ↑ (⃗ r) − n ↓ (⃗ r): Despite assuming a slowly varying electron density, the LSDA approach has been quite useful for calculating the electronic properties of different types of solid-state systems, both in bulk and near surfaces and interfaces where the density is far from being uniform.Also, the LSDA predictions for the magnetic properties (magnetic moments and interactions, as we discuss in Sec. 3 later on) are quite accurate in many cases, for example, for elemental transition metals Fe, Co and Ni and rare earths as well as for weakly correlated transition metal compounds (for more detailed discussion see [10]).Moreover, LSDA correctly predicts that all other elements are non-magnetic in a sense that they do not have intrinsic moments, while Pd is close to being magnetic due to the so-called Stoner instability [15].This instability is triggered when the electronic density of states at the Fermi level n(E F ) is large enough, so that it becomes favourable for electrons to attain a finite spin polarization which reduces the total energy (more discussion of fundamentals in [16,17]).The Stoner model is useful for understanding the origin of magnetism in solids in simple terms.It contains the kinetic energy of electrons, which increases when a finite magnetization is induced due to increase of Fermi energy, and magnetic exchange energy which is negative and modelled as E m = −m • I, where m is the total magnetization and I is the Stoner parameter that can be calculated from density functional theory, which was done, for example, for 3d and 4d elements in [18].The competition between the two energies determines whether a system is magnetic or not, and the Stoner term −m•I is a simple model of complex electronic correlation effects which, together with quantum statistics for the wavefunction (more details in Sec.3), can lead to the emergence of magnetism.The condition for this is n(E F ) • I > 1 which marks the onset of the aforementioned Stoner instability where the gain in magnetic energy outweighs the kinetic energy growth.
A natural way to try to improve the local-density approximation is to take into account the density variations in the first order by including the density gradient in the DFT potential, leading to the generalized-gradient approximation (GGA) [19,20].The underlying idea is to expand the exact DFT potential as a functional of electronic density in a series where the zero-order term is the LDA expression (7) and the first-order term is GGA.Next term contains the Laplacian of density (∇ 2 n) and forms the basis of the meta-GGA method [21].The process of increasing the accuracy of DFT by constructing more complex density functionals is usually illustrated by the Jacob's ladder (Fig. 3a) where each higher step represents a better approximation to the exact DFT functional.In principle, DFT approximations with higher-order terms should provide more accurate results than LDA.However, higher-order terms have to be constructed such that several properties of the exact DFT functional (so-called sum rules [22]) are satisfied.
Even though GGA is built in this way, it does not always improve the theory predictions for the structural and magnetic properties compared to LDA, which probably has to do with error cancellation.For example, while LDA in general underestimates the lattice parameters of most systems, GGA usually overestimates it by a similar amount.This kind of error can be critical when studying systems the functional properties of which are sensitive to structural details, such as ferroelectric where the predicted stability of ferroelectric polarization can change dramatically depending on the DFT approximation [23].Possible ways to improve the DFT description in such cases are i) HSE06 hybrid functional [24] and ii) PBEsol functional [25].In case of HSE06, a portion of exact exchange energy (around 25%) is added to the PBE functional, which also improves the theory prediction for the electronic band gap [26].In case of PBEsol, the predicted lattice parameter lies in-between the LDA and GGA predictions and is usually very close to the experimental value.This is useful for modelling the mechanical properties or phenomena induced by pressure, because the PBEsol predictions match better the measurements at zero pressure.It should be also noted that the computational complexity of PBEsol is comparable to usual PBE and is much lower than that of HSE06.
In practice, the Kohm-Sham equations (5) of DFT are solved in a certain basis, for example, plane waves (Quantum Espresso code [27,28]), projector augmented waves (VASP code [29,30]), linear muffin-tin orbitals (RSPT code [31,32]) etc.Moreover, one can differentiate between all-electron and pseudopotential implementations of DFT codes where the core states are either treated on the same footing as the valence electrons or are excluded from equations by working with pseudized wavefunctions which show a smoother behavior near the atomic cores.In general, all-electron codes are considered the most accurate, while the accuracy of pseudopotential codes depends on the quality of pseudopotentials, which are generated based on all-electron results.There can be some differences in DFT theory predictions depending on the particular implementation aspects discussed above, and a detailed analysis can be found in [33].

Beyond-DFT approaches
For systems where electronic correlations play a significant role, it may be required to include further corrections in the DFT functional.Lack of correlations in the original LDA and GGA functionals can, for example, lead to false predictions of metallic behavior or a considerably underestimated value of the band gap (see e.g.chapter 4 in [34]).This is especially the case for transition-metal and rare-earth compounds with localized d or f electrons.In these systems, electron interactions correspond to energy scale between 1 − 10 eV.In contrast, the exchange splitting in LSDA which produces a finite spin polarization is of the order of 1 eV (see for example Table II in [35]) which corresponds to the Hund's exchange and explains the failure of LSDA (and GGA) for moderately and strongly correlated systems.The DFT+U approach [36][37][38] can be very useful in such cases, since it explicitly introduces additional Coulomb repulsion for chosen states in the spirit of the static mean-field Hubbard model.
The total energy in the DFT+U method contains the usual DFT part (LDA or GGA) and two additional terms, the Coulomb repulsion term and double-counting correction.
a) The Couloumb repulsion term reads (according to [39]): Here, nγ 1 γ 2 are on-site occupation numbers expressed in the combined orbital and spin space and U γ 1 γ 3 γ 2 γ 4 are the matrix elements of Coulomb interaction which, in the basis of atomic-like orbitals, can be written in terms of a few parameters (see page 13 of chapter 6 in [40]).
In Eqn.(9), there are terms that add an energy penalty when two electrons of opposite spin occupy the same orbital and atomic site.This favors half-filled orbitals and leads to finite spin polarization induced by electronic correlations.In case DFT already predicts a magnetic state, the U -related term leads to enhancement of magnetic moments and possibly can open or increase the electronic band gap.For example, DFT calculation for antiferromagnetic insulating FeO oxide predicts a metallic system (Fig. 3b), while the DFT+U method correctly predicts an insulating behavior (Fig. 3c) with a band gap of 2 eV according to the band structure (Fig. 3d) and increases the Fe moment, bringing it closer to the measured value of 1.9 µ B .Similar improvement can be achieved for other antiferromagnetic oxides (NiO, CoO etc.) and many further systems (chapter 4 in [34]).The choice of the U parameter can be made based on the agreement between theory and experiment in terms of predicted structural, electronic and magnetic properties or based on independent calculations using constrained randomphase approximation [43] (calculation examples in [44][45][46][47][48]).Usually, the U parameter is around a few electron volt for moderately correlated systems and can approach 10 eV for strongly correlated systems, such as oxides (Table I in [36]) and rare-earth elements and compounds [49][50][51].
There are also other terms in Eqn. 9 which correspond to the Hund's coupling between different orbitals of the same atomic site and favor parallel alignment of spins on those orbitals.More detailed discussion of these two kinds of interactions will be done in Sec. 3 where we will see how one can derive an effective spin model from the electronic Hamiltonian.
b) The double counting term is necessary, because the approximate exchangecorrelation functional in DFT already contains some portion of true electronic correlation and one has to subtract this contribution before adding the Coulomb repulsion term ( 9) to avoid double-counting (DC) correlations.There is no unique way of introducing DC corrections and its choice depends on the nature of correlated orbitals.Two widely discussed limiting cases are the fully localized (FLL) and around mean-field (AMF) limits.The FLL double-counting correction is discussed in [39]: The around mean-field limit is considered, for example, in [52][53][54], and its simplified form reads: where l is the orbital quantum number.In both expressions (10) and (11), the quantities n ↑ and n ↓ are related to the occupation numbers in the Coulomb term (9) through the definitions n ↑ = Tr(n ↑ γ 1 γ 2 ) and n ↓ = Tr(n ↓ γ 1 γ 2 ), and n = n ↑ + n ↓ .In general, DFT+U corrections make the exchange-correlation potential orbitaldependent and, in case of FLL (Eqn.10), favor integer orbital occupations.This leads to shifts of band energies by an amount proportional to the U parameter, such that orbitals with n > 1/2 are shifted downwards in energy, while orbitals with n < 1/2 are shifted in the opposite direction (see discussion in [55]).
The DFT+U method outlined above addresses the electronic correlations on the static mean-field level in the spirit of Hartree-Fock approach, where quantum operators are replaced by their average values.A significant improvement over this approximation is offered by the dynamical mean-field theory [56][57][58] which also takes into account the temporal fluctuations of orbital occupations on different atomic sites.In this theory, it is possible to describe a metal-insulator transition induced by electronic correlations even in the paramagnetic phase.This is in contrast to the static DFT+U approach where long-range magnetic order might be necessary in the calculation in order to obtain an insulating behavior.Combination of DFT and DMFT is a powerful method of studying realistic material models and has been extended over the last two decades to include also the non-local correlations effects.Since DFT+DMFT is not in the focus of this review, the interested reader is referred to the reviews [56][57][58] and book [40], for example.

Origin of magnetism
While we talked a lot about magnetism in electronic correlated systems in the previous section, we did not discuss much its physical origin.This section will focus on this topic in the context of magnetic interactions in general, which can appear in such systems and can lead to long-range magnetic order.
From the point of view of classical physics, it is difficult to understand why some systems can be become magnetic in the presence of external field, not to mention permanent magnets like Fe, Co and Ni which have finite total magnetization even in the absence of external field.The problem with classical description was pointed out by Niels Bohr in 1911 (see publication in collected works [59]) and later on by Hendrika Johanna van Leeuwen in 1921 [60], and their results are nowadays referred to as Bohr-Van Leeuwen theorem, which states that the magnetization vanishes in any classical system placed in external magnetic field.The starting point is the classical Hamiltonian for N interacting electrons (see, for example, discussion in [17]): which contains the vector potential ⃗ A i ≡ ⃗ A(⃗ r i ) and interactions between all the electrons.When calculating the classical partition function Z N , which involves integration over the phase space of the ⃗ p i and ⃗ r i variables, one can substitute variables ⃗ A i and show that the partition function, in fact, does not depend on the vector potential.Since the average magnetization M = k B T ∂ ∂H ln Z N is proportional to the derivative of Z N over magnetic field, the magnetic moment vanishes exactly (M = 0).
In order to explain magnetism of electrons in solids, one has to take into account their quantum-mechanical behavior and replace (12) with: where ⃗ p i and ⃗ r i are now the momenta and position operators which do not commute.As discussed in [17], this quantum Hamiltonian can be derived as an approximation to Dirac's equations.From the later, one also concludes the existence of electron spin which leads to additional terms proportional to Pauli matrices ⃗ σ in Eqn.(13).The Hamiltonian also contains a coupling between the electron spin and orbital moment ⃗ l i .
So far, our discussion addresses an electron system in applied magnetic field.However, even without it some systems are permanent magnets due to the presence of magnetic interactions.It is interesting to see where such interactions can come from, considering that the underlying Hamiltonian ( 13) is not spin-dependent for the simple case of H = 0 and ζ = 0.The well-known answer is that quantum statistics of electronic wavefunction due to the Pauli exclusion principle together with electrostatic Coulomb repulsion between electrons produces the so-called exchange interaction.The idea is that two electrons sitting on the same orbital of the same atomic site cannot have the same spin, otherwise they have to be spatially separated by some average distance R 1 and occupy different orbitals.This configuration has, in fact, lower Coulomb energy e 2 /R 1 compared to the case of two electrons with opposite spin, which can be in the same orbital and, accordingly, at a smaller average distance from each other R 2 < R 1 (Fig. 4a).Because of this, the energy of the system becomes spin-dependent and can be effectively described by the exchange interaction which favors a ferromagnetic alignment of spins on the same site, known as the Hund's coupling.

Direct exchange
Let us extend this consideration onto the case of two electrons occupying two orbitals, each on a different site (Fig. 4b).This situation is quite different, because now electrons can hope between the two sites leading to six possible electronic configurations: |↑↓, ⟩, | , ↑↓⟩, |↑, ↓⟩, |↓, ↑⟩, |↑, ↑⟩, and |↓, ↓⟩, meaning that either both sites are half-occupied or one of them is empty, while the two electron spins can be parallel or antiparallel (detailed discussion of this model is found, for example, in [34], chapter 7 by Erik Koch).It is more illustrative to write the Hamiltonian for this problem in the second quantization, instead of the first-quantization adopted in Eqn. ( 13): The essential components of this model are i) electron hopping with energy t and ii) Coulomb repulsion energy U between two electrons on the same orbital/site.Diagonalization of this Hamiltonian in the basis of six two-electron states mentioned above leads e.g. to a triplet eigenstate represented by |↑, ↑⟩, |↓, ↓⟩ and superposition (in the leading order of perturbation theory).Thus, in contrast to Hund's coupling, direct exchange between two sites favors an antiferromagnetic spin alignment.
If we want to build a model of this system with only states that have a fixed spin S = 1 2 on each site, then we disregard the intermediate states with double occupancy of one of the sites and downfold the original Hamiltonian onto the subspace of states |↑, ↓⟩, |↓, ↑⟩, |↑, ↑⟩, and |↓, ↓⟩, as shown in detail in [34].After rewriting the resulting second-quantization expression in terms of spin operators, this procedure leads to an effective spin Hamiltonian: This downfolding from the full electronic space to a spin subspace filters out the high-energy states with double occupancy but still takes into account their effect on other states leading to effective spin-spin interactions, so the effective Hamiltonian ( 15) is spin-dependent while the original one ( 14) is spin-independent.The resulting spin model (Eqn.15) has the form of the Heisenberg model, which was derived for the first time by Werner Heisenberg in 1928 [61] for quantum S = 1/2 systems, followed by numerous discussions in the literature, for example, in the context of single-band Hubbard model for extended systems [62,63].

Beyond-Heisenberg models
The procedure of deriving effective spin model from the electronic Hamiltonian described above is general and can be applied to a variety of electronic systems.Didactic discussion of different aspects of this procedure from the point of view of perturbation theory is given in [64] (accompanying software is provided in [65]), where strongly correlated electrons on a square lattice, for example, are considered and the electron hopping between the neighboring sites (t ≪ U ) is treated as a perturbation.For S = 1/2 systems with Hamiltonian like Eqn. (15) generalized to periodic crystal lattices, it can be shown that the effective spin model contains two kind of terms: where the summations run for i ̸ = j ̸ = k ̸ = l.The first term here corresponds again to the Heisenberg model, while the additional term represents a four-spin interaction.So already in this simple case one obtains beyond-Heisenberg terms in the spin model which describe more complex excitations of the electronic system and include multi-site electron hopping processes.Different orders of the perturbation theory for S = 1/2 systems bring corrections to the exchange parameters J ij and K ijkl in Eqn. ( 16) but do not produce any further types of spin-spin interactions.
Important to note is that the spin model ( 16) derived from the electronic Hamiltonian is a quantum spin model with S = 1/2, so ⃗ S 1 and ⃗ S 2 are operators with commutation relationships common for angular momentum operators.
For higher-spin systems S ≥ 1, additional contributions to the spin model are possible, such as the biquadratic and three-spin interactions for S = 1 case: Due to the quantum nature of spin Hamiltonian, biquadratic exchange does not appear for S = 1/2 systems, because it can be reduced to a usual bilinear Heisenberg interaction due to the properties of spin-ladder operators (see for example [refs]).Similarly, on-site anisotropy of the form K U ( ⃗ S i • ⃗ e) 2 vanishes exactly for S = 1/2 systems, since this terms is equivalent to applying twice the ladder operators Ŝ+ and Ŝ− to |m z = −1/2⟩ and |m z = +1/2⟩ states, respectively.In general, terms up to order ( ⃗ S i • ⃗ S j ) 2S are allowed in the quantum spin Hamiltonian for spin number S.
Adding spin-orbit coupling to the original electronic Hamiltonian in the form λ( ⃗ L i • ⃗ S i ) and including it in the perturbation theory expansion would lead to new types of magnetic interactions in the effective spin model, such as the Dzyaloshinskii-Moriya (DM) interaction ) and its higher-order variants, e.g.
).The bilinear DM interaction H DM can exist in magnetic systems where the inversion symmetry is broken for bonds connecting magnetic moments.Such symmetry breaking can be especially strong near surfaces or interfaces, and well-known examples are transition metal multilayers, e.g.Pd/Fe/Ir(111) [66], Pd/Co/Pd [67][68][69][70], Pt/Co/Ta [71], Ir/Fe/Co/Pt [72] and exchange-biased multilayers [73], to name just a few.The DM interaction favors non-collinear magnetic orders with a certain chirality defined by the DM vectors ⃗ D ij and, historically, was proposed in [74,75] to explain the phenomenon of canted magnetism in α-Fe 2 O 3 , MnCO 3 and CoCO 3 .For certain crystal symmetries, it can even stabilize topogical magnetic textures with a size below a hundred nanometers (for reviews see [76][77][78][79][80]).More details about the requirements on symmetry are given in Sec. 4 of this review.

LKAG first-principles approach
While the approach described above provides, in principle, the values of magnetic interactions, these are rather approximate, because they are obtained within the perturbation theory where t/U is assumed to be a small parameter (t -electron hopping, U -strength of electronic correlations).More accurate approach was suggested in the LKAG paper [81] where the idea is to calculate the energy change ∆E ij due to a small perturbation of the reference magnetic state.More specifically, this perturbation includes an infinitesimal canting of two spins on sites i and j (Fig. 2b).In the nonrelativistic case, the energy change ∆E ij is then proportional to the corresponding Heisenberg exchange parameter J ij .The later can be calculated from perturbation theory, but this time the advantage is that the canting angle is, indeed, a small parameter, in contrast to t/U , so the LKAG approach is more accurate and is applicable, in principle, to any magnetic system with intrinsic moments.Using the language of Green functions, the Heisenberg exchange in the LKAG approach can be expressed as follows: As we see, the Heisenberg exchange depends on the electron Green function for two sites (i and j) which is reminiscent of electron hopping parameter t discussed on previous pages, and is proportional to the spin splitting ∆ of electronic states on these sites.Both quantities ( Ĝij and ∆) contain the information on the electronic structure, making the LKAG approach material-specific, and are reprensented by matrices in orbital space.Correspondingly, the trace in Eqn. ( 19) is equivalent to summing over different orbitals.However, it is possible to analyze orbital-resolved contributions by skipping the trace operator and internal matrix multiplications, leading to expression [82,83]: which is written now in terms of summation over Matsubara frequencies ω n .
The LKAG approach is a widely used and efficient method for calculating magnetic interactions in electronic systems, both periodic and non-periodic.It allows to study the long-range character of interactions while using the minimal chemical unit cell, in contrast to the total energy mapping method which always requires constructing large supercells leading to high computational costs.Detailed discussion of magnetic interactions, ways to calculate them and examples for solid state systems can be found in the comprehensive review [84].In the original LKAG paper [81], the authors also discussed how to estimate the Curie temperature and spin stiffness, that can be directly compared to experiment, and calculations for Fe, for example, showed a good agreement between theory and experiment.Furthermore, the exchange parameters obtained in the LKAG approach correspond to energy variations at zero temperature and, for that reason, naturally describe the magnon excitations with a high accuracy (for a review, see [85]).
Regarding the long-range character of Heisenberg interaction, the LKAG approach has been useful for analyzing the mechanism of magnetic exchange in transition metals Cr, Mn, Fe, Co, and Ni.In [82,83], orbital decomposition of Heisenberg exchange revealed oscillations of t 2g -related contributions as ∼ sin(k F r)/r 3 according to the Ruderman-Kittel-Kasuya-Yosida mechanism [86][87][88] and short-range character of e grelated contributions (Fig. 6).The t 2g contributions are affected by the properties of the Fermi surface; for example, the Fermi wavevector k F determines the oscillation period of the Heisenberg exchange over distance, while the e g contributions depend on the electron hopping parameter resembling the double exchange.Such analysis makes the connection between the electronic properties and magnetic interactions clearer and helps to understand and predict how the magnetic interactions change in response to structural and chemical variations.
The LKAG formula (19) has been extended for the relativistic case to provide information on the Dzyaloshinskii-Moriya interaction and for correlated systems [89][90][91][92], such that, for example, the J xy component of the generalized exchange tensor looks like this: The important difference here is the presence of frequency-dependent self-energy Σ i (iω n ), which can be calculated within the dynamical mean-field theory approach [56][57][58].Furthermore, the Pauli matrices σx and σy describe the non-collinear character of the spin configuration due to the canting of the considered spins along the xand y-directions.From the off-diagonal components like J xy ij one can calculate the Dzyaloshinskii-Moriya interaction ⃗ D ij and symmetric anisotropic exchange Γij , while the Heisenberg exchange is obtained from the diagonal components J αα ij , so that the full exchange tensor is written as follows (see discussion in [93]): and the effective spin model is represented by the generalized Heisenberg model: where ⃗ S i are unit vectors showing the direction of the local spin axis on site i.This means that the J ij parameters calculated in the LKAG formalism already contain the magnitude of the magnetic moments on individual sites.This should be kept in mind when comparing these J ij values to the results obtained using other methods.
It should be noted that the magnetic interactions in the LKAG approach are evaluated for a certain spin configuration called reference state.Very often it makes sense to use the lowest energy spin configuration as the reference state or a magnetic state which is close to the ground state, if the later is too complicated to treat within DFT (e.g.large-wavelength spin spiral).If the reference state is non-collinear and different from the ground state, then it is necessary to apply a constraining field to stabilize such a configuration.However, in that case, for DFT-based calculations of exchange interactions one has to bear in mind that the constraining field is different from the gradient of the total magnetic energy obtained from DFT [94].This is reflected in the way that the electronic system is mapped onto an effective spin model, as elaborated upon in recent works [95,96].

Basic equations
Although effective spin models, discussed in the previous section, are a powerful tool for studying the time-dependent behavior and temperature-dependent properties of magnetic systems based on atomistic spin dynamics equations (ASD, for details see Sec. 5), there are limits to the size of spin systems that can be modelled on modern computing architectures.For large systems, the simulation time scales almost linearly with the total number of spins N spins and, in practice, such atomistic simulations are feasible for up to a hundred million spins.For 2D systems described by (n×n) supercells (N spins ∼ n 2 ), this gives an excellent opportunity to study complex spin textures, such as domain walls, spin spirals, skyrmions etc. on the length scale of up to a few hundred nanometers.However, for 3D systems the total number of spins grows faster with the system dimensions (N spins ∼ n 3 ), which reduces significantly the size of systems that can be modelled atomistically, in terms of simulation time and required memory.
In order to address the magnetic properties of larger systems, one can switch to the micromagnetic description where the underlying atomic structure is completely disregarded in the simulation and the magnetization is treated as a continuous vector field ⃗ m(⃗ r) defined at each point of space, not just at atomic sites.Of course, in the actual simulation ⃗ m(⃗ r) is discretized on e.g. a 3D grid (n × n × n), and the simplification here, compared to ASD simulations, is related to the fact that each of the n 3 micromagnetic elementary cells usually covers around a few nanometers of the system in each spatial direction.Thus, it requires less micromagnetic cells than actual atomic spins to model a system of a given size.
Further below, the basic equations and definitions of the micromagnetic approach are derived in a didactic way, followed by details about numerics and concrete examples from literature.

Heisenberg interactions.
Let us first derive the micromagnetic model corresponding to the Heisenberg exchange (first term in Eqn. 1).Following the methodology discussed in Ref. [97], one can start from the atomistic picture of spins on sites i and j and write ⃗ S j = ⃗ S i + ( ⃗ S j − ⃗ S i ).The difference in the brackets by definition is equal to ⃗ m(⃗ r j ) − ⃗ m(⃗ r i ), since the magnetization in different regions is described by the micromagnetic function ⃗ m(⃗ r) which can be expanded in a Taylor series: Here, ⃗ R ij is the vector in real space connecting spins on sites i and j, and ⃗ m ≡ ⃗ m(⃗ r i ).The magnetic energy for the Heisenberg exchange can be written now as follows: After writing out the expression in the brackets, we see that the first term is a constant energy contribution, since it is proportional to ⃗ m 2 = 1 (constant length of magnetic moment vectors is assumed).The second term with a 1 st -order derivative can be rewritten: and is equal to zero, since m α m α = 1.
For that reason, the first non-vanishing term in Eqn.(25) reads: Usually, only the diagonal terms (α = β) are considered in micromagnetics and, for equal diagonal elements A xx = A yy = A zz = A (e.g. in cubic systems), the magnetic energy density is written as follows: The energy is proportional to the so-called spin stiffness A: which characterizes the overall strength of magnetism in the system.For example, for ferromagnets (A > 0) one expects that the Curie temperature increases with the spin stiffness A. We also see that the whole complexity of magnetic interactions J ij between different spin neighbors is boiled down to just one number A, which reflects the simplifying character of the micromagnetic approach.However, the spin stiffness A is obtained from atomistic interactions J ij which are calculated for a given crystal structure and chemical composition from the first principles of quantum mechanics.This means that the micromagnetic model is also obtained from first principles and reflects the properties of a concrete material, which is a great advantage of the multiscale approach.
Dzyaloshinskii-Moriya interactions.Now let us derive the micromagnetic energy density ε DM originating from the DM interactions, following the same methodology as outlined above for the Heisenberg exchange.The atomistic spin model (second term in Eqn. ( 1)) is the starting point: Again, we replace ⃗ S i with the micromagnetic function ⃗ m ≡ ⃗ m(⃗ r) and expand ⃗ S j up to the 1 st -order term: where ⃗ R ij is the distance vector between spins i and j.Substituting this into Eqn.(30) and noticing that ⃗ m × ⃗ m = 0 leads to a micromagnetic energy density: where one can define the spiralization matrix D ≡ D αβ (α, β = x, y, z) as The spiralization matrix can have nine non-zero components and reflects the crystal and magnetic symmetries.For systems with chiral crystal structure, such as cubic B20 compounds MnSi [98], FeGe [99][100][101][102] and Fe 1−x Co x Si [103], the spiralization matrix is diagonal (D xx = D yy = D zz = D) and this scenario is referred to as the so-called bulk DM interaction.By substituting the matrix elements D αβ = D δ αβ one can obtain the micromagnetic energy: This type of DMI favors the formation of Bloch skyrmions (Fig. 7a), which are localized magnetic objects with topologically non-trivial winding of atomic spins (further discussion in Sec.5.2).
Different scenario takes place for transition metal multilayers, such as Pd/Fe/Ir(111) [66] and Pd/Co/Pd [67][68][69][70], which have the C 3ν crystal symmetry.In that case, the symmetry of the DM interaction vectors leads to non-zero off-diagonal elements D xy = −D yx = D, while other elements are zero.This scenario corresponds to the so-called interfacial DM interaction with a micromagnetic energy which stabilizes Néel skyrmions (Fig. 7b): This result agrees with the Lifshitz invariants expected for the C 3ν crystal symmetry, as discussed, for example, in [104] (Eqn.8), [105] (Eqn.6) and [106] (Table I).The same micromagnetic expression is obtained also for polar crystals, such as lacunar spinels GaV 4 S 8 [107] and GaV 4 Se 8 [108], which at low temperature both have the C 3ν symmetry.
There is another class of materials, Heusler compounds, which have D 2d crystal symmetry.Magnetic compounds of this type, Mn 1.4 Pt 0.9 Pd 0.1 Sn [111,112] and Mn 2 NiGa [113], show the so-called anisotropic DM interaction characterized by a DM matrix with the only non-vanishing components D xy = D yx = D. We notice here that both off-diagonal elements have the same sign, as opposed to the interfacial DMI discussed above.This kind of structure of DMI leads to the micromagnetic energy (notice sign changes compared to Eqn. 35): which favors antiskyrmions (Fig. 7c).In terms of the spin winding, the later are inbetween the Néel and Bloch skyrmions, because the spins are winding perpendicular to the radius vector from the center of the antiskyrmion in one direction and parallel to it in other directions.Interestingly, antiskyrmions have been also observed in Fe 1.9 Ni 0.9 Pd 0.2 P which has S 4 symmetry [114].
While concrete simulations using atomistic spin dynamics or micromagnetics (see Sec. 5) are necessary to predict the stability of different topological magnetic objects  [109], figure g) -from [110], and figure h) -from [107] (Fig. 7a-c) for given material parameters and external conditions (temperature, applied field etc.), it is often useful to look at the Dzyaloshinskii-Moriya micromagnetic matrix D αβ (which is computationally easier to obtain than to perform dynamic simulations) to see what kind of topology can be supported by the DM interaction.
It is important to note one can combine the atomistic spin dynamics (ASD) description and micromagnetic approaches to address phenomena on different length scales in the same system within the same simulation [97,[115][116][117].For example, one can describe more accurately the interaction of magnetic textures (domain walls and skyrmions) with atomistic defects (impurities, dislocations etc.) by applying ASD in the region around the defect and micromagnetics (MM) in the rest of the system.The two regions are connected by an interface region which can be constructed in different ways to provide as seamless transition between the two regions as possible [97,115,118].In recent work [117], this approach was demonstrated, for example, for a skyrmion moving near a triangular defect with locally larger anisotropy under the influence of spin-transfer torque (STT).In case of weak STT, the skyrmion is pinned at the defect in the ASD-MM simulations, while pure MM calculation results in the annihilation of the skyrmion.At larger STT, the skyrmion goes past the defect in both methods, but the MM simulations show some reduction of skyrmion size which is not seen in the more accurate ASD-MM approach.Another example in that work [117] is concerned with a linear dislocation which, similarly to the previous example, can pin a skyrmion, if STT is not strong enough.The important difference is that the skyrmion can go around the dislocation several times, depending on the magnitude of STT and damping constant, meaning that at some point it moves against the direction of STT.This unusual behavior predicted by the ASD-MM simulations cannot be reproduced by purely MM approach.In general, the ASD-MM approach combines the advantages of both methods: i) the accuracy of the ASD and ii) the larger system size that can be simulated micromagnetically.

Numerical calculation of micromagnetic parameters
Despite the seeming simplicity of formulas ( 29) and ( 32) defining the micromagnetic parameters A and D, their numerical evaluation for real systems based on the magnetic interactions calculated from first principles can be a challenge [93,110,[119][120][121].In particular, magnetic interactions in metallic systems have a long-range character.This might not be apparent at first glance, as illustrated for ferromagnetic Fe bcc in Fig. 8a, where it seems that the Heisenberg interaction J ij is insignificant for distances above ∼ 10 Å.However, it should be kept in mind that the average number of neighbors grows with distance R ij and, moreover, factors of R 2 ij and R ij enter the expressions for the spin stiffness and spiralization matrix (Eqns.29 and 32).As we discussed in Sec.3.4, certain contributions to the Heisenberg exchange in ferromagnetic transition metals come from the t 2g − t 2g orbital interactions which have an RKKY character and oscillate with distance as ∼ sin(k F r)/r 3 (Fig. 6).This leads to a slow convergence of the sums ( 29) and ( 32) over different neighbors with respect to the real-space cutoff distance.As shown in Fig. 8c, the spin stiffness A oscillates significantly (green curve) with varying cutoff R even at larger distances of 5 lattice constants (a), which prevents an accurate determination of the spin stiffness parameter.In fact, the estimate of A based just on the nearest-neighbor Heisenberg exchange is close to experiment (A exp = 320 meV • Å2 ) only by lucky chance (first green point in Fig. 8c).
More than 20 years ago, a solution to this problem was suggested in [122] where the idea is to introduce an exponential factor in the definitions of the micromagnetic parameters: Here, parameter µ is chosen to be positive, so that the exponential factors lead to faster decaying terms in the sums (37) and improve the numerical convergence with respect to the real-space cutoff.Converged results are obtained starting from a certain value µ c , which is usually around 1.0−3.0(depending on the system), and based on a set of A(µ) and D(µ) calculated at different µ ≥ µ c one can extrapolate to the limit µ → 0. The extrapolation can be done using a 3 rd -order polynomial or exponential functions of µ, as discussed in the literature, for example, for lacunar spinel GaV 4 S 8 [107] and B20 compounds [93].
For the simpler, yet illustrative, example considered in this review, Fe bcc, this procedure gives 298 meV • Å2 for the spin stiffness (see the dashed line crossing the y-axis at µ = 0 in Fig. 8d) which is quite close to the measured values of (280 − 330) meV • Å2 [123,124].In contrast, the direct summation of Eqn. 29 for spin neighbors within a distance of up to 5 lattice constants gives the estimate around 140 meV • Å2 which is more than a factor of two smaller.Another example is the B20 compound FeGe for which very different theory estimates of the spin stiffness and micromagnetic DM parameter were obtained in the literature (see discussion in [121]).
These micromagnetic parameters show again large oscillations as functions of the realspace cutoff distance (see Fig. 7 in [121]) and the theory estimates are quite far from the measured values.The reason for this discrepancy is so far unclear, but one may speculate that higher-order magnetic interactions may play a role in this system as well as the presence of multiple sublattices [121].
The method outlined above is one way of determining the micromagnetic parameters from first principles and relies on the real-space representation, as can be seen from formulas (29) and (32).There are other methods as well, for example, the one discussed in [81,125], which is suitable for periodic crystalline systems and based on the k-space representation.The idea is to calculate the energy δE = D αβ q α q β of a spin spiral with a small q-vector using the multiple-scattering theory which leads to the expression for the spin stiffness: where T is the scattering path operator in the multiple-scattering theory [126].Accurate calculation of D αβ depends now on the number of k-points used for the Brillouin zone summation in this expression, instead of the real-space cutoff in Eqns.29 and 32.Another advantage of this formulation is that one can describe disordered alloys within the coherent-potential approximation.Overall, this method provides estimates of the spin stiffness which are similar to the real-space method discussed above and which lie within the range of measured values.Similar approach for calculating the spin stiffness in the k-representation has been developed recently using the transport theory [127].

Basic equations
Deriving the effective spin or micromagnetic models of magnetic systems is just one of the steps towards the actual modelling which is supposed to address the magnetic ground state at given conditions as well as the magnetization dynamics when these conditions are varied in time.This is the motivation of this final chapter of the review and we will consider here the fundamental aspects of magnetization dynamics and, more specifically, examples of different magnetic states that are observed in systems with Dzyaloshinskii-Moriya interaction.Theory successes in modelling such systems from first principles and the perspectives on predicting new systems with non-collinear or even topologically non-trivial magnetism will be discussed as well.
Similarly to the way that an effective spin model can be derived from the electronic problem, the equations describing the dynamics of such spin models can be derived from the time-dependent quantum-mechanical Kohn-Sham equations.The key ingredient of this derivation (discussed in detail in book [10]) is a term proportional to σ • ⃗ B eff which describes the spin-dependent part of the Hamiltonian and contains the effective magnetic field due to intrinsic magnetic interactions and external field.Another important aspect is the assumption that the magnetization density is well localized and locally collinear around each atom, such that its dynamics can be characterized by a single vector ⃗ m i indicating the local magnetization direction.This atomic moment approximation is quite accurate for most systems, with a few exceptions being fcc Pu [128] and Cr monolayers [129], and allows to represent approximately the complexity of the spin density distribution in a given system by a discrete set of magnetic moments ⃗ m i on different atomic sites.Comprehensive and didactic discussion of the derivation of atomistic spin dynamics equations and many other aspects of magnetization dynamics are covered in paper [9] and book [10].
Both in the atomistic and in micromagnetic models, the magnetization dynamics can be described by the Landau-Lifshitz-Gilbert (LLG) equation [130,131]: which describes the time-evolution of magnetic moments ⃗ m i (t) on different atoms, molecules or other effective magnetic entities (numbered with index "i") or in different micromagnetic regions.The system-specific parameters here are the Gilbert damping constant α, which characterizes the energy dissipation, and the effective field ⃗ B i .The later contains contributions from the external magnetic field as well as the Heisenberg and Dzyaloshinskii-Moriya interactions, dipole-dipole energy, on-site anisotropy etc.At finite temperature T , random field ⃗ B i,f l proportional to √ αT is added to the effective field ⃗ B i to include the fluctuation effects in the simulation.The connection between the damping and fluctuations in the spin system characterized by the quantities α and ⃗ B i,f l is not accidental and is a consequence of the fundamental fluctuation-dissipation theorem [132,133].The resulting stochastic differential equations are solved using the method of Langevin dynamics [134] by rewriting the LLG equation with the fluctuating field as a Fokker-Planck equation (solution methods for these equations can be found, for example, in [135]).
It is the accuracy of calculating the effective field ⃗ B eff , which mostly determines the predictive power of the atomistic spin dynamics simulation for a given system.The ⃗ B eff field is determined from the spin Hamiltonian as − ∂H ∂ ⃗ m i [130].Correspondingly, one would get different expressions for the atomistic and micromagnetic cases, where the starting points are the Heisenberg model (1) and micromagnetic functional (2).Different contributions to the effective field are summarized in Table 1 for these two approaches.Higher-order spin-spin interactions highlighted in Sec.3.3 would make further contributions to the effective field and can be evaluated using the general gradient formula mentioned above.It should be noted that, in the approximation of constant-moment length (| ⃗ m i | = 1), the component of the effective field parallel to the magnetic moment ⃗ m i does not exert any torque and vanishes from the LLG equation (39).Another important aspect of first-principles spin dynamics is concerned with the concept of the constraining field, which, as we already mentioned in Sec.3.4, is used to Table 1.Summary of the effective field related to the Heisenberg, Dzyaloshinskii-Moriya (DM), dipole-dipole interactions and uniaxial anisotropy for the atomistic and micromagnetic pictures.

Origin
Atomistic picture Micromagnetic picture stabilize a non-collinear spin configuration different from the magnetic ground state.In a recent work [94], the main conclusion is that the constraining field is not always equal to the gradient of the Hamiltonian, if the later contains mean-field-like contributions or is based on density functional theory.This can also have important consequences when the magnetic exchange parameters are determined from the constraining field approach [94][95][96].

Prediction of magnetic properties and textures
Based on the effective spin model, one can calculate the corresponding thermodynamic quantities using either the Monte Carlo (MC) stochastic sampling or the atomistic spin dynamics (ASD) described by the LLG equation.In case of MC approach, a series of quasi-random perturbations of the spin configuration is generated using the Metropolis algorithm (more details in [136], pages 17-25) and the thermodynamic properties are calculated as averages over this data set, which is supposed to represent the variety of system configurations and their probabilities to a good accuracy.The later depends on the number of Monte Carlo samples.In Fig. 9b, the MC prediction for the temperature-dependent magnetization M (T ) and susceptibility χ(T ) is shown for bulk FeGe compound (structure in Fig. 9a), which is one of the skyrmionic compounds [99-102, 137, 138].The simulation is based on the spin model obtained in [110], which includes the Heisenberg and DM exchange interactions between spin neighbors up to a distance of 3 lattice parameters (∼ 14.1 Å).The magnetization M (T ) decreases monotonically with temperature and approaches zero around T C ∼ 260 K, while the magnetic susceptibility reaches the maximal value at that temperature, which is a typical behavior of ferromagnets.Similar result is obtained in ASD simulations with the same exchange interaction parameters (Fig. 9c).The theoretical estimate for the critical temperature is close to the experimental value T C ∼ 280 K [99][100][101]137].
It is important to note that in the classical spin simulations the decrease of magnetization at low temperatures is linear, which is in disagreement with the quantum Figure 9. a) Crystal structure of the skyrmionic B20 compound FeGe and the DM vectors for the nearest-neighbor bonds (reproduced from [110]).Temperaturedependent magnetization and susceptibility are obtained from the b) Monte Carlo and c) atomistic spin dynamics (ASD) simulations using the spin model from [110].d) Magnetization and e) heat capacity of ferromagnetic hcp Gd as functions of temperature obtained from Monte Carlo simulations based on different types of statistics (classical, quantum and mixed); plot a) is reproduced from [110] and plots d) and e) are reproduced from [139]; experimental data points are from Ref. [140].
picture where magnon excitations reduce M (T ) by a quantity ∼ T 3/2 in case of 3D systems.This artifact of classical simulations also leads to a constant finite heat capacity at low temperatures, which reflects the classical equipartition theorem.However, in the quantum picture the heat capacity should approach zero as T 3/2 .Ways of correcting this inconsistency were suggested, for example, in [141][142][143] where one of the ideas is to rescale the simulation temperature based on the average magnon energy.
The ASD simulation with this pure quantum statistic predicts well only the lowtemperature part of the M (T ) curve [141][142][143], while the classical statistics is reasonably accurate at higher temperatures, around and above the critical temperature T C .In a recent work [139], the authors proposed a combination of both approaches using a mixed statistics which interpolates smoothly between the quantum and classical limits.That is, in the Monte Carlo simulations at low temperature the probability W qt (∆E, T ) of changing the state of the system, accompanied by energy change ∆E, is defined using quantum statistics [141][142][143], while at high temperatures near and above the magnetic transition the probability W cl (∆E, T ) is described by the Boltzmann distribution.Inbetween the two limits T = 0 and T ∼ T C , the probability is a linear interpolation with a factor α = 1  and g) are reproduced from [110] and figures f) and h) are reproduced from [144].
result, the mixed-statistics Monte Carlo simulation predicts a temperature dependence of the magnetization M (T ), on the example of ferromagnetic hcp Gd, which is in a much better agreement with measurements (Fig. 9d) and a correct behavior of the heat capacity (Fig. 9e).
Let us discuss how the atomistic spin dynamics (ASD) and micromagnetic approach can be used to search for the magnetic ground state of a given system.In both formalisms, a possible strategy is to perform the annealing simulation which starts from a random magnetic configuration at high temperature and the magnetic dynamics is simulated using the LLG equation (39) while the temperature is gradually lowered down to zero (illustration in Fig. 10a).Such simulations usually result in a local energy minimum which can bare many features of the actual magnetic ground state.We demonstrate this procedure on the concrete example of Mn monolayer on W(001) surface.Literature studies [144][145][146] indicate a particularly strong DM interaction in this layered system induced by the Mn/W(001) interface (similarly to Mn/W( 110) system [145,147]) with large spin-orbit coupling due to W and inversion symmetry breaking.For the nearest-neighboring spins, the DMI is almost 50% of the Heisenberg exchange.Also, the uniaxial on-site anisotropy is very large ∼ 2.5 meV per Mn site, according to ab initio estimates from [110].Using the first-principles values of the magnetic interactions and other parameters from that work [110], we compare further below the resulting ground state found in ASD.For the smallest considered simulation cell of dimensions (25 × 25) the ASD simulation predicts a single-domain spin-spiral phase (Fig. 10b) with a small spatial period, similarly to experimental result [144].When a larger (100 × 100) simulation cell is used, the spin structure splits into several domains oriented either in [110] or [1][2][3][4][5][6][7][8][9][10] directions.Both domain types have the same energy and, for that reason, can coexist in the system.The distribution of domains depends on the initial spin configuration, which is chosen randomly, and the cooling rate.Importantly, for lower cooling rate, i.e. longer simulation time at each temperature, the number of domains decreases (example with just two domains in shown in Fig. 10c) and the required number of time steps in the simulation depends on the cell size.For that reason, single-domain state is achieved faster in smaller (25 × 25) cell, while the larger (100 × 100) cell splits into a multi-domain state for the same cooling rate.For comparison, similar spin configuration is predicted by Monte Carlo (MC) simulations at 13 K in [144] shown in Fig. 10h in comparison with the spin-polarized scanning transmission microscopy image (Fig. 10f) at the same temperature.The magnetic state predicted by the MC simulations consists again of spin-spiral domains oriented in [110] and [1][2][3][4][5][6][7][8][9][10] directions, in accordance with measurements.
Apart from spin spirals, the DM interaction can also stabilize topologically nontrivial magnetic states, for example, skyrmionic phases with either isolated skyrmions or skyrmion lattices.Being originally proposed as topological defects in high-energy physics by Tony Skyrme in 1961 [148], skyrmions have been identified also in other research fields, e.g.liquid crystals [149] and superconductors [150].Skyrmions are compact magnetic objects formed by atomic spins which are winding around a certain direction, such that the spins in the center of the skyrmion are opposite to those at its rim and the intermediate spins rotate smoothly in-between the center and the rim (Fig. 7ac).This winding is associated with a non-trivial topology which can be characterized by the topological charge: where vector ⃗ m = ⃗ m(⃗ r) describes the direction of magnetization at a given point ⃗ r in space.For Bloch and Néel skyrmions (Fig. 7a,b) the topological charge equals +1, while it is −1 for antiskyrmions (Fig. 7c).Several material classes with different types of topological magnetic textures are known (see e.g.reviews [76][77][78][79][80]) and searching for new systems of that kind is an on-going effort in the research community.
Theoretical simulations based on first-principles approaches can provide useful insights into the physical mechanism of stability of topological magnetism and even make quantitative predictions of materials properties.In the following, we consider a few literature examples illustrating this statement.
The first example is the nanoskyrmion lattice observed in Fe monolayer on Ir(111) surface.Spin-polarized scanning transmission tunneling microscopy experiments showed  c) are reproduced from [151] and figure d) -from [66] and e) -from [152].
that the Fe moments in this system form a square lattice of skyrmions with a period of ∼ 1 nm [151].Interestingly, this skyrmion lattice is incommensurate with the underlying hexagonal lattice of Fe sites.From the analysis of experimental data, one can conclude that the observed magnetic state is a superposition of two spin spirals with Q-vectors ⃗ Q 1 and ⃗ Q 2 which are at 92.2 o to each other and have the same magnitude 0.277 × 2π/a.Theoretical analysis in [151] based on an effective spin model ( 16) with an addition of DM interaction revealed that the nearest-neighbor four-spin interaction favors the double-Q magnetic state combining the ⃗ Q 1 and ⃗ Q 2 spirals.This magnetic superposition can lead to either skyrmions or antiskyrmions, which in the non-relativistic case are degenerate in energy, but due to the specific symmetry of the DM interaction in this system skyrmions are lower in energy.Using the magnetic interaction parameters calculated from first principles, the authors quantified the effect of the Heisenberg, fourspin and DM interactions on the formation of the nanoskyrmion lattice and obtained a good agreement with experiment in terms of the magnitude of the ⃗ Q 1 and ⃗ Q 2 vectors that minimize the magnetic energy.
It is interesting to see how one can drastically modify the ground state of this layered nanoskyrmion system by simply placing a single Pd layer on top of Fe.Experimentally [66], it is known that this changes the ground state to a spin spiral with a period around 6−7 nm which can be transformed to a skyrmion lattice of similar length scale in applied magnetic field around 1 T (Fig. 11d).Theoretical calculations in [152] revealed that the Fe/Ir (111) interface still dominates the DM interaction, which determines the chirality of spin spirals, but the presence of the Pd layer enhances the ferromagnetic Heisenberg exchange between the nearest-neighbor Fe spins by more than a factor of 2, compared to the Fe/Ir(111) system.This is based on spin-spiral energies calculated for both systems (Fig. 11e).Furthermore, there are competing antiferromagnetic interactions for more distant neighbors which favor a spin spiral state even without the spin-orbit coupling.At zero field, both factors lead to stability of spin-spiral state instead of skyrmions.Another statement from that work is that the agreement with experiment in terms of spiral wavelength, skyrmion size and critical magnetic fields for these two phases is especially good for the hcp stacking of the Pd layer on top of Fe layer.On the other hand, the fcc stacking would lead to deeper energy minimum of a spin spiral (Fig. 11e) leading to higher magnetic field required to stabilize a skyrmion lattice.It is worth noting that the strong spin-orbit coupling (SOC) due to the Ir (111) surface is essential for the skyrmion stability, as demonstrated by calculations for Fe/Pd(111) multilayer where SOC is an order of magnitude weaker and is not enough to stabilize even spin spirals (Fig. 11e).
These literature studies of Fe/Ir (111) and Pd/Fe/Ir (111) systems demonstrate that the multiscale approach can be useful for understanding the physical mechanism of the observed magnetic phenomena and how they react to variations of the material properties.With some adjustment of theory approximations and models, the agreement between theory and experiment can be achieved even on the quantitative level.
Let us now consider another example where theory predicts an interesting magnetic phenomenon which is later on confirmed experimentally.Theory analysis of 2D isotropic magnets with DMI in [154,155] showed that not only skyrmions but also antiskyrmions can be stabilized, but in an extremely narrow range of applied magnetic field.In a more recent work [153], micromagnetic simulations revealed that in the 3D generalization of this model antiskyrmions become more stable due to the dipolar field contribution leading to easy-plane anisotropy and the possibility of non-collinearity in the direction perpendicular to the skyrmion plane, which enhances the DMI energy.However, the skyrmion-antiskyrmion pairs are less stable than the isolated topological magnetic particles.Furthermore, skyrmionium has been found to be a stable solution in the 3D simulations (Fig. 12a).These theory predictions were confirmed successfully by experiments in the same work [153] done on thin FeGe samples with a thickness of ∼ 70 nm.Lorentz transmission electron microscopy (TEM) images (Fig. 12b) show clearly how the original zero-field helical spiral state transforms into skyrmions and antiskyrmions (Fig. 12c) in applied magnetic field (skyrmionium is also stabilized at a slightly stronger field) and how the pairs of these topological particles annihilate when in close proximity to each other (see futher figures in [153]).The simulated TEM images (Fig. 12a) are in a good agreement with the measurements, which helps to identify the nature of topological magnetic textures observed in experiment.It turns out that size effects related to the thickness of the FeGe sample are important for the formation of antiskyrmions, which can be stabilized if the thickness does not exceed the spiral wavelength.Thicker samples, on the other hand, only show usual skyrmions [153].

Spin-lattice dynamics
It is worth noting that the LLG equations (39) are derived in the adiabatic limit where the spin degrees of freedom are assumed to be much faster than the ionic degrees of freedom, which is true in many cases.Material-specific spin dynamics in this approximation with parameters obtained from first principles has been considered in [9,[156][157][158][159][160][161], while some of the first works [162][163][164] on spin dynamics included also the coupling between the spin and ionic sub-systems.
More recently, methodology for simulating coupled spin-lattice phenomena is suggested in [165] where the idea is to combine the LLG equations (39) and molecular dynamics simulation in the harmonic-phonon approximation.The full Hamiltonian of the system, discussed in that work [165], consists then of the pure spin (H S ) and lattice (H L ) parts and the coupling term H LS .The spin part is based on the original Heisenberg model where the exchange tensor J αβ ij is expanded in a Taylor series as a function of small displacements ⃗ u of different atoms: This expression contains the pure spin part H S and a part of the spin-lattice coupling in the form of cross-terms depending both on the magnetic moments and ionic displacements.The lattice contribution to energy contains the kinetic energy of the nuclei and takes into account the dependence of the ionic forces on the magnetic configuration by expanding the force constants into a Taylor series with respect to small deviations of magnetic moments ⃗ m i from the equilibrium state: Combining both expansions ( 41) and ( 43) allows to write the spin-lattice coupling in the following general way: while the pure spin and lattice parts are defined in this way: Also, the procedure for determining all the coupling coefficients in this expression is proposed in [165].The idea is to consider a supercell where one of the atoms is displaced and the magnetic interactions are calculated as functions of the displacement, making it possible to calculate the first and second derivatives in Eqn.(44).Alternative ways to calculate these coefficients, for example, using a model for the exchange parameters were proposed in [166][167][168][169]. Evaluation method for the spin-lattice coupling coefficients using the Green function method, similar to the LKAG approach [81], and embedded cluster method have been proposed recently in [170] and [171], respectively, and tested for several representative systems.For example, analysis of spin-lattice coefficients in [170] indicated the importance of Dzyaloshinskii-Moriya interactions in Fe bcc, induced by lattice vibrations, for the transfer of angular momentum between the spin and lattice subsystems.
Modelling the coupled spin-lattice dynamics at finite temperature requires solving the equation of motion both for the spins and for the lattice sites simultaneously, leading to an extended system of equations compared to Eqn. (39): The fields acting on the spin and lattice degrees of freedom are calculated in a similar way based on Hamiltonian H total = H S + H L + H LS defined by Eqns.44-46: Additional contribution to these fields comes from fluctuations at finite temperature which can be approximately described by white (Gaussian) noises with correlation functions (µ, ν = x, y, z): B f l i,µ (t)B f l j,ν (t ′ ) = 2D M δ ij δ µν δ(t − t ′ ) (50) The important parameters here are D M and D L which characterize the strength of fluctuations but, on the other hand, through the fluctuation-dissipation theorem [132,133] are intimately related to the energy dissipation in the spin and lattice subsystems: where α and ν are the spin and lattice damping parameters.
It should be noted that the possibility of continuum formulation of the equations of coupled spin-lattice dynamics has been discussed in recent work [172].The idea is similar to the micromagnetic approach (Sec.4 of this review) which introduces a continuous vector field ⃗ m(⃗ r) describing the magnetization distribution.For the spinlattice case, one can also define the mechanical energy in terms of the strain and stress tensors ε αβ (⃗ r) and σ αβ (⃗ r), using the well-established elasticity theory [173], as well as cross-terms describing the coupling between the strain tensor and magnetization and their gradients.These cross-terms can be obtained systematically starting from the atomistic model (Eqns.[44][45][46] and expanding the continuous variables in real space. Getting back to the more recent work [165], it is interesting to see how combined spin-lattice atomistic simulations, using all the aforementioned equations ( 41)-( 52) work in practice.The general methodology developed in that work is demonstrated, for example, for elemental ferromagnet Fe bcc.Specific form of the spin-lattice coupling corresponding to magnetostriction mechanism is considered in that work and all the Hamiltonian parameters are obtained fully from first principles using simulation cells with several thousands of atoms.From these calculations, an important conclusion was that one has to take into account both two-site (i, j = k) and three-site (i, j ̸ = k) terms and that the magnitude of the coupling coefficients ⃗ A ijk decays with distance slower than the Heisenberg exchange parameters J ij .The coupled spin and lattice dynamics with these first-principles parameters showed, for example, that the magnetization of Fe is somewhat reduced due to spin-lattice effects, especially near the magnetic transition temperature (Fig. 13a).Magnon frequencies, which are calculated from the spin-spin dynamic structure factor, also show interesting effects due to spin-lattice coupling, such as the softening by a few percent for temperatures in the range 300 − 700 K.It should be noted that the effect of thermal fluctuations that reduce the average magnetization lead to much larger magnon softening already in the pure spin dynamics simulations (see Fig. 17a in [165]).
The magnitude of spin-lattice effects in Fe bcc is small but it can be more pronounced in systems close to structural or magnetostructural phase transition where the magnetic interactions and possibly even the macroscopic magnetic order are sensitive to structural details.As mentioned in [165], possible candidate systems could be Laves compound YCo 2 , Invar Fe-Ni alloy, and fcc phase of Fe, where the later is known to host a variety of magnetic phases.Recently, another class of materials has been suggested as candidates for strong spin-lattice coupling, namely, the quasi-twodimensional VOCl [174][175][176] and similar compounds with Sc, Ti, Cr and Fe instead of V, and Br instead of Cl [177][178][179][180][181]. Theoretical calculations in [182] showed that this interesting physics is preserved also in the single-layered version of these systems.
Another candidate for large spin-lattice coupling is CrI 3 , which is one of the first discovered two-dimensional magnets and has topological magnon states [184].The later are due to the magnon gap ∼ 4 meV, also observed in theoretical simulations [92] but underestimated in magnitude.Spin-lattice coupling in this system has been discussed in terms of phonon spectra being sensitive to the magnetic structure [185] as well as magnetic interactions changing in response to structural distortions [183].It appears that the frequencies of some of the Raman-active phonon modes in CrI 3 are different for the ferro-and antiferromagnetic configurations (Fig. 5ab in [185]).This effect can be attributed to the difference of the electronic bands and character of the band gap for the two examples of possible magnetic configurations, as indicated by density functional theory calculations accompanying the experimental findings in [185].In another work [183], it was shown that Heisenberg interaction can even change the sign from FM to AFM for sufficiently large Cr displacements ∼ 0.2 Å (Fig. 14b,e).Similar effect under lattice strain is predicted in [186].Although average thermal  displacements in the real system are not expected to reach this value (see Fig. 8 in [183]), the magnitude of magnetic interactions will still change considerably during the coupled spin-lattice dynamics already at room temperature.The high complexity of CrI 3 magnet is also due to the sensitivity of Heisenberg and Dzyaloshinskii-Moriya interactions to displacements of non-magnetic ligands, not just the Cr sites.These spin-lattice effects can be understood in terms of the competition between AFM coupling of t 2g orbitals and FM coupling of t 2g and e 2g orbitals [183] (Fig. 14d).

Conclusions
We have discussed in this review how magnetic systems can be modelled on different length scales using the multiscale approach, starting from the electronic properties on the scale of individual atoms and electron hopping between them, proceeding with effective atomistic spin models representing magnetic excitations and textures in the system for a scale of up to several hundred nanometers and finishing with micromagnetic description of magnets in the continuous approximation where system sizes of several micrometers can be modelled.The important point is that all these different levels of multiscale modelling are connected to each other, because larger-scale models are derived from smaller-scale models, meaning that the model parameters are obtained from first principles and inherit the electronic structure feature of a given material.This improves the predictive power of the multiscale approach, as we have discussed on several examples from literature where the studied systems show interesting magnetic textures (spin spirals, skyrmions etc.).Furthermore, recent advances in modelling coupled spin-lattice phenomena were highlighted as well.The current successes and further development of the multiscale approach will ensure that many more magnetic systems and unusual phenomena will be discovered theoretically, leading to new applications and fundamental knowledge.

Figure 2
Figure2.Schematic illustration of the multiscale approach for modelling magnetic systems.First, the electronic properties are calculated from the first principles of quantum mechanics, leading to electron model (details in Sec. 2).From this, an effective spin model can be obtained, characterized by Heisenberg J ij and Dzyaloshinskii-Moriya D ij parameters that describe the interactions between spins ⃗ S i and ⃗ S j on different atomic sites (details in Sec. 3).Finally, micromagnetic description of magnetic texture (example of a typical simulation result is shown) is derived from the atomistic spin model of the previous step (detals in Sec. 4).The orders of magnitude of the relevant length scales and some of the material properties are listed for each model, together with the basic equations in the simplified form.Model parameters are highlighted by the green color.

Figure 3 .
Figure 3. a) Jacob's ladder of DFT exchange-correlation funtionals with progressing accuracy.Reproduced from thesis [41] and adapted from the talk "Basics of DFT" of K. Burke and L. Wagner during the ELK-2011 conference.Comparison of the DFT and DFT+U predictions for b,c) the density of states and the d,e) band structure of the antiferromagnetic Mott insulator FeO (reproduced from [42]).

Figure 4 .
Figure 4. a) Illustration of the Pauli exclusion principle and exchange interaction for two electrons on the same atomic site.b) Comparison of magnetic energies of antiparallel (singlet state |s⟩) and parallel (triplet state |t⟩) spin configurations of two electrons on different sites (each with one orbital).For the antiparallel state the energy is lowered due to existence of electron hopping (with energy t) between the sites and intermediate state |↑↓, ⟩ with both electrons on one site.The energy penalty for having both electrons on the same site and orbital equals U .

Figure 5 .
Figure 5. a) Electronic and b) effective spin models for a periodic crystal lattice (reproduced from [64]).a) The electronic Hamiltonian corresponds to the Hubbard model and contains electron hopping terms, Coulomb interactions between electrons on the same site in the same orbital (U ) or different orbitals (U ′ ) as well as the Hund's coupling J H . b) The spin model derived from the Hubbard model includes, for example, the bilinear two-site (J) and four-site (K) interactions, three-spin exchange (Y ) and the biquadratic interaction B.

Figure 6 .
Figure 6.a) Orbital decomposition of Heisenberg exchange in magnetic transition metals; b) different contributions from e g and t 2g orbitals.c) RKKY character of Heisenberg interaction between the t 2g orbitals in Fe bcc and short-range character of interactions involving the e g orbitals.Reproduced from [82].

Figure 7 .
Figure 7. Schematic illustration of a) Bloch and b) Néel skyrmions and c) antiskyrmion that are stabilized by the d) bulk, e) interfacial and f) anisotropic Dzyaloshinskii-Moriya (DM) interactions (the corresponding DM matrices are shown as well).Examples of systems with such interactions are g) B20 compound FeGe, h) multiferroic lacunar spinel GaV 4 S 8 and i) Heusler compound Mn 2 NiGa.The relevant crystal structures and symmetries of each compound class are shown as well.Figures a) and b) are reproduced from [109], figure g) -from[110], and figure h) -from[107]

Figure 8 .
Figure 8. a) Heisenberg interaction J ij in Fe bcc for spin neighbors at different distances R ij .b) Heisenberg interaction multiplied by R 3 ij vs the distance R ij ; the oscillating character of the plotted quantity indicates a significant contribution of the RKKY interaction.c) Spin stiffness A of Fe bcc vs the cutoff distance for the summation in Eqn.37 with different µ values.d) Spin stiffness A of Fe bcc vs the exponential decay parameter µ (see Eqn. 37).Separate data point on the y-axis indicates the result of direct summation with µ = 0, while the dashed curve shows the extrapolation to µ → 0 limit based on the converged sums at finite µ-values.

Figure 10 .
Figure10.a) Schematic illustration of the annealing procedure in atomistic spin dynamics (ASD) simulations.Relaxation of the spin structure is done at gradually lowered temperature until T = 0 K is reached.b) Resulting spin structure (singledomain spin spiral) for Mn/W(001) layered 2D system with uniaxial anisotropy simulated using a (25 × 25) cell.d) Spin configuration obtained in a larger simulation cell reveals coexisting spin-spiral domains oriented in[110] and[1][2][3][4][5][6][7][8][9][10] directions.c) Allowing for longer simulation times and more temperature steps in annealing leads to larger spin-spiral domains.e) Crystal structure of Mn/W(001) layered system with the DM vectors shown for the nearest-and next-nearest-neighbor bonds.f) Experimentally measured magnetic state[144] showing slightly distorted spin-spiral domains.g) Top view of the Mn monolayer structure and DM vectors.h) Magnetic structure predicted by Monte Carlo simulations at 13 K. Figures e) and g) are reproduced from[110] and figures f) and h) are reproduced from[144].

Figure 11 .
Figure 11.a) Crystal structure of Fe monolayer on Ir(111) surface.b) Schematic illustration of the nanoskyrmion lattice found in this system.c) Magnetic energy calculated for the double-Q state as a function of the Q-vector magnitude; contributions of the Heisenberg exchange, four-spin and DM interactions to the total energy are shown.d) Schematic representation of the magnetic phases observed in Pd/Fe/Ir(111) multilayer at increasing applied magnetic field.e) Calculated spin-spiral energy for various multilayers as functions of the Q-vector.Figures b-c) are reproduced from[151] and figure d) -from[66] and e) -from[152].

Figure 12 .
Figure 12. a) Isosurfaces (m z = 0) of the spin textures in FeGe system with dimensions (512 × 512 × 50) nm 3 from micromagnetic simulations showing skyrmion, antiskyrmion, skyrmion-antiskyrmion pair and skyrmionium.The figure below is the calculated over-focused Lorentz TEM image for these spin textures.b) Measured Lorentz TEM images of FeGe at 95 K in zero field and H = 200 mT (sample size: 1 µm × 1 µm).c) Combined experimental and theoretical phase diagram showing the skyrmionic, antiskyrmionic and helical spiral phases of FeGe; the symbols are the experimental data points; the red-shaded region represents the skyrmion lattice phase and the vertical dashed line at 287 K shows the Curie temperature.All figures are reproduced from [153].

Figure 13 .
Figure13.a) The reduced magnetization of Fe bcc at different temperatures modelled with pure spin (solid line) and spin-lattice (dashed line) dynamics; the inset show the scaling of the Binder constant with the system size, which is used for more accurate estimate of the ordering temperature (T C = 850 K). b) The relative softening of the magnon frequencies due to spin-lattice interaction at different temperature; magnons are considered along a high-symmetry path in the Brillouin zone.Reproduced from[165].

Figure 14 .
Figure 14.a) The crystal structure of CrI 3 monolayer.b,e) The variation of Heisenberg interactions as functions of Cr displacement along different directions.d) Different orbital contributions to the Heisenberg interaction and the direct exchange between two Cr moments (grey orbitals).Figures are reproduced from[183].
Figure 14.a) The crystal structure of CrI 3 monolayer.b,e) The variation of Heisenberg interactions as functions of Cr displacement along different directions.d) Different orbital contributions to the Heisenberg interaction and the direct exchange between two Cr moments (grey orbitals).Figures are reproduced from[183].