Modelling long-range interactions in multiscale simulations of ferromagnetic materials

Atomistic-continuum multiscale modelling is becoming an increasingly popular tool for simulating the behaviour of materials due to its computational efficiency and reliable accuracy. In the case of ferromagnetic materials, the atomistic approach handles the dynamics of spin magnetic moments of individual atoms, while the continuum approximations operate with volume-averaged quantities, such as magnetisation. One of the challenges for multiscale models in relation to physics of ferromagnets is the existence of the long-range dipole-dipole interactions between spins. The aim of the present paper is to demonstrate a way of including these interactions into existing atomistic-continuum coupling methods based on the partitioned-domain and the upscaling strategies. This is achieved by modelling the demagnetising field exclusively at the continuum level and coupling it to both scales. Such an approach relies on the atomistic expression for the magnetisation field converging to the continuum expression when the interatomic spacing approaches zero, which is demonstrated in this paper.


Introduction
There are multiple ways of describing the physics of magnetic materials. At the smallest scale, the spin and orbital movements of electrons are modelled by electronic structure calculations. At a larger scale, the rapid subatomic variations are averaged out and the interaction of spin magnetic moments of individual atoms is simulated, often by using parametrised interactions obtained from a smaller scale. The interaction of atomic spins is described by a system of coupled non-linear ordinary differential equations (ODEs). At the macroscopic scale, 1 non-linear partial differential equations (PDEs) are used to describe the evolution of volume-averaged quantities. The choice of a computational approach depends not only on the scale of application, but also on the required computational efficiency. The atomistic models, although relatively accurate, are prohibitively expensive to solve, whereas continuum models are computationally efficient but may lack certain accuracy.
In contrast to targeting a single scale, multiscale modelling strategies potentially offer the accuracy of atomistic models combined with the efficiency of the macroscopic models. All multiscale models can be categorised into sequential (one-way coupling) and concurrent (two-way coupling) methods. The concurrent multiscale models, in turn, can be separated into two groups-the partitioned-domain and the hierarchical approaches [26]. The latter is referred to as the upscaling approach in this paper.
In the partitioned-domain approach, the entire physical domain is split in regions represented by the atomistic and the continuum models, separated by an explicit interface. The information exchange takes place at the interface and the major challenges for these models are handling non-local atomistic interactions and averaging fast atomistic variations at the interface. In the upscaling approach [1], atomistic models are solved locally to find unknown macroscopic quantities in an initially incomplete macroscopic model. The macromodel is then evolved over the entire computational domain. The upscaling strategy, which is considered in this paper, follows the general framework of the heterogeneous multiscale method (HMM) and uses a two-way coupling between the atomistic and macromodels, where the atomistic simulations use the macrodata as initial and boundary conditions, while the macromodel uses the information coming from the local computations of the atomistic model.
The domain partitioning approach is intended for studying the cases of an interaction of a magnetic structure, e.g. a domain wall or a skyrmion, with an isolated heterogeneity, e.g. a crystallographic defect. In this case, the region of interest is modelled using the atomistic approach, while in the rest of the computational region, the continuum model is used. Such approach relies on the continuum model to be welldefined, i.e. derivable from the atomistic model up to a small course-graining error, which might be neglected far from the region of interest. This is the case for relatively small temperatures and for homogeneous materials. The upscaling approach, on the other hand, is intended for cases when the material microstructure is heterogeneous, but representable using periodically stacked representative volume elements (RVEs). In this case, the continuum model is not well-defined and must be obtained by upscaling. The upscaling strategy is also applicable in other cases, when the continuum model is not well-defined, e.g. magnetic structures under a high temperature or subjected to a high-frequency external fields. In the rest of the introduction, a short overview of applications of the partitioned-domain and the upscaling approaches in relation to multiscale problems for ferromagnetic materials is given.
Construction of multiscale models for magnetic materials is a rapidly developing field and a number of partitioned-domain techniques have been proposed in the past [3,15,19,20,24,25]. For the overview and the comparison of various methods, the reader is referred to exhaustive review articles [21,22] discussing partitioneddomain methods in general and [18] discussing an application of multiscale models to magnetic materials in particular. As mentioned above, a major challenge for the partitioned-domain approach is constructing an interface without introducing surplus artefacts into simulations. In [24], the problem of high-frequency wave reflections from the atomistic-continuum interface has been addressed by introducing additional numerical damping, while in [25], a way of handling non-local interatomic interactions at the interface, by introducing a transition zone with partially coarse-grained interactions, has been suggested.
In terms of upscaling approaches, a way of constructing a macroscopic model of ferromagnetic materials, which is fully coupled to an atomistic model, has recently been reported in [7], where an analysis of the dynamics of a single particle and a chain of particles subjected to a high-frequency external field has been given. The extension of the method to problems at elevated temperatures has been addressed in [6]. In the case of a non-zero temperature, the macroscopic magnetisation vector field, which is the volume average of atomistic spin magnetic moments, has temperature-dependent length. In [6], it has been shown that the upscaling method captures accurately the reduced magnetisation length at the macroscopic scale.
In the case of modelling magnetic materials, there is an additional challenge for multiscale models-the existence of the long-range dipole-dipole interactions between the atomic spins [2]. These interactions cannot be handled in the same way as the short-range interatomic interactions, since this would require unreasonably large padding and/or transition zones in the partitioned-domain approach and unreasonably large microscopic domains in the upscaling approach. Such treatment of long-range interactions diminishes all advantages of multiscale approaches. Therefore, these interactions should be handled in a conceptually different way-using a continuum approach. To the best knowledge of the authors, this idea originates from [15,16], where it has been conceptually introduced for the partitioned-domain approach in application to the magnetisation dynamics. The aim of this paper is twofold: (a) to present a way of incorporating the long-range interactions into the latest partitioned-domain approach based on [25], which has a number of distinctive features compared with the initial partitioned-domain approaches [15,16], and (b) to incorporate the long-range interactions into the upscaling formalism for magnetisation dynamics, which has been recently proposed in [6,7]. This paper is organised as follows. In Section 2, mathematical models governing the behaviour at the atomistic and the continuum scales are described and a convergence study is carried out to quantify the approximation errors in relation to modelling the long-range interactions. In Section 3, the multiscale models based on the partitioned-domain and the HMM frameworks are presented. Finally, numerical results are provided in Section 4 to demonstrate the accuracy of the proposed methods.

