Atomistic-continuum multiscale modelling of magnetisation dynamics at non-zero temperature

In this article, a few problems related to multiscale modelling of magnetic materials at finite temperatures and possible ways of solving these problems are discussed. The discussion is mainly centred around two established multiscale concepts: the partitioned domain and the upscaling-based methodologies. The major challenge for both multiscale methods is to capture the correct value of magnetisation length accurately, which is affected by a random temperature-dependent force. Moreover, general limitations of these multiscale techniques in application to spin systems are discussed.


Introduction
In recent years, atomistic spin dynamics (ASD) became an essential tool for understanding experimentally observed behaviour of magnetic materials [17]. In this approach, the dynamics of spin magnetic moments of individual atoms is modelled. In contrast to ASD, the classical approach of understanding ferromagnets, which originates from the first half of twentieth century, is micromagnetics [10]. It operates with the volume-averaged quantities, such as magnetisation. The area of application of these two descriptions is somewhat different: ASD aims at describing the phenomena at the nanometre scale, while the micromagnetic theory is usually applicable at scales larger than the micrometre. Due to relatively high computational cost of ASD, it is useful to construct multiscale methods that combine ASD and micromagnetics in a single unified framework. An extended background and discussion of the purposes of such multiscale models can be found in [4,27].
At finite temperatures, the rotational motion of the spin magnetic moments, which have the constant length, is described by a system of Langevin equations. The presence of temperature is modelled by a Gaussian noise term. In this case, each individual computational realisation is affected by a random force. However, the physically meaningful quantities of interest are the statistical averages of the magnetic moments, which have reduced lengths due to the cancellations in the averaging. For an accurate description of the dynamics, the reduction in the magnitude of the magnetisation or spin magnetic moments must be carefully accounted for. A traditional way of computing these statistical averages is the Monte Carlo method, see e.g. [25], which is computationally expensive, as large number of replicas of the solution is required to compute statistically meaningful quantities. An alternative way is to derive explicit equations for the evolution of statistical averages of the spin magnetic moments. One such approach is to use the Landau-Lifshitz-Bloch (LLB) equation, which describes the dynamics of the statistical averages, see e.g. [19]. A similar technique (a variant of the LLB) is proposed in [8] where the length (given by the so-called power law) is embedded into a macroscopic equation similar to LLB. These approaches typically contain some restrictive closure arguments, , where E denotes an ensemble average and m is the spin magnetic moment, which are not necessarily true in general [32]. To avoid such restrictive arguments, one may either resort to fully atomistic simulations, or design multiscale numerical algorithms which do not use these assumptions.
One of the most common types of concurrent multiscale methods is the partitioned domain approach [31]. Within this method, the computational domain is split into several regions with different material descriptions, e.g. atomistic and continuum descriptions [26]. The biggest challenge for these approaches is the construction of an error-free coupling at the interface separating the domains. There are several initial attempts of using partitioned domain atomistic-continuum multiscale models in application to micromagnetics [20,22,3,27].
Another common type of multiscale modelling is based on upscaling, where the microscopic scale effects are systematically upscaled to an initially incomplete macroscopic model describing the macroscopic scale behaviour. In this case, the microscopic description has to be consistent with the current macroscopic variables, which leads to a two-way coupling between the microscale and macroscale models. There are several multiscale frameworks, such as the heterogeneous multicale method (HMM) [1,16] or the equation-free approaches [23], which use this idea. The upscaling-based strategy 1 that is used in this paper rests upon the HMM methodology, where the microscale model is typically known, but expensive to solve over the entire macroscopic domain of interest, while the macroscale model is incomplete, as it lacks certain macroscopic quantities, which might be time-dependent. These missing quantities are then computed by carrying out the microscopic simulations in small spatial and temporal domains (sometimes referred to as representative volume elements). The methodology is efficient in the case of problems, which possess an inherent scale separation. The efficiency of the method is due to the fact that the size of the microscopic domains may be chosen comparable to the smallest length scale present in the system.
The aim of the present study is to describe the magnetisation dynamics using efficient atomistic-continuum multiscale formalisms, in which number of degrees of freedom is significantly reduced in comparison to a full atomistic simulation. Moreover, the overall ambition is to design algorithms, which do not suffer from restrictive closure arguments, and, at the same time, capture the correct dynamic behaviour. The microscale model, considered in the present work, is a discrete system of atomistic particles interacting at a finite temperature. For the purpose of this study, it is assumed that an atomistic model provides the exact material behaviour, although this may not necessarily be true in reality. In a multiscale method, a microscale model has to be consistently coupled to a macroscale model, which is addressed in this work under two different multiscale frameworks; the domain-partitioning and the HMM. Another aim of this work is to demonstrate the shortcomings of the standard continuum approach for a number of specific cases, which motivates the necessity of the multiscale modelling. 1 Here, the term "upscaling-based strategies" is used to refer to the class of methods such as HMM or the equation-free methods, where a two-way coupling between the micro and macro scales is necessary for an accurate description of the macroscopic phenomena.
In Section 2 of this paper, a number of issues in relation to atomistic-continuum transition at non-zero temperatures are discussed, and the domain partitioning approach is presented. In Section 3, the detailed mathematical formulation of the HMM approach for the case of micromagnetics is given. Afterwards, in Section 4, various multiscale examples are considered.

Towards multiscale modelling
The consistency between atomistic and continuum descriptions is required in any type of multiscale coupling. However, the construction of a deterministic continuum model that corresponds to the stochastic atomistic description is a complicated problem. Several issues related to estimation of the macroscopic quantities of stochastic atomistic systems are highlighted in subsequent subsections. The problem of selecting an appropriate continuum model is further discussed in subsection 2.4.