Atomistic spin dynamics
At the atomistic scale, the mathematical model is the atomistic Landau-Lifshitz-Gilbert equation [8,13], which is given by

2)
where γ is the gyromagnetic ratio, λ is the phenomenological 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 tensor 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. Finally, the term H a,i is a demagnetising field, which originates from long-range dipole-dipole interactions between spin magnetic moments, and is given by [2] where r ij is the vector connecting atoms i and j , and V a is the volume occupied by a single atom. 2 Parameters λ, μ, J ij and K a can be computed from electronic structure calculations [12] and are considered to be constant for a given material, and μ 0 is the permeability of free space.

Continuum models for magnetisation dynamics
In this section, two different continuum models are presented. In the first part, a wellknown continuum model from the classical micromagnetic theory is provided. It is followed by an alternative continuum model based on the upscaling approach.

Continuum model from classical micromagnetism
At the continuum scale, the magnetisation dynamics is modelled by the following non-linear partial differential equation [2,9]: where M is the normalised magnetisation field and β L and α L are the same coefficients as used in (2.1). At zero temperature, exchange tensorial 3 parameter A e can be obtained directly from the atomistic parameters: where r ij is the vector connecting atoms i and j . The sum is evaluated over all atoms with which atom i interacts. In (2.9), tensor A e is assumed to be spatially constant. Since the anisotropy term is local, the same anisotropy tensor K a is used in the continuum and the atomistic equations. The demagnetising field is denoted as a vector field H c (x) for all x ∈ in , where in ⊂ R 3 defines the interior region of a given magnetic material. The outer region is defined as out = R 3 / in . This field is given by where U is the solution of the following PDE [2]: Here, the coefficient c L is introduced to preserve the physical meaning of the magnetisation as the 'density' of spin magnetic moments, while the continuum equation is based on the normalised magnetisation, i.e. |M| = 1.
The continuum magnetisation field M is equal to the normalised ensemble average of the volume average of the atomic spin magnetic moments m i . The normalisation is introduced due to the nature of the LLG equation, as the unit-length vectors are usually used in the formulations. In general, such atomistic-continuum transition introduces an error to the solution, which is dependent on interatomic spacing and on the magnetisation gradient [25].
It must be noted that at finite temperatures, the continuum model must be modified. These modifications differ depending on the approach and are discussed in [6]. However, it has been shown that even with the modifications, the continuum model cannot approach the atomistic model with a predefined accuracy at finite temperatures, i.e. there is always a finite temperature-dependent error. This is one of the reasons for introducing an alternative continuum model based on upscaling.

A continuum model based on upscaling
The basic idea behind upscaling approaches is to start by assuming a macromodel, in which certain quantities are unknown and must be obtained from a given microscopic model. The form of the macromodel usually requires some knowledge about the physical laws that govern the evolution of macroscopic variables. A macromodel in the form of has been proposed and analysed in [7]. In this macromodel, term F is an unknown quantity, which is then upscaled using the local microscopic equation. On the other hand, macroscopic variable M is obtained by using a suitable time discretisation. While designing such upscaling strategies, one important issue is the synchronisation of the microproblems using the coarse-scale variables, which is achieved by assigning suitable initial and boundary conditions for the microproblems. Note that the longrange field does not appear explicitly in the macromodel (2.13), but the macromodel should capture the effect of the long-range interactions. Later, in Section 3.2, it is shown that this can be achieved by including the long-range field in the microscopic models without computing the computationally expensive atomistic dipole-dipole interactions.

Quantification of the errors associated with approximations of the demagnetising field
It is well-known that the atomistic demagnetisation field (2.6) is consistent with the continuum demagnetisation field (2.10). This can be seen from the representation of the continuum PDE for U using the Green's function [2], subsequent application of the divergence theorem and substitution into expression for H c , which gives The atomistic expression (2.6) can be seen as a discretisation of the integral expression for H c .
Since construction of the multiscale models with the long-range interaction relies on the consistency between the atomistic and the continuum expressions for the demagnetisation field, the aim of this section is to demonstrate this consistency numerically and to investigate the convergence rates of the errors related to (a) the approximation error, quantifying the error between the continuum field H c and the atomistic field H a,i , and (b) the geometric error, which is due to neglecting the far-field particles in the computation of H a,i . For the sake of comparison, a brief derivation of a 1D solution for the magnetic potential is provided below.

Analytical solution for the magnetic potential in 1D
Consider domain in = [0, R] × R 2 , which is bounded in x 1 direction and is infinite in x 2 and x 3 directions. Assume also a magnetisation vector field M = M(x 1 ), which is a function of the first coordinate only. Then it follows that This is a one-dimensional problem, i.e. U in and U out depend only on x 1 , since magnetisation M depends only on x 1 . Hence, the Laplace operator reduces to the second derivative by x 1 and the following solution is obtained: For U out to be bounded as x 1 → ∞, it is necessary that ∂ x 1 U out (x 1 ) = b 1 = 0. This implies that U out is a constant. Moreover, since there is a jump in the derivative of U at interfaces x 1 = 0 and x 1 = R, one obtains Substituting these equalities into equation (2.15) gives c 1 = 0. Hence, U in becomes This leads to It must be noted that in the derivation above, U out does not necessarily decay to zero, but remains bounded for all x 1 . However, the derivative ∂ x 1 U out must decay to zero at infinity for a well-defined solution. This is related to the one-dimensional nature of the problem. Namely, the Green's function in 1D does not decay to zero. In a three-dimensional problem, the solution itself also must to decay at infinity.
Analytical solutions are also available in a few specific three-dimensional domains with a uniform magnetisation over the domain [2].

A qualitative comparison between H c and H a,i
In what follows, the aim is to demonstrate that terms H a,i and H c , which are given by (2.6) and (2.10), respectively, match qualitatively. Domain in = [0, 1] × R 2 is considered and magnetisation M = M(x 1 ) is assumed to be dependent on the first coordinate only. From (2.17), the continuum solution is given by For the purpose of this example, c L = 1 is assumed. To compute the atomistic field H a,i via equation (2.6), first, a truncation of domain in is required. Therefore, the following hyper-rectangle is defined: Next, the atomistic lattice is defined as a uniform discretisation of R , where the step size of the discretisation is chosen to be equal to the interatomic distance a. Namely, The atomistic field is computed over the lattice (2.18), which contains (N 1 + 1)(N 2 +1) 2 number of atoms. It is also assumed that the magnetic moments are given by Figure 1 confirms a qualitative match between the atomistic expression (2.6) for H a,i and the continuum equation (2.10) for H c . It is evident, from equation (2.6), that the computed atomistic field H a,i depends on the truncation length R, as well as the choice of the atomic distance a. This dependency results in a small difference between the atomistic and the continuum fields depicted in Fig. 1. The next subsection focuses on studying the convergence rates for the errors with respect to a (the approximation error) and R (the geometric error). Note that the x 2 and x 3 components of the continuum field are equal to zero and hence are not included in the figure.

The approximation error
Given a magnetic body ⊂ R 3 filled with a number of atoms, the continuum field H c (x i ), x i ∈ is obtained as the limit of the atomistic field H a,i when a → 0, i.e.
In particular, a sequence of atomic distances given by a k = 2 k a min , where k = 0, 1, . . . , 5 and a min = 0.00625 is used, and the following differences are recorded.
The convergence rate can then be obtained by computing The values of E k and s k are shown in Table 1, and a first-order convergence rate with respect to the interatomic spacing for the approximation error is observed for this specific example.

The geometric error
When computing the atomistic field H a,i , one often has to deal with large computational geometries relative to the atomic distance a. When the size of the magnetic body is large, computation of H a,i via the summation formula (2.6) is unreasonably expensive. One strategy can be ignoring the atoms, which are located far from atom i, implying a truncation of the computational geometry. In this subsection, the goal is to understand the decay of the error, which arises from truncating the computational geometry. For this, the domain [0, 1] 2 × [−R, R] is uniformly discretised using the interatomic distance a = 0.1. For the computations, it is assumed that The convergence is studied for R k = 2 k , k = 0, 1, . . . , 4, and the errors are recorded. To find the rate s of the geometric error, which is assumed to be O(R −s ), the following ratios are computed.
The results summarised in Table 2 show a second-order convergence rate for the geometric error.

Energy-based and force-based coupling
In the multiscale partitioned-domain coupling approach considered in this paper, the entire computational region is split into two subregions-the atomistic and the continuum domains. There is a 'sharp' atomistic-continuum interface separating these two regions, as illustrated in Fig. 2. Increase of k corresponds to the increase of the truncation radius, i.e. for k = 0, R = 1, while for k = 3, R = 8. Here, s k is the numerical approximation of the convergence rate with respect to the cutoff radius It is well-known that all partitioned-domain methods can be separated into two conceptually distinct groups-the energy-based coupling and the force-based coupling [26]. In the energy-based coupling, the total energy functional of the system is written and forces or torques are derived from it, in the case of modelling deformation or magnetism, respectively. The continuum and the atomistic equations are discretised in time and are advanced together as a single unified system. Thus, in the absence of damping and when an energy-conserving time-stepping method is used, the total energy of the system is conserved. In the force-based coupling, the continuum and the atomistic regions are advanced separately, while exchanging the boundary conditions via the padding atoms and the interface. This ensures that the correct solution is transferred between the regions; however, the total energy of the system is not well-defined in this case.
The continuum region can be discretised using either the finite-element method (FEM) or the finite-difference method (FDM). FEM has an advantage of creating an interface that is conforming to the atomistic lattice, which allows constructing an energy-conserving coupling method [25]; however, it has a more complex implementation than FDM. In the case of FDM, the coupling requires boundary conditions for the continuum region, which can be obtained from additionally constructed 'padding' nodes, where the solution is obtained by volume-averaging of the atomistic solution.
Most energy-based methods have a disadvantage of having numerical artefacts at the atomistic-continuum interface, which are referred to as 'ghost-forces' in the case of modelling deformation or 'ghost-torques' in the case of magnetism, which emerge due to a non-local interatomic interaction. These artefacts can only be removed by employing complex methods of constructing transition zones at the atomisticcontinuum interface [23,25]. Thus, the energy-based coupling requires modification of atoms close to the interface (into transition atoms) to remove ghost-forces/torques. The force-based coupling, on the other hand, requires construction of pad atoms to provide the boundary conditions for the atomistic region, which is simpler both computationally and in terms of implementation. The disadvantage of the force-based methods is the absence of the well-defined total energy of the system. Thus, the energy-based and the force-based methods have a somewhat different scope. For the discussion of the energy-based vs the force-based methods, the reader is referred to [26].