2.1.
Modelling approaches. In this subsection, atomistic and continuum modelling approaches are summarised. For simplicity, a 1D spatial arrangement is considered. Full descriptions can be found in [18] and [2], respectively. In what follows, the vector quantities are always denoted by bold face letters, e.g. m ∈ R 3 , and the Euclidian 2-norm is written as |m| = m 2 x + m 2 y + m 2 z .
2.1.1. Atomistic spin dynamics. In the atomistic approach, the dynamic behaviour of spin magnetic moments of individual atoms is described by the atomistic Landau-Lifshitz-Gilbert equation [9,18]: where γ is the gyromagnetic ratio, λ is the phenomenological (Gilbert) damping constant, m i is the direction of spin magnetic moment, µ is the length of spin magnetic moment, J ij are constants of Heisenberg exchange interaction between atoms i and j, K a is the anisotropy constant, p a is the anisotropy axis (uniaxial anisotropy case is considered) and H e is the external field. Thermal excitations are taken into account by adding a stationary stochastic field with the following statistical properties: where ρ and ν are the Cartesian coordinates of h i , k B is the Boltzmann constant and T is temperature. More details regarding the stochastic term can be found in [11]. Parameters λ, µ, J ij , K a and p a can be computed from electronic structure calculations [17], and are considered to be constant for a certain material. At the continuum scale, the dynamics of the magnetisation is modelled by the continuum version of the Landau-Lifshitz-Gilbert (LLG) equation [2,14]: where M is the normalised magnetisation field and β L and α L are equal to atomistic parameters in (2.1). At zero temperature, exchange parameter can be obtained directly from the atomistic parameters: where r ij is the distance between atoms i and j, and the sum is over all atoms with which atom i interacts (also A e is assumed to be spatially constant in this equation). Since the anisotropy term is local, the same anisotropy parameters K a and p a are used in the continuum and the atomistic equations. The accuracy of these choices is discussed in the following subsections. At finite temperatures, the continuum model has to be modified. These modifications differ depending on the approach and are further discussed in Section 2.4.
The continuum magnetisation field M is supposed to be equal to some spatial average of atomistic m i . However, due to the nature of the LLG equation, unitlength vectors are used in the formulations above. In general, atomistic-continuum transition introduces an error to the solution, which is dependent on magnetisation gradient, for details see [27].

2.2.
Direct numerical simulations of macroscopic quanteties for a 1Dsystem, using the atomistic model. In this subsection, averaged quantities of 1D atomistic spin systems are investigated. The results highlight some problems with the standard continuum approach, which cannot capture the behaviour of some 1D atomistic systems.
The parameters of the dynamic deterministic continuum model, which is presented above, can be constructed out of the parameters of the dynamic deterministic atomistic model, which results in a certain solution error that depends on the magnetisation gradient [27]. In this case, the first step towards building a consistent continuum model that describes material behaviour at finite temperatures is the replacement of the stochastic atomistic description with a deterministic atomistic description with effective 2 parameters, from which the continuum parameters can be determined using known relations, e.g. (2.8). In this subsection, it is shown 2 Effective parameters are such parameters that when substituted into an alternative model, result in an alternative solution, which is close to (or exactly the same as) the original solution provided by the original mode. In the present context, this means that the ensemble average of the solutions given by the stochastic atomistic model is equal to the solution of the deterministic atomistic model, which includes effective parameters.
parameter µ J K a p a λ material 1 7.63µ B 1.28 · 10 −21 11.86 · 10 −24 e y 0.1 material 2 7.63µ B 1.28 · 10 −21 0 -0.1 Table 2.1. Material parameters that were used in the simulations. Here, µ B = 9.274 · 10 −24 is the Bohr magneton. that this approach results in an inaccurate model: such effective parameters should depend not only on temperature, but also on the state of magnetisation and other properties of the system, such as boundary conditions. At finite temperatures, the mean-field approach [2] is usually used, in which it is assumed that the statistical averages of atomistic spins interact with the same exchange constant (J ij ) as atomistic spins at 0 K, while the length of spins is decreased due to averaging [2], i.e. |m i | becomes less than 1. From a numerical point of view, this can be implemented by replacing m i with m * i s i , where |m * i | = 1. If this assumption had been always valid, the only transition from the stochastic atomistic model at finite T to the deterministic atomistic model would have been estimation of s i . However, s i depends on the properties of the atomistic system, which is shown below. In subsection 2.2.2, it is demonstrated that the mean-field approach can give quantitatively different results in comparison to a direct numerical simulation.
It was already shown elsewhere that the effective exchange parameter depends on temperature [7]. In subsection 2.2.3, it is shown that it also depends on the state of magnetisation, therefore it is impossible to estimate an effective exchange parameter that is valid in a general case even for a particular temperature and a particular system.

2.2.1.
Dependence of macroscopic quantities on system size. When material is modelled at an atomistic level, macroscopic quantities, which are extracted from the atomistic simulations, can depend on the size of the atomistic system. For example, in [24], the dependence of the magnetisation length, which was extracted from the atomistic simulations, on the size of the system was analysed. Such dependence exists for all types of boundary conditions, including periodic, but it is commonly accepted that there is a convergence with respect to the system size. In [24], in addition to material anisotropy, a weak external field was used to stabilise the atomistic system and align it in a certain direction. However, the presence of the external field has a far more profound effect on the macroscopic quantities, and there is a significant non-linear dependence of the magnetisation on the value of the external field. In Figure 2.1, it can be seen that in the case of 1D arrangement of atoms (with periodic boundary conditions) not only the magnetisation length 3 depends on the magnitude of the external field, but also the time required to reach the equilibrium. Moreover, when the external field is absent, convergence (within given time interval and for given range of system sizes) was not observed, although the material uniaxial anisotropy was taken into account. When an external field is present, convergence is observed.  Parameters that were used in the simulation are summarised in Table 2.1, material 1, and Table 2.2, set 1, and are taken from [18]. Since this paper is methodologyoriented, units are ommitted for all quanitites, which is common in numerical analysis. Only nearest-neighbour interaction was taken into account. The implicit mid-point method [15] was used for time-stepping with a time step ∆tD = 0.005. The initial directions of m i were selected to be e y . The relatively low magnetisation value for this temperature can be explained by the 1D nature of the problem: each atom interacts only with two neighbours.