Including the long-range interactions
In the case of magnetic materials, the long-range dipole-dipole interactions are by definition non-local. However, they cannot be handled in the same way as, for example, the exchange interactions, which can be truncated at a relatively short range. In the energy-based approach, the transition zone that removes the ghost-torques must be of the same width as the interaction distance [25], which makes it impractical to create such a zone. In the force-based approach, the width of the region with padding atoms must also be larger than the interaction distance, which also renders it impractical. Therefore, the only solution is to use a conceptually different handling of these interactions.
The idea of handling the long-range dipole-dipole interactions in multiscale models of magnetic materials is to model these interactions using an exclusively continuum approach. Since the atomistic and the continuum domains occupy different spatial domains, an auxiliary computational mesh that covers the entire physical region must be introduced, as illustrated in Fig. 2. The continuum equation for the magnetic potential (2.12) is then solved on this auxiliary mesh. Within the region of this auxiliary mesh that covers the atomistic region, the atomistic demagnetising field H a,i is equated to the continuum demagnetising field H c , while the magnetisation M in (2.12) is, in turn, equated to the atomistic solution m i , which might require interpolation and/or volume-averaging that is discussed below. This idea relies on the convergence of the atomistic expression for the demagnetising field (2.6) to the continuum expression (2.10) as a → 0.
There are differences, however, depending on whether the approach is energybased or force-based and whether FEM or FDM is used for the continuum region. In the energy-based approach and the FEM-discretised continuum with the finiteelement mesh refined down to the atomistic lattice, Fig. 2a, an auxiliary mesh can be constructed, which contains nodes that exactly coincide with the continuum nodes within the continuum region and that exactly coincide with the atomistic positions within the atomistic region. This ensures that interpolation of the solution and the demagnetising field is avoided. Thus, when an energy-conserving time-stepping method is used, such discretisation is energy-conserving.
In the force-based approach and the FDM-discretised continuum with the structured mesh, Fig. 2b, the auxiliary mesh should also be structured and should be an extension of the continuum mesh to the entire computational domain. In this case, within the region of the auxiliary mesh that covers the atomistic region, the magnetisation M in (2.12) is obtained by volume-averaging of the atomistic spin magnetic moments, while the atomistic demagnetising field H a,i is obtained by interpolation of the continuum demagnetising field H c .
It is also possible to have the force-based approach and the FEM-discretised continuum. In this case, since the conservation of the total energy becomes irrelevant, the interpolation can be used, which means that the mesh used for the solution of the equation for the demagnetising field can be arbitrary. The final combination of the energy-based approach and the FDM-discretised continuum is somewhat strange and probably does not have a practical purpose, as the construction of error-free interface coupling is not straightforward, i.e. the specific interaction between each interface node and surrounding atoms the minimises the ghost-forces/torques must be derived in this case.

The scheme used in the numerical examples
In the numerical examples of this paper, the continuum region is discretised using FDM. The regions are coupled using a variant of the force-based coupling, modified to be used together with the implicit time-stepping.
To provide the boundary conditions for the atomistic region, the padding atoms are constructed, see Fig. 2. The solution at the padding atoms is obtained by bilinear interpolation of the continuum solution with subsequent normalisation. The normalisation is introduced to preserve the length of the spin magnetic moments.
To evolve the continuum solution using an implicit time-stepping method, the continuum mesh is extended to the entire computational domain. The solution at the continuum nodes, which overlap with the atomistic region, is obtained by a normalised weighted average of the atomistic solution inside the box with side x centred at the node, where x is the continuum mesh size. For all atoms inside the box, the weight is assigned as the area of the intersection of the box with side a centred at the atom and the box with side x centred at the node. The normalisation is introduced to preserve the nodal length of the vector field solution.
Furthermore, the auxiliary mesh for solving (2.12) is introduced. It coincides with the extended continuum mesh, which is discussed above. Since the nodes of the meshes coincide, magnetisation M in (2.12) is taken to be equal to the continuum magnetisation.
To reduce the high-frequency wave-reflection from the atomistic-continuum interface, additional numerical damping is added to atoms close to the atomisticcontinuum interface [24,25]. This damping acts as a low-pass filter for the waves travelling from the atomistic region to the continuum region, as the solution is 'attenuated' to an average solution within a certain window. Due to a dispersive nature of the spin waves, the damping is non-linear and depends on time derivative of the solution. The analysis of the dynamics of the damping layer and the exact form of the modification can be found in [24].
Following the force-based coupling methodology, the time stepping is performed separately for the atomistic and the continuum regions. The implicit mid-point method [10] is used to solve the equations in time. Within a particular time step, the continuum region (extended to the entire computational domain) is solved first to obtain the current time-step values, which includes the solution of the equation for the demagnetisation field. This gives the solution at the padding atoms at the current time step. The atomistic region is subsequently solved using the padding atoms as the boundary conditions and the demagnetisation field at the current time step. Finally, the solution at the continuum nodes, which overlap with the atomistic region, is overwritten by volume-averaging of the atomistic solution.

Heterogeneous multiscale methods
Recently, an HMM approach has been formulated in application to multiscale problems arising in micromagnetism. First, in [7], a multiscale method has been proposed to simulate the coarse-scale dynamics of a chain of atomistic spins. The atomistic spins were subjected to a high-frequency external field and a mathematical investigation of the convergence rates in relation to the coupling/upscaling errors, originating from a micro-macro coupling, was given for a simplified setting of a single spin.
At finite temperatures, the atomistic LLG (2.1) also includes a white noise term. The noise term results in fluctuations of the magnetic moment vectors m i . The macroscopic quantities of interests, in this case, are the expected values of the local averages (in space and time) of the magnetic moments. By taking the inner product of the (2.1) with m i , it is easy to see that the length of each individual moment m i (t) is equal to |m i (t)| = 1, ∀t. However, due to the thermal fluctuations, the statistical averages acquire reduced lengths, i.e. |E[m i ](t)| < 1. This has been the major reason for the development of the finite-temperature HMM-based model [6]. In both HMM-based algorithms, for zero and for non-zero temperature, the modelling of the long-range interactions has not been considered, as the main ambition has been to model the local terms and the temperature effects accurately, in [7] and [6], respectively. In Sections 3.2.2 and 3.2.3 of this paper, an extension of the algorithms from [6,7] is presented, after introducing the mathematical tools and notations in Section 3.2.1. In particular, it is demonstrated that the microproblems associated with both multiscale methods must be modified in a suitable way to capture the correct macroscopic dynamics in the presence of the long-range interactions.