2.2.2.
Dependence of macroscopic quantities on temperature. As discussed above, there is a strong dependence of the macroscopic quantities on temperature. There are several ways to quantify this dependence, e.g. the mean-field approach [2]   is made, which is not valid in general [32]. Moreover, due to its limitations, the analytical approach cannot provide either time required to reach the thermodynamic equilibrium or dependence of macroscopic quantities on the system size.
In Figure 2.1, in addition to numerical results, analytical estimates are shown. Mean-field estimate of m was obtained by solving the following equation [2]: where L is the Langevin function and N n is the number of nearest neighbours, which in this case is 2, see the appendix for a derivation of this formula in the presence of an external field and exchange interaction. The mean-field magnetisation length is significantly larger than the length of the magnetisation that is directly obtained from modelling of the dynamics of spin systems, as seen in Figure 2.1. Moreover, the mean-field m is less affected by the magnitude of the external field.

2.2.3.
Non-linear dependence of the effective exchange coefficient on the magnetisation gradient. In most cases, the macroscopic quantities of a stochastic spin system appear to be non-linearly dependent on external conditions. The non-linear dependence with respect to the external magnetic field was demonstrated in previous subsections. In this section, the importance of boundary conditions is highlighted. Specific boundary conditions are used here to create a magnetisation gradient, which influences effective exchange coefficient in a non-linear way. The approach is somewhat similar to the "domain wall stiffness approach" [7,21].
In the case of mechanical behaviour of 1D atomistic chain, a relevant macroscopic characteristic is the stiffness of the chain. It can be obtained numerically by applying forces to the boundary atoms and calculating the displacements of atoms. This method can be extrapolated to spin systems, where the interatomic interaction is Heisenberg exchange. An exchange stiffness of the entire 1D spin system is introduced. The spin system is twisted by the external magnetic field applied only to the boundary atoms.
To illustrate the approach, the 1D chain of N atoms interacting via Heisenberg exchange only with the nearest neighbours is considered. The system does not have magnetic anisotropy, K a = 0; and open boundary conditions are used (atoms 1 and N interact only with one neighbour). External field H e = He y is applied to atom 1, external field H e = H (e x sin φ + e y cos φ) is applied to atom N , while external field H e = 0 is applied to all other atoms. The problem is schematically illustrated in Figure 2.2.
At 0 K, the equilibrium configuration for this problem can easily be determined: where θ is the solution of the following equation: This means that there is a non-linear dependence between the magnitude of the applied field, H, the applied field divergence angle, φ, and the spin moments divergence angle at the equilibrium, θ. At finite temperature, the equilibrium configuration can only be found by numerical simulation. After performing a number of simulations, it was observed that the orientation of the spin magnetic moments, χ i , of the ensemble average of the solution depends linearly on i, as well as in the deterministic case; however, the divergence angle, θ was different. Therefore, the behaviour of such material can be approximated by the deterministic model with a different (effective) parameter. When this deterministic model with an effective exchange constant J * , is subjected to the same conditions as the stochastic model with the real exchange constant J, the resulting divergence angle of spin magnetic moments at the equilibrium, θ, is the same. Hence the effective parameter J * can be determined by performing numerical simulations at finite temperature, extracting the angle θ from the ensemble average of the solution and using equation (2.11), in which J is replaced by J * .
In Figure 2.3a, the effective exchange coefficient between two spins, J * , normalised by J is plotted. It is clearly seen that the effective exchange constant depends non-linearly on the angle between the magnetic fields applied to the boundary atoms. This implies a non-linear dependence of the continuum exchange constant on the local magnetisation gradient in the case of finite temperatures.
Problem parameters are summarised in Table 2.1, material 2, and Table 2.2, set 2. A time step ∆tD = 0.005 was used. The initial directions of m i were selected to be according to equation (2.10) with real J, hence J * J −1 = 1 at t = 0. A chain of 8 atoms was considered. Angle θ was determined by a linear fit of angles χ i , which are obtained by projecting the ensemble average of the solution on xy-plane. The following value of H was selected: which minimises the statistical error of calculation of J * for small divergence angles. The angle φ was varied: values of 135 • , 90 • and 45 • were used. The coefficient J * J −1 shows the change of the exchange parameter when the constrained (i.e. with the external field applied to the boundary atoms) stochastic system is replaced by the constrained deterministic. In the case of unconstrained system (i.e. no external field), if the mean-field assumption is to be believed, the stochastic system can be replaced by the deterministic system with the same exchange parameter, but different magnetisation length. After re-normalisation m i = m * i s, |m * i | = 1, this replacement becomes mathematically equivalent to changing the exchange parameter, i.e. s = J * J −1 . For an unconstrained system with the same number of atoms, see Figure 2.3b, this coefficient is significantly larger than for the constrained system. This confirms that if greater accuracy is required, the effective parameters, which are used in the deterministic LLG equation, should not be precalculated constants (for a given temperature) but should be functions of the local magnetisation and the local magnetisation gradient.
The simulation in Figure 2.3b was performed with the same material parameters as used for the constrained system; however, without an external field; N s = 100 number of realisations was used. The coefficient s = J * J −1 in this simulation was calculated as the length of particle-averaged spin magnetic moment. The mean-field value was obtained using equation (2.9).

2.3.
Partitioned domain approach. The application of partitioned domain approach to the modelling of magnetisation dynamics is covered in [27,28]. The key idea behind the methodology is the separation of the computational domain into two regions: atomistic and continuum, which are connected at the interface. In the case of nearest-neighbour interatomic interaction, neither transition region, nor padding atoms (atoms behaviour of which are obtained by interpolation of the solution at neighbouring continuum nodes) are necessary. In the 1D case, the interface consists of a single point and the coupling is straightforward: the interface point interacts with the nearest continuum node as a continuum node and it interacts with the nearest atom as an atom (in the case of nearest-neighbour interatomic interaction). However, when the dynamic material behaviour is modelled, the introduction of an additional damping region at the interface is necessary [13]. This additional numerical damping is applied to the atoms in the proximity of the interface and acts as a low-pass filter that absorbs high-frequency waves that are otherwise reflected from the interface (due to the difference in discretisation parameters) and contribute to numerical error. The subsequent subsection follows the analogous section from [27].
2.3.1. Damping band. The equation describing the dynamics of numerically-damped atoms is the following [28]: where γ Di is the strength of the damping for atom i and m Ai is the normalised average of the solution over a certain region. The damping strength increases gradually with atomic positions approaching the interface and is proportional to the time derivative of m i : where g D is a constant determining the damping strength, r R i is the distance to the closest atom outside of the damping band, r P i is the distance to the interface atom/node and a is the interatomic spacing.
Within the damping band, the solution is attenuated to an average value, m Ai , which is calculated within a certain window surrounding atom i. The width of the window is denoted as s A . An additional average quantity m * Ai is introduced: Coefficients c are fitted to all m j within the averaging window using least squares. Finally, m Ai is obtained by normalisation:

2.4.
Selection of the macroscopic scale model for the multiscale framework. Due to the highlighted issues, the selection of an appropriate model for the macroscopic scale in the case of domain partitioning is ambiguous. First of all, based on the type of the problem, the nature of the macroscopic model has to be chosen, which can be quasistatic 4 , dynamic deterministic or dynamic stochastic. Quasistatic models are suitable for the cases when an equilibrium magnetisation distribution around an isolated atomistic region, which can be fully dynamic, has to be determined. The example of such coupling for the case of mechanical behaviour (coupling molecular dynamics and elastic finite-element continuum descriptions) is the CADD method [29]. Since it might be unphysical to model all the waves generated within the atomistic region, which eventually diffuse and are converted to heat over large distances, there is an additional advantage of using a quasistatic continuum model as it avoids such problems.
In cases when dynamic phenomena should be modelled within the continuum region, e.g. domain wall movement, one of the dynamic models has to be used. In this case, while waves are transmitted within the continuum, atomistic regions are idle and thereby reduce computational efficiency; however, this issue cannot be avoided for such type of problems. The most simple model, which can be used for the dynamic continuum, is the standard Landau-Lifshitz-Gilbert equation, as presented in section 2.1.2, with the length of magnetisation rescaled due to temperature dependency (note that |M| = 1, but parameters in front of M are rescaled). However, such an approach cannot describe a local change of the magnetisation length due to change of its orientation with respect to the external field or anisotropy direction, which might be important for some applications. This problem can be handled partly by assuming space-and time-dependence of the magnetisation length, and adding an additional PDE describing its evolution. In this case, such PDE should 4 Quasistatic models describe systems, evolution of which is slow enough for the system to remain in internal equilibrium. In the case of micromagnetism, quasistatic models are simply static equations (i.e. without time derivative of magnetisation), however with time-dependent boundary conditions. be either derived from the atomistic equation or compiled based on empirical observations.
The last option for the continuum region is the stochastic Landau-Lifshitz-Bloch equation [19,13]. In this approach, due to the stochastic nature and the structure of the dynamic equation, local changes of the expected value of magnetisation length are permitted. The advantage of the stochastic LLB is that it can naturally describe domain walls with nonconstant magnetisation length [13]. However, in other areas of physics, such as modelling the mechanical behaviour of solids, stochastic descriptions of materials are not typically used at the macroscopic scale [31], since stochastic effects diminish with the system size.