Averaging kernels
In this subsection, the basic mathematical tools and notations for the HMM algorithms in Sections 3.2.2 and 3.2.3 are introduced. The HMM algorithms developed in [7] and [6] are based on the notion of upscaling, where a local average of small scale features in the atomistic solution (2.1) is computed and used in a macroscale model. The local averaging takes place in small domains in space and time. In practice, the spatial size of the averaging is comparable with the size of a few interatomic distances, i.e. η = ma, m ∈ Z + . The temporal averaging, however, takes place on a domain of size τ = O(ε), where ε is a time scale, at which the microscopic dynamics undergoes some variations. For averaging, the space K p,q of averaging kernels (weight functions) is introduced. The space K p,q consists of functions K, which have compact support in [−1/2, 1/2], and Applying a kernel K ∈ K p,q to an ε-periodic function f ε (t) = f (t/ε), where f is 1-periodic, with an average defined asf := 1 0 f (s) ds results in the following arbitrarily high convergence rates, see e.g. [4,5], where K τ (·) := 1 τ K(·/τ ) is a scaled kernel, q is the smoothness parameter associated with K, and Note that a constant kernel belongs to the space K 1,−1 , i.e. q = −1, and, therefore, the corresponding error becomes O(ε/τ ), in view of the estimate (3.1). In general, smoother the kernel K (higher q), higher the convergence rates become. Moreover, if τ is an exact integer multiple of ε, then the constant C in the above estimate is zero and the averaging is exact. For a numerical verification of the convergence rate in the above estimate and more general results for non-periodic integrands, see [4,5,11].
These local averaging kernels will be used in the description of the HMM algorithms below.

HMM at zero temperature
To present the numerical method, a 1D case is considered. The extension of the algorithm to higher dimensions is self-evident, and skipped in the exposition. First, it is assumed that the full atomistic system consists of N = (r + )L + 1 number of magnetic moments that are located on a set of discrete points in 1D, i.e.
, where a represents the interatomic distance and r ∈ Z + and ∈ N are two non-negative integers. Moreover, the magnetic moments are supplied with periodic boundary conditions (BCs) and are the solutions of the atomistic LLG model is a high-frequency external field oscillating with the wavelength ε in time. The index i in H ε e,i is to allow for spatially non-uniform external fields. For the exchange coefficient, a nearest-neighbour interaction is assumed, i.e.
The macroscopic variable is defined as the local average of (2r + 1) microscopic magnetic moments. To define the macrovariable, assume that the coarse grid is given by

K η (x I (r+ )+j − x I (r+ ) ) K τ * m I (r+ )+j (t)
where η = (2r + 1)a is the size of the local spatial averaging domain and τ > ε is that of a temporal averaging. From formula (3.2), it is evident that between two consecutive macroscopic points, a total number of magnetic moments is skipped while averaging.

Macromodel
The macromodel takes the form Upscaling The last step is to upscale the quantity F I (t * , MĨ ) in (3.3) by Note that instead of a full atomistic simulation over the entire computational domain, a fewer number of atoms (2r + 1 atoms) are coupled together in the microproblem (3.4). Moreover, the boundary atoms and the initial data of the microproblem are forced to be equal to the coarse-scale variables, to synchronise the microscopic model with macrovariables. In (3.4), the quantity H c is computed by solving the equation (2.12) on the macroscopic grid. Moreover, in the computation of H c , the macrosolutions M are used in the right hand side of (2.12). It is worth mentioning that the damping term is not included in the micromodel, since it is modelled at the macroscopic level, see equation (3.3). This is similar to the HMM algorithm from [7], where the convergence of the macroscopic solutions to the exact coarse-scale solutions has been proved in the absence of the long-range field. The main novelty of the current algorithm is the replacement of the microscopic long-range interaction field with the continuum long-range field H c , which can be efficiently approximated using a standard finite difference/element method on the macroscopic domain. This approach leads to a tremendous gain in computational cost due to the fact that the atomistic computation of the long-range field is avoided, which would otherwise require the atoms in a given microscopic domain to communicate with the atoms located in neighbouring microscopic domains over a large macroscopic geometry.

Remark 1
In principle, the macro, (3.3), and the micro, (3.4), problems can be discretised by any convergent time-stepping method. But if certain discrete conservation properties are required, a special care must be given to the choice of the method. The particular choice of the numerical methods, used for the simulations in this paper, can be found in the numerical results section; see also [14] and the references therein for a review about time stepping methods in micromagnetism.

Remark 2
In general, the macroscopic quantities are much smoother in time and space, as they do not 'see' the variations at atomic scales. Hence, in computations, the macroscopic model (3.3) is discretised using a time step t, which is much larger than a time step δt used for a discretisation of the micromodel (3.4).