An upscaling-based approach
In subsection 2.4, various ways of selecting a macroscopic model and their limitations were mentioned. In this section, an alternative strategy based on HMM is proposed to compute approximate upscaled macroscopic quantities (statistical averages of the atomistic spins given by (2.1)). The approach should be seen as complementary to the domain partitioning approach that is described in subsection 2.3, whereby the continuum equations (needed for the domain partitioning) are obtained computationally. Before proceeding with the description of the method, few mathematical tools and notations are presented.
3.0.1. Averaging kernels. A key ingredient of the HMM is a compression step, which typically uses some general purpose averaging kernels. To illustrate the idea, f ε (t) = f (t/ε) is assumed, where f is a 1-periodic integrable function. A local average of this function can be computed by selecting a small region I η = [−η/2, η/2] at first and then averaging f ε over I η afterwards. This operation, however, results in the following error: This averaging can be done more accurately using the space of averaging kernels K p,q which consists of symmetric functions K, such that K(t) = K(−t), K is compactly supported in [−1/2, 1/2], and d (q+1) K dt q+1 ∈ BV (R), where BV is the space of functions with bounded variations in R, and the derivative is understood in the weak sense. Moreover, the parameter p ∈ Z + represents the number of vanishing moments, i.e. K has p vanishing moments if R K(t)t r dt = 1, r = 0, 0, 0 < r ≤ p.
To simplify the notation, a scaled version of the kernel, denoted by K η (t) := η −1 K(t/η), and the shorthand notation are also introduced. Using a kernel K ∈ K p,q , it is possible to achieve higher order convergence rates and improve the first order approximation (3.1), see e.g. [5,6]  for theoretical results where the following inequality can also be found: In Figure 3.1, the decay of the error for kernels with different p, q pairings is shown.
In the remaining part of this section, the averaging kernels, which are described above, are used to design an upscaling-based strategy for an atomistic-continuum coupling.
3.0.2. HMM. The problem setting is as follows: the atomistic spins interact with each other with a certain exchange parameter J ij and each spin magnetic moment is subjected to an external field H ε e,i (t) = H e,i (t, t/ε), varying at a slow and a fast scale, where ε represents the size of the fast variations. In this setting, a time step small enough to resolve the ε-scale variations of the external field is needed to perform an atomistic spin dynamics simulation. Moreover, to compute a statistical average of the particle magnetisation, a large number of samples of an already expensive atomistic problem have to be computed. Note that the above structure for external field H ε e,i (t) is chosen only as an example, and the same methodology would work equally well for other types of external fields, which have a scale separation of some form, with no conceptual change in the algorithm.
On the other hand, an upscaling-based multiscale method aims at reducing the degrees of freedom, in time and space, by introducing a discretised macroscopic magnetisation field which is related to the local spatial and temporal averages of the spin magnetic moments. From a modelling point of view, one additional requirement is to include the effect of the thermal noise on the magnetisation.
To define the variables, the atomistic spin magnetic moments m ε i are introduced. These magnetic moments are located on ((r + ℓ)N + 1) discrete points defined by , where a is the distance between the particles, r, N ∈ Z + are positive integers and ℓ is a non-negative integer. The superscript ε on m ε i indicates the dependence of the spin magnetic moments on the high-frequency variations in the external field H ε e,i . Furthermore, a set of N + 1 macroscopic variables M I located on the macro grid {X I = I(r + ℓ)a} N I=0 are introduced. In Figure 3  For simplicity, it is assumed that the atomistic particles are globally constrained by periodic boundary conditions, which makes the above definition ofM I (t) welldefined also for the endpoints; namely I = 0 and I = N . The local meansM I are formed by taking the local averages (not only in time but also in space) of 2r + 1 particles. To compute M I , a number of realisations ofM I are needed, which results in a reduction of the macroscopic magnetisation length.
The key idea behind the HMM algorithm for finite temperatures is separate calculation of the directions and the lengths of the expected values of spin magnetic moments, which areM I and s(t), respectively. The macroscopic quantities are introduced in the following way: M I := s(t)M I , whereM I is equal to the local average (3.2) but computed at zero temperature 5 , and s(t) is a variable which describes the effect of the temperature in the system.
An algorithm for computingM I for a single spin at zero temperature (namely when a single spin is subjected to a high-frequency field and the local averaging takes place only in time) was proposed in [4]. At zero temperature, s(t) = 1 holds, and therefore M =M, and the algorithm from [4] is presented below (index I is omitted here). 5 AlthoughM I is computed at zero temperature, the directions of spins at macroscopic nodes are still influenced by the temperature in a nonlinear way via boundary conditions for the microscopic problem. This is evident from the finite-temperature HMM algorithm to be introduced shortly.
Algorithm 1: A single spin subjected to a deterministic oscillatory field.
Step 1. Macro problem: the macroscopic model is incomplete and is given by where F is the missing data in the model and M is the macro solution to be computed.
Step 2. Micro problem: to compute F(t a , M) and close the macro problem, the micro problem is solved first where I + τ := (t a , t a + τ /2] and I − τ := [t a − τ /2, t a ), and τ > ε.
Step 3. Upscaling: the unknown parameter F(t a , M) in the macro model is computed by Remark 1. Starting from a spin magnetic moment m ε (0) with unit length, it was shown in [4], that the length of the macroscopic magnetisation M satisfies Namely, the length of the macroscopic magnetisation is equal to one, up to an upscaling error.
Remark 2. In [4], it is shown that when H ε (t) := H(t, t/ε), where H(t, ·) is a 1periodic smooth function, the error between the HMM (algorithm 1) and the exact macroscopic dynamics is given by Macro error Upscaling error where r 1 and r 2 represent the order of accuracies for a macroscopic and a microscopic solver. Here, ∆t and δt are the macroscopic and the microscopic step sizes. The total computational cost (of HMM) to simulate the dynamics until a time T is given by T τ / (∆tδt). Let β := 1/ (q + 3). Then the computational cost to achieve an error tolerance TOL = O(ε 1−β ) becomes On the other hand, a direct numerical simulation (DNS) gives an error (δt/ε) r2 (with a computational cost T /δt). Therefore, the cost to achieve the same error tolerance becomes The parameter β can be made small by taking large values for q, which does not increase the computational cost of HMM. Moreover, It is clear from the above estimates that the computational cost of the HMM decreases when higher order methods are used for the macro and the micro solvers, while cost of the DNS can not get better than O(ε −1 ).
In algorithm 1 above, the damping term appears only in the macro model (but not the micro problem). However, it was shown (numerically and analytically) that the HMM still captures the correct macroscopic dynamics and hence the correct damping effects, see [4] for the numerical analysis and a motivation for the algorithm.
The aim here is to present a generalisation of algorithm 1 to a system of particles interacting at non-zero temperature. To describe the multiscale algorithm, the notations are adopted. Note the differences between H ε det,i and H sto,i . The term H ε det,i is deterministic but oscillatory, while H sto,i is stochastic and includes the filtered external field K τ * H ε e,i (t * ). The notations used for the terms, other than H ε e,i , in To model the macroscopic magnetisation dynamics of a system of magnetic particles at finite temperature, a macro model, a micro model, an upscaling step and a length scaling procedure are required (instead of solving the continuum equation directly). In Figure 3.3, a schematic description of the method is illustrated. The entire multiscale algorithm is presented below.
Algorithm 2: A chain of magnetic particles at nonzero temperature.
Step 1. Macro problem: the macroscopic model is deterministic and given by d dt