HMM at non-zero temperature
An extension of the zero temperature algorithm from [7] to non-zero temperature was introduced in [6]. Modelling the long-range interactions requires yet another set of modifications to the algorithm from [6]. To describe these modifications, let be 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 superscript ε in the term H ε e is to denote that the external field has high-frequency variations.
The model at non-zero temperature requires a modification of the zero temperature model. In particular, an additional step is needed to capture the reduction in the length, which arises from taking statistical averages of atomic moments. The precise algorithm (in the presence of the long-range interactions) is given below.
Macromodel With a slight deviation to the algorithm at zero temperature, the macromodel at nonzero temperatures is given by (3.8) whereĨ := {I − 1, I, I + 1}, F I is the missing data in the model and s I (t)M I (t) is the ultimate macrosolution, which models the coarse-scale dynamics. In particular, M I (t) has unit length, up to an upscaling error, and represents the direction and s I (t) ≤ 1, which is computed below, accounts for the reduction in the magnetisation length.
Upscaling The quantity F I (t * , sĨ MĨ ) is computed by where η = (2r + 1)a. where τ f > τ r , and τ r is the time it takes to reach the thermal equilibrium, andM is defined similarly as in the microproblem. Then, with η = (2r + 1)a, the following is computed

Computation of s I (t )
(3.12) Note that in the zero and the non-zero temperature algorithms, the long-range continuum field appears in the microproblem. Moreover, for the non-zero temperature HMM, it is necessary to include the long-range continuum field in the length scaling procedure, equation (3.11), as well. This is due to the fact that the magnetisation length is also influenced by the long-range field, see e.g. [2] for a mathematical motivation.

Partitioned-domain example: domain wall kinetics in a 2D structure
The advantages of the partitioned-domain mutiscale technique are revealed in cases when the atomistic resolution is required locally, while the rest of the computational domain is homogeneous and can be approximated with sufficient accuracy by the continuum model. One such example is the domain wall kinetics in a material with local defects. In [25], such a problem was considered and the performance and the advantages of the partitioned-domain multiscale technique were demonstrated. However, in [25], the domain wall was created using only the exchange and the anisotropy terms in the LLG equation. In this paper, a similar example is considered, however, in which the domain wall is created by the exchange and the demagnetisation terms. The field-induced movement of the domain wall in the presence of a void in the magnetic structure is investigated. From the physical point of view, the void in the material can correspond to a microcrack of the sample or to impurity atoms.
The domain wall in a material with the 2D (111) fcc stacking of atoms is considered. The material contains a hexagonal void with the side of na, where a is the lattice spacing and n is ranging from 3 to 6. The major effect that is observed in the simulations, is the blocking of the domain wall by the void of size 6a, while for the size of the void up to and including 5a, the domain wall is only slowed down by the void.

Computational setup and model parameters
Since all quantities are considered to be dimensionless, β L = 1 was used. The damping was selected to be α L = 0.1. Atoms were selected to be arranged according to 2D (111) fcc stacking in the xy-plane. Lattice spacing was taken to be a = 1/64. The exchange coefficients were taken to be J ij /μ = (2/3)a −2 , which gives the continuum exchange tensor A e = A e I with A e /μ = 1, where I is the identity tensor. The anisotropy term was not taken into account either, K a = 0. The coefficient that defines the magnitude of the demagnetisation field was taken to be c L = 2π 2 , which gives the approximate width of the domain wall w D = 1. The external field was applied in the x-direction, H e = H x e x , where H x = 5 was assumed.
The width and the height of the computational region were taken to be x L = 4 and y L = 1, respectively. Neumann boundary conditions were used at x = 0, x = x L , y = 0 and y = y L for both LLG and magnetostatic equations. The computational region was partitioned into the atomistic and the continuum subregions. The continuum discretisation step was taken to be x = 4a. The atomistic subregion was located in the centre of the computational region. The width and the height of the atomistic subregion were taken to be x A = y A = 48a = 0.75. The atomistic subregion contained a hexagonal void with a side of na, where n ∈ {3, 4, 5, 6}, the centre of which was located at x = 2 and y = 0.5. Time step t = 10 −2 was used.
Within the atomistic region, the behaviour of atoms close to the interface was modified with the additional numerical damping. Parameters of the damping region, the optimal values of which depend on the width of the damping region and the difference between discretisations, were taken from [25], where the same difference between the discretisations of the regions was used. The damping strength and the width of the averaging window was taken to be g D = 625, s A = 3 x. The width of the damping band was selected to be 16a, as the large width of the damping band ensures that the atomistic solution is not contaminated by wave reflections.
The following initial conditions for the domain wall were used: where x 0 is the position of the centre of the domain wall. The initial position of the centre of the domain wall was taken to be x 0 = 1. To the left of the domain wall, the magnetisation is oriented primarily in +z-direction, while to the right, it is oriented in −z-direction. Within the domain wall, the magnetisation changes from +z-direction to −z-direction and has a non-zero component in +y-direction.
Since the domain wall moves during the simulation, to analyse the results, it is important to obtain from the simulations the exact position of the centre of the domain wall as a function of time. The domain wall centre is defined as the curve, along which M z = 0. Since the computational solution, M, is defined at the grid points, an auxiliary quantity ζ = arccos (M z ) − π/2 is calculated at each grid point and linearly interpolated between the grid points. Thus, for each y = y 0 , the domain wall centre along x-axis, x 0 , is found by solving ζ = 0.