5)
where r ′ = r + ℓ, and F I is the missing data in the model and M I is the macro solution to be computed. Moreover, I ′ := {I − 1, I, I + 1} represents the index set of the nearest macro solutions. The computation of s I (t) is described in Step 4 below.
Step 2. Micro problem: to compute F I (t a , M) and close the macro problem, the micro problem is solved at first d dt m ε Ir ′ +j (t) = β L m ε Ir ′ +j (t) × H ε det,Ir ′ +j (t), t ∈ I ± τ , j = −r, . . . , r, m ε Ir ′ +j (t a ) =m(x Ir ′ +j ), (3.6) m ε Ir ′ −r (t) =m(x Ir ′ −r ), m ε Ir ′ +r (t) =m(x Ir ′ +r ),  The input-output relation for the multiscale algorithm to approximate the macroscopic dynamics at finite temperature. One step of the method (assuming a suitable time-stepping) is depicted. If the macro solver is explicit, the steps are followed in the precise order starting from a given M(t). If the macro solver is implicit, the steps are repeated until the solution M(t + ∆t) (before the length scaling step) converges. Once the convergence is achieved, the length scaling procedure takes place.
Step 3. Upscaling: the parameter F I (t a , M I ′ ) is computed by where η = 2ra.
Step 4. Length scaling: the macro solution M I (t) is replaced by s I (t)M I (t)/|M I (t)|, where s I (t a ) is computed by solving the following stochastic LLG equation for j = −r, . . . , r and t ∈ (t a , t a + τ f ] where τ f > τ r , and τ r is the time it takes to reach the thermal equilibrium. Then, with η = 2ra, the following is computed Remark 3. Note that the parameter s I (t) is allowed to vary in time and space. This approach is more general than the mean-field approximation (MFA) since spatially nonuniform magnetisation lengths are automatically captured here, whereas MFA assumes a uniform length everywhere.
Remark 4. The fact that the HMM algorithm, which is described above, uses local averages (in time and space) of the microscopic solution, makes it possible to upscale the influence of other kind of microscopic variations such as defects (modelled by J ij = 0 for certain i, j pairs) or fast variations in the material properties (modelled by rapid spatial change of J ij or p a ). The resulting macroscopic averages, however, need to be smooth enough so that the macro discretisation resolves it. If not, the atomistic picture needs to be retained, and the partitioned domain approach would be more appropriate.
Remark 5. In equation (3.9), the usage of kernels K η (centred at x Ir ′ ) for spatial averaging is important. At finite temperatures it holds that s I (t) < 1. However, since the boundary conditions for m in (3.8) are fixed, the length of the temporal average of the boundary atom is always 1. The usage of a kernel resolves this issue, as K η vanishes towards the boundaries of the micro domain, see e.g. Figure 3.1. A more secure way of reducing the boundary error is to use a damping band near the boundary, see the equation (2.13), where a damping layer is added to reduce the high frequency reflections from the boundary. However, these damping bands did not seem to be necessary for the simulations presented in the next section.

Numerical results and discussions
4.1. Domain partitioning.

4.1.1.
Damping region at the atomistic-continuum boundary. In [28], a special damping region for domain partitioning multiscale method was suggested. This damping region reduces reflections of high-frequency waves from the atomistic-continuum interface, thereby reduces the numerical error of the solution within the atomistic region of the mutiscale system. The damping band acts as a low-pass filter.
This concept was previously tested only in deterministic cases [28,27]. In the case when the atomistic region is stochastic, the main requirement for the damping band is not to affect the thermodynamic characteristics of atoms. In the spin dynamics approach, each atom is coupled to the heat bath, which allows using a rather simple damping for a certain layer of atoms, while preserving correct thermodynamics in neighbouring layers. A numerical experiment designed to demonstrate this is presented below.
The dynamic behaviour of the multiscale system, in which the atomistic region is embedded into the continuum region, is simulated. The details regarding the multiscale coupling method and the damping region at the interface can be found in [27]. Material parameters (µ, J, K a , p a , λ) that were used in the simulation are summarised in Table 2.1, material 2. Atomistic region was located in the centre of the geometry, the entire length of which was L. Spatial step of ∆x = 4a was used in the continuum model, where a = 0.3636 · 10 −9 is the interatomic spacing, as well as the Neumann boundary conditions. The atomistic model consisted of N atoms, N D out of which were numerically damped atoms (see [27]), i.e. N D /2 damped atoms at left and right edges of the atomistic region. Numerical damping parameter was selected to be g D γ 0.5 = 41.53 and the averaging window width was s A = 3∆x. The parameters are summarised in Table 4.1, set 1. The initial state of the system was uniform alignment of all spins. Since this system had only exchange In Figure 4.1, the evolution of cumulative distribution function of interatomic interaction energy (ε i ) is shown. When system is at the equilibrium, the cumulative distribution of energies should be close to the Boltzmann distribution, W (e) = 1 − exp (−e/ (k B T )). As seen from the results, the non-equilibrium initial configuration evolved into equilibrium configuration with the correct energy distribution. This example illustrates that the damping band does not affect the thermodynamic characteristics of the spin system, while reducing wave reflections form the interface. Since it is possible to adjust the parameters of the damping band independently of the geometry and control the error resulting from wave reflections, such damping region can be very useful for coupling dynamic atomistic region to quasistatic continuum region as well.
4.1.2. 1D system, switching of the uniform external field. One of the tests of consistency between atomistic and continuum regions is to check for an identical response of the regions to a non-stationary external field under equilibrium conditions. To illustrate this, a multiscale system, in which the atomistic region is embedded into the continuum region, is again considered. The system was subjected to an initial applied external field H e = 0. When the system reached equilibrium, tγ = 10, the field was changed to H e = H e e z , H e = 2, which initiated a precessional motion of the spin magnetic moments. Material parameters (µ, J, K a , p a , λ) that were used in the simulation are summarised in Table 2.1, material 2, while other simulation parameters are summarised in Table 4.1, set 3. The Neumann boundary conditions were used in the continuum region.
In Figure 4.2, the behaviour of the magnetisation of the system at thermodynamic equilibrium is shown. It can be seen that the spin moments within the atomistic region rotate with the same angular speed as the continuum region magnetisation, since the average m z component of the atomistic spins is always close to the continuum value.
Although this is a trivial example, it illustrates that the basic requirement of consistency between models, which are used in the multiscale system, can be fulfilled in certain cases without extension of the continuum description. In this example, parameters of the continuum model corresponded to 0 K, and consistency was observed since low temperature was imposed. This is not always the case, as demonstrated in the next example. As was already discussed, the continuum model should be selected depending on the problem to which the multiscale method is applied and sometimes even the basic models provide the solution with an acceptable numerical error. 4.1.3. 1D system, uniform magnetisation gradient. The domain partitioning approach also has limitations. These limitations are not related to the coupling between atomistic and continuum regions, but to the nature of the domain partitioning methodology itself. The problem concerns finite-temperature dynamics, where the macroscopic model is dynamic.
Since the atomistic region is stochastic, excitations of different length scale appear in it. Although the probability of large scale excitation, i.e. several spin atomic moments moving in approximately uniform direction, is small, such excitation creates long wave length spin waves, which propagate within the atomistic region. These waves can propagate into the continuum region. In this case, when the ensemble average solution is calculated, a continuum magnetisation with decreased length is obtained. Such long wave length waves cannot be filtered out by the damping band, since it would make modelling dynamics pointless: the single purpose of the damping band is to filter out the waves which cannot be represented by the continuum. In this setting, it is impossible to distinguish between "useful" long-range excitations, which have to modelled, and "unphysical" ones.
To demonstrate this, a multiscale system, in which the atomistic region is embedded into the continuum region, is considered. Magnetisation at the edges of the continuum domain was fixed, i.e. the stationary Dirichlet boundary conditions were used in the continuum region. Material parameters (µ, J, K a , p a , λ) that were used in the simulation are summarised in Table 2.1, material 2, while other simulation parameters are summarised in Table 4.1, set 2. The evolution of the system was simulated until the equilibrium state, which corresponded to tγ = 10 4 , was reached. The final result is obtained by calculating the ensemble average.
In Figure 4.3, the state of the system at the equilibrium is shown. It can be seen that the expected value of the magnetisation length 6 in the continuum is decreased, see Figure 4.3b, due to waves which are generated within the atomistic region and propagate into the continuum region. The consequences of this issue can be less profound if larger mesh difference is used. Therefore, HMM should not suffer from this issue due to a relatively large separation of scales. However, this problem can be completely eliminated only by using a quasistatic continuum model. In this case, all waves originating from the atomistic region will be either damped within the damping band region or reflected from the interface.

Upscaling via HMM.
In what follows, three different applications of the HMM for non-zero temperatures (Algorithm 2), presented in Section 3, as well as its limitations are described. The first example deals with a single magnetic particle which interacts with a high-frequency external field. As a single particle is studied, i.e. no particle interaction is involved, and the spatial averaging in the algorithm is not used. The second example concerns the behaviour of a number of interacting particles. In this example, the number of atomistic particles is chosen equal to the number of macroscopic nodes. The microscopic model (3.6) as well as the stochastic LLG equation (3.8) are not solved locally (as described above) but globally over the entire domain, and the upscaling step (3.7) is performed only in time. The third example corresponds to the case of interacting particles and disparate spatial scales for micro and macro models. Hence all parts of the algorithm (as described in the Section 3) are utilised, i.e. multiscaling both in time and space is considered. In the last part, the limitations of the HMM are illustrated. 4.2.1. HMM: Single spin. In this first example, a single magnetic particle is subjected to a high frequency field which is nonzero in the z-direction. Namely, (4.1) H ε e = (1 + cos(0.43t) + cos(2πt/ε) 2 )e z , while the parameters in (2.1) are set to K a = J ij = 0. Initially, the direction of the spin magnetic moment is set to 1 √ 3 (e x + e y + e z ). When thermal equilibrium is reached, the solution has to be in the direction of the applied field (z-direction). The external field is varying rapidly over time, and hence it influences the evolution of the length of magnetic moment 7 . In Figure 4.4, it can be seen that the HMM captures the evolution computed by a direct numerical simulation of the spin length. In the same Figure, the HMM solution at zero temperature (with s(t) = 1 or D = 0) is shown to emphasise that if the length scaling part (3.8) of the HMM is not used, then the method does not follow the correct dynamics, see Figure 4.4c. For these simulations, an implicit mid-point rule with 35 macroscopic time steps is used but other numerical integration methods may also be used without any conceptual change in the algorithm, see e.g. [12,15,14] for a review of numerical methods in micromagnetism. The macroscopic magnetisation is evolved until t = 5, and the micro problems are solved over an interval of size τ = 5ε. The physical parameters are chosen to be β L = 1, α L = 10 and D = 0.2. I=0 is equal to that of the microscopic variables {m ε i } 9 i=0 , and the multiscaling is performed only in time. Namely, the micro problem (3.6), and the stochastic LLG equation (3.8) are solved globally (with the periodic boundary conditions), and the upscaling step (3.7) uses averaging only in time. At first, the algorithm was applied to a case with an external field, which is aligned along the z-direction: H ε e,i (t) = (1 + cos(0.43t) + sin(0.73t) + cos(2πt/ε) 2 )e z , i = 0, . . . , 9. This field is uniform in space and oscillatory in time. Similar to the example of a single spin, the HMM solution was compared to a direct numerical simulation. As the external field (4.2) is uniform in space, all the particles have the same statistical behaviour. Figure 4.5 shows that the magnetisation length is accurately captured for all the components. A decrease in the magnitude of the external field leads to a decrease in the magnetisation length as the behaviour of the system becomes more dominated by the thermal noise (leading to a disordered system). Note the decrease in the magnetisation length in Figure 4.5c as the external field, shown in Figure 4.5d, is decreasing. The notation M ν (t, q) in the figure is used for the HMM  solution at the spatial point x = q, where ν can be x, y, or z coordinate. For these simulations, an implicit mid-point rule with 35 macroscopic time steps is used at the macroscopic scale and the final simulation time is t = 5.1. The Heun method [30] is used for timestepping in the micro problem. An interval of size τ = 5ε is used for the micro problem. The physical parameters are chosen to be β L = 1, α L = 10, and D = 0.2.
When the external field is nonuniform in space, the magnetisation length should be nonuniform as well. Figure 4.6b clearly shows that the magnetisation converges to different equilibrium states for two different points in space (x = 0.2 and x = 0.6). The corresponding external fields at these two points are shown in 4.6c. In these simulations 20 macroscopic time steps are used, and the final simulation time is   Figure 4.5. The HMM is shown to capture the correct magnetisation dynamics for this example of a non-uniform external field. Note that the expected value of the exact reference solution (direct numerical solution) is computed by averaging 70 replicas of the solution. In Figure 4.6a, a small deviation between the HMM solution and the reference solution is observed. This is mainly due to the fact that only 70 replicas are used while approximating the exact expected value. These deviations decrease upon computing the expected value using a larger number of copies. 4.2.3. HMM: Chain of spins, multiscaling in time and space. In this example, the case of a chain of magnetic particles, which are subjected to the external field (4.2), is considered. The atomistic particles are defined on a fine grid {x i } 99 i=0 , whereas the macroscopic magnetisations are defined on a coarse grid {X I } 9 I=0 . Moreover, each macroscopic variable corresponds to an average of 11 magnetic moments. Here, the HMM uses upscaling both in time and space.   The number of time steps that were used in the macro model was 35, while other numerical parameter values used for these simulations are η = 0.1 (or r = 5 which means that the macroscopic variables are computed by averaging 2r+1 = 11 spins), ε = 0.001, τ = 5ε, and the physical parameters are β L = 1, α L = 10, D = 0.2, J ij = 1 for |i − j| ≤ 1, and J ij = 0 when |i − j| > 1, and µ = 1. The macroscopic solver uses an implicit mid-point rule with 40 time steps and the dynamics up until t = 5.2 is shown. Similar to previous examples, the Heun method is used as the microscopic solver.