Results
In Fig. 3, the field plot of the z-component of magnetisation for the case of n = 5 is illustrated. It can be seen that when the domain wall approaches the void, the thickness of the domain wall decreases locally. The region of the domain wall that is located in the upper half of the 2D plate slows down, while the the lower part of the domain wall moves past the void (t = 0.8 and t = 1.0). Afterwards, the upper part Finally, an equilibration process is observed (t = 1.8 and t = 2.0). Thus, the void causes oscillations in the structure of the domain wall. Moreover, in the field plots, it can be seen that the void acts as a 'gradient concentrator', i.e. the region around the void creates higher gradient in comparison with the regions further from the void.
From the physical point of view, the energetically preferred states of spin magnetic moments are approximate alignments of spins either along e z or opposite to e z . The intermediate states of spins, which correspond to the domain wall, have higher energy. When the domain wall passes through the void, a region of the wall is absent (due to the void) and, thus, a number of 'high-energy' spins are absent. Therefore, to minimise the total energy, the preferred states of the spin system are such that the void covers spins with the highest energy locally. This explains why the domain wall and the gradient lines in Fig. 3 tend to stick to the void.
The states of spin magnetic moments around the void when the domain wall passes the void are shown in Fig. 4. Although the initial structure of the domain wall is of the Bloch-type, i.e. the spins have zero M x component, it is clearly seen that as the domain wall interacts with the void, spins acquire a non-zero M x component. This is the mechanism by which the domain wall is slowed down by the void in this example. Since the external field is aligned with e x , the torque that acts on the spins with significant M x component is small, which decreases the angular velocity of the spins in the centre of the domain wall and the thereby leads to the decrease of the speed of the domain wall propagation. Moreover, at t = 0.7, it is seen that the upper part of the domain wall is inclined, as opposed to the lower part of the domain wall,

Partitioned-domain example: multiscale modelling error, 1D example
The proposed multiscale technique obviously has a modelling error, which is the error due to the representation of the demagnetisation field using the discretised continuum model. This error is proportional to c L , which can be understood as the parameter defining the magnitude of the demagnetisation field. This is demonstrated in this section using an example of a 1D domain wall moving form the continuum region into the atomistic region.
As was shown above, in the 1D case, H c = −c L M x e x . This means that the demagnetisation field acts similar to the anisotropy, but with the negative sign. Moreover, due to the structure of the LLG equation (2.7), CM can be added to effective field H , where C is an arbitrary constant, without influencing the solution of the LLG equation. This is used in the example below.

Computational setup and model parameters
The same β L and α L as in the example above were used. Atoms were selected to be arranged in a 1D chain. Lattice spacing was taken to be a = 1/64. The exchange coefficients were taken to be J ij /μ = a −2 , which gives the continuum exchange parameter A e /μ = 1. Biaxial anisotropy was used, The introduction of parameter c L into the anisotropy is explained below. The coefficient that defines the magnitude of the demagnetisation field, c L was varied from 0.4 to 12.8. These parameters give the approximate width of the domain wall w D = 1. The external field was applied in z-direction, H e = H z e z , where H z = 1 was taken.
It can be seen that in the above presented setup, the exact continuum solution of the moving domain wall should not be affected by parameter c L , because However, in the computational solution, there will be a numerical error due to different treatments of the demagnetisation field and the other interactions. Also, it should be noted that negative anisotropy in e z -direction was chosen. Although this is not a realistic scenario, selection of a negative coefficient is mathematically allowed. It allows creating a setup, for which the exact continuum solution is independent of c L . The length of the computational region was taken to be x L = 6. Neumann boundary conditions were used at x = 0, x = x L for the LLG equation. The 1D analytical solution of the magnetostatic equation, H c = −c L M x e x , was used at the continuum scale to isolate the modelling error, i.e. not to introduce an error due to numerical solution of (2.10) and (2.12). The computational region was partitioned into the atomistic and the continuum subregions. The continuum discretisation step was taken to be x = 4a. The atomistic subregion was located in the centre of the computational region. The length of the atomistic subregion was taken to be x A = 128a. Time step t = 10 −2 was used. Same g D , s A and the width of the damping band as in example above are used.
The following initial conditions for the domain wall were used: M = e x sin θ cos ϕ + e y cos θ cos ϕ + e z sin ϕ, where x 0 is the position of the centre of the domain wall. The initial position of the centre of the domain wall was taken to be x 0 = 1.5. The simulations were run until t end = 4 and compared at that point. This roughly corresponded to the domain wall being in the centre of the atomistic region. The reference solution with c L = 0 was used. The solution in the atomistic region was used for the error calculation. The error was defined as the L 1 -norm of the difference between solutions, divided by 3N, where N is the number of atoms of the atomistic region.

Results
In Fig. 6, the dependence of the multiscale modelling error on c L is shown. This is the error due to numerical multiscale treatment of the demagnetisation field, compared with the reference case, where the analytical expression for the demagnetisation field is used. The error is proportional to c L .

A chain of magnetic particles
Here, the HMM algorithm at zero temperature is applied to a chain of magnetic particles. It is assumed that N = 101 magnetic moments are located on a one-dimensional Fig. 6 The dependence of the multiscale modelling error on c L for the 1D domain wall propagation example lattice consisting of points {ia} 100 i=0 , where a = 0.01 is the atomic distance and the initial configuration is given by m i (t)| t=0 = cos(2πia)e x + sin(2πia)e y , i = 0, 1, . . . , 100.
Moreover, the high-frequency external field with f ε (t) = 1 + cos(0.43t) + cos 2 (2πt/ε) and ε = 0.01 is used in the simulations. Note that all the atomic particles are under the influence of the same external field and hence the spatial dependency is omitted. The continuum demagnetisation field H c in the microscopic model (3.4), is given exactly by in this one-dimensional setting, cf. the derivations in Section 2.3.1. The parameters β L = α L = μ = 1, and the nearest neighbour exchange interactions with J = 1, see Section 3.2.2, are chosen for the simulations. Moreover, the anisotropy is assumed to be zero. The problem is discretised by a midpoint rule on the macroscopic and microscopic scales with the temporal step sizes t = 0.01 on the macroscale and δt = ε/10 on the microscales. Moreover, the temporal microscopic box is τ = 5ε. On the spatial dimension, 10 macroscopic points are used to describe the magnetisation dynamics whereas the atomistic chain consists of a total of 101 particles. Figure 7 demonstrates the evolution of the macroscopic dynamics and compares it with the solution of the full atomistic model. In the atomistic simulation, it is assumed that H a,i = −m x e x , which is justified by the convergence of the atomistic long-range field to the continuum counterpart as a → 0, see Section 2.3.3. As time increases, the magnetisation vectors are pointing to the e z direction, i.e. the direction of the given external field, and it is shown in the picture that the correct macroscopic dynamics are captured even though the macroscopic discretisation parameters under-resolves the scales of atomistic variations both in time and space. Note that the presence of the long-range field has a clear effect on the dynamics of the magnetisation, cf. the numerical results in [7] in the absence of the long-range field.

Modelling error
The error, which arises from having the long-range field in the modelling, is studied. The error between the HMM solution and the atomistic solution is recorded with varying values for the coefficient c L . The very same numerical parameters as in the previous example are chosen in the simulation. It is observed that the decay rate of the error is O(c L ), see Fig. 8.

A chain of magnetic particles at nonzero temperature
In this section, the aim is to show an example, where the existence of the longrange interaction has an effect over the macroscopic magnetisation length at elevated temperatures. For a numerical study, N = 101 atomistic particles located on a  Fig. 7 The HMM solution, using 10 macroscopic points, is compared with an atomistic simulation using 100 magnetic moments. A good match between the two solutions is observed  Fig. 8 The error between the HMM solution and the full atomistic simulation is depicted for an increasing values for the long-range field. A first-order convergence, in terms of c L , is observed in the simulation one-dimensional lattice consisting of points {ia} 100 i=0 , where a = 0.01 is the atomic distance, are considered. The initial configuration is uniform and given by m i (t)| t=0 = 1 √ 3 e x + e y + e z , i = 0, 1, . . . , 100, equipped with periodic boundary conditions. The anisotropy is assumed to be zero, and the external field is pointing in the x direction, i.e.
The atomistic particles are subjected to a thermal noise with a standard deviation of D = 0.01, and the parameters β L = μ = 1, α L = 10, with the nearest neighbour interactions for the exchange coefficient, similar to previous examples, are chosen for the simulation. In Fig. 9, a fully atomistic simulation (with 100 realisations) with and without the long-range field is compared against the HMM solution up to a final time T = 5. It is shown that the long-range field has an effect on the length of the magnetisation and that the HMM accurately captures this long-range effect. In Fig. 10, the behaviour of the classical mean-field approach is shown. The mean-field approach is a well-known approach to predict the length of the ensemble averages, , see e.g. [2,6] for the derivations of the mean-field formulas, which can be adapted to the present setting in an obvious way. It is observed that unlike HMM, the mean-field approach deviates from the true atomistic simulation in the presence of the longrange interaction. This is due to the fact that the mean-field approach suffers from the mentioned restrictive closure arguments, which do not hold in general since the

Conclusions
This paper demonstrates a way of incorporating the long-range dipole-dipole interactions between the atomistic spin magnetic moments into the existing atomisticcontinuum coupling methods based on the partitioned-domain and the upscaling strategies. This is achieved by modelling the demagnetising field exclusively at the continuum level and coupling the continuum demagnetising field to the atomistic solution. This approximation relies on the atomistic expression for the magnetisation field converging to the continuum expression when the interatomic spacing approaches zero. It has been demonstrated that for both multiscale strategies, the modelling error is O(c L ), where c L is the coefficient defining the magnitude of the demagnetising field. Moreover, the present article includes numerical results addressing the convergence of the atomistic long-range field to the continuum long-range field and geometric errors involved in the atomistic simulations of the demagnetisation field. Within the framework of partitioned-domain methods, it has been discussed how to account for the long-range interactions within the energy-based methods and the force-based methods. In both approaches, an auxiliary continuum mesh is constructed that covers the entire computational domain, overlapping with the atomistic region and coinciding with the continuum mesh of the continuum region. The equations for the demagnetising field are then solved on this auxiliary mesh. Within the atomistic region, the atomistic demagnetising field is taken to be equal to the continuum demagnetising field.
The computational examples of this paper attempted to highlight cases when the proposed multiscale approaches excel in terms of efficiency. The effect of the voidaffected kinetics of the domain wall has been modelled using the partitioned-domain approach. Only a small region around the void has been modelled at the atomistic scale, while the rest of the 2D structure has been modelled using the continuum model with a coarser resolution. The domain wall structure itself has been the result of the demagnetising field. The partitioned-domain methodology allows resolving fine-scale details of the interaction of the domain wall and the void, while replacing the solution far from the void with a close continuum approximation to increase computational efficiency.
The main novelty and the advantage of the proposed algorithm for the upscaling strategy are the avoidance of the long-range atomistic communication between the magnetic moments that are located in different microscopic boxes. Nevertheless, the macroscopic effect of the long-range interactions is captured accurately in the multiscale formalism. This leads to a significant computational gain not only in comparison with a naive computation of the dipole-dipole interactions (which scales quadratically with respect to the number of particles) but also when compared with more efficient multiscale methods such as the fast multipole method [17,19], which is a linear scaling algorithm. The fact that the long-range field is included in the microscopic model using the continuum field H c allows obtaining sub-linear scaling computational costs with respect to the atomistic degrees of freedom. The accuracy of the method is also demonstrated using an example of a chain of magnetic particles as well as examples at elevated temperatures.
The multiscale methods presented in this article can also be applied to problems with moving atomistic regions. One possibility to tackle this problem is to include an adaptive algorithm at the macroscopic scale, which handles the movement of atomistic regions on the fly. One can then resort to the microscopic model, once the atomistic locations are dynamically determined. Creating dynamically changing boundaries of the atomistic regions is challenging only from the programming point of view, as mathematically no additional modifications of the methods are required.

Funding Information Open access funding provided by Mälardalen University.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.