4.2.4.
Limitations of HMM. In subsection 4.2.3, the dynamics of the magnetisation was simulated until t = 5.2. Here, precisely the same example with a longer simulation time, t = 8, is considered. When t ≈ 6 the external field (4.2) oscillates around zero, see Figure 4.8d. In this regime, the system becomes dominated by the noise, as external field is too small to stabilise the system, and a large number of magnetic moments is required to capture the correct statistics, see Figure 2.1. This fact leads to a breakdown of the HMM strategy around t ≈ 6, see Figure 4.8a-c. If the number of particles in the microscopic simulations is large, the cost of the HMM  will not differ from that of a full atomistic simulation. Hence, switching to the full atomistic simulations is preferable, in a short temporal regime around t ≈ 6, and once the system is stable enough, the HMM algorithm may be employed again.

Conclusions
In this article, application of atomistic-continuum multiscale computational strategies to modelling finite-temperature magnetisation dynamics was discussed. There is a number of issues that arise due to additional noise term appearing at the atomistic scale, which takes into account non-zero temperature condition. In particular, when the ensemble average behaviour of the atomistic model (ASD) is obtained, the result cannot be accurately described using simple continuum models, as they rely on different restrictive assumptions regarding magnetisation length. This inconsistency creates problems for coupling techniques.
Two major coupling techniques were discussed in this article: domain partitioning and upscaling via HMM. These techniques aim at modelling somewhat different cases, but the approaches are complementary to each other. Domain partitioning is suitable when the difference between spatial scales, which are relevant for modelling certain phenomena, is not large. Moreover, as it is explicitly required to use  a specific continuum model at the continuum region in the case of domain partitioning, the application of this approach is limited to relatively low temperatures, since at high temperatures, non-linear dependency of the continuum parameters on magnetisation gradient, which is not taken into account in established models, plays a significant role. The selection of a continuum model for the case of domain partitioning was further discussed in Section 2.4.
It was shown that by using a damping band it is possible to filter out some high-frequency noise at the atomistic-continuum interface. However, by design, this damping band cannot filter out low frequency waves, which originate in the atomistic region and some of which are the result of thermal fluctuations. This leads to inhomogeneous in space reduction of magnetisation length in the continuum region, when ensemble averages are calculated. Such behaviour is certainly unphysical and is the major problem for domain partitioning approach when the continuum region is dynamic. Depending on the magnitude of the temperature and timescales involved in the modelling, this artifact can be neglected, but it fully disappears only when the continuum model is quasistatic.
Upscaling-based approaches are well-suited for problems with large separation of scales. Such a strategy is proposed under the HMM framework. A microscale model, describing the dynamics of spin magnetic moments in spatially confined regions, is used to obtain parameters needed to close an initially incomplete macroscale model. The method uses a two-way coupling between a macroscopic and a microscopic model. The macroscopic solver uses large step sizes in time and space, leading to a significant reduction of the degrees of freedom to be computed in comparison to a full atomistic spin dynamics simulation. An analysis for the computational cost is given in Remark 2 for a single spin at zero-temperature. It is shown that the HMM requires a much lower computational cost in comparison to a direct numerical simulation. Moreover, similar arguments can be safely made when interactions between the atomistic particles are present. The main computational gain comes from the fact that the atomistic model is solved only locally both in time and space, while a direct numerical simulation entails a full resolution of the microscopic scales over the entire temporal and spatial domain. The method has been observed to capture the correct magnetisation dynamics in three different examples.
The upscaling-based approach, which is presented here, is complementary to the domain partitioning approach in the sense that the continuum model, which is required for the domain partitioning, can be obtained by the upscaling, rather than using other physically restrictive continuum models. While using such upscaling methods, the limitations should also be taken into account. For systems which are dominated by noise (weak external field and/or anisotropy), a large number of atomistic spins is required for an accurate representation of macroscopic scale quantities. This is a case, where HMM would require a computational cost similar to that of a full atomistic simulation.

6.2.
A chain of interacting magnetic particles subjected to an external magnetic field. When a number of interacting magnetic particles m i are subjected to a uniform external field H e , the total energy 8 is given by In a mean-field approximation, the contributions of individual spins m i to the total energy are separated and for each m i all neighbouring moments are replaced by their statistical averages in the following way: Moreover, a simple case of H e = H e e z and p a = e z is usually considered, which means that all statistical averages of all spins should be aligned along e z . In this case, it is additionally assumed that the statistical averages of all the spins are equal, i.e. m z i = m z j . Using the derivation from the previous section, m z i is obtained: