Modelling long-range interactions in multiscale simulations of ferromagnetic materials

TL;DR

This paper introduces a multiscale modelling approach for ferromagnetic materials that effectively incorporates long-range dipole-dipole interactions by modeling the demagnetising field at the continuum level and coupling it with atomistic and continuum models.

Contribution

It presents a novel method to include long-range interactions in atomistic-continuum multiscale models of ferromagnets by modeling the demagnetising field at the continuum level.

Findings

The atomistic magnetisation converges to the continuum magnetisation as interatomic spacing decreases.

The proposed coupling method accurately captures long-range dipole interactions.

The approach enhances multiscale simulations of ferromagnetic materials.

Abstract

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.

Tables

Table 1. Table 2.1. The approximation error. Decrease of k 𝑘 k corresponds to the decrease of the interatomic spacing, i.e. for k = 4 𝑘 4 k=4 , a = 0.1 𝑎 0.1 a=0.1 , while for k = 0 𝑘 0 k=0 , a = 0.00625 𝑎 0.00625 a=0.00625 . Here, s k subscript 𝑠 𝑘 s_{k} is the numerical approximation of the convergence rate with respect to the interatomic spacing.

$k$	$4$	$3$	$2$	$1$	$0$
$E_k$	$0.1605$	$0.1019$	$0.0619$	$0.0324$	$0.0165$
$s_k$	$0.6555$	$0.7173$	$0.9344$	$0.9713$

Table 2. Table 2.2. The geometric error. Increase of k 𝑘 k corresponds to the increase of the truncation radius, i.e. for k = 0 𝑘 0 k=0 , R = 1 𝑅 1 R=1 , while for k = 3 𝑘 3 k=3 , R = 8 𝑅 8 R=8 . Here, s k subscript 𝑠 𝑘 s_{k} is the numerical approximation of the convergence rate with respect to the cutoff radius.

$k$	$0$	$1$	$2$	$3$
$E_k$	$0.5992$	$0.1963$	$0.0538$	$0.0138$
$s_k$	$1.6096$	$1.8663$	$1.9545$

Equations

$$\frac{\mathrm{d}}{\mathrm{d}t}{\mathbf{m}}_i=-{\beta }_{\mathrm{L}}{\mathbf{m}}_i\times {\mathbf{H}}_i-{\alpha }_{\mathrm{L}}{\mathbf{m}}_i\times \left({\mathbf{m}}_i\times {\mathbf{H}}_i\right),\left|{\mathbf{m}}_i\right|=1,$$

$${\beta }_{\mathrm{L}}=\frac{\gamma }{1+{\lambda }^2},{\alpha }_{\mathrm{L}}=\frac{\gamma \lambda }{1+{\lambda }^2},$$

$${\mathbf{H}}_i=\frac{1}{\mu }\left(\sum _j{J}_{ij}{\mathbf{m}}_j\right)+\frac{1}{\mu }{\mathbf{K}}_{\mathrm{a}}\cdot {\mathbf{m}}_i+{\mathbf{H}}_{\mathrm{e}}+{\mathbf{h}}_i+{\mathbf{H}}_{\mathrm{a},i},$$

$$⟨h_{i\rho }\left(t\right)⟩=0,⟨h_{i\rho }\left(t\right)h_{j\nu }\left(s\right)⟩=2D{\delta }_{ij}{\delta }_{\rho \nu }\delta \left(t-s\right),$$

$$D=k_{\mathrm{B}}T\frac{\lambda }{\mu \gamma },$$

$${\mathbf{H}}_{\mathrm{a},i}={\mu }_0\mu \left(-\frac{{\mathbf{m}}_i}{3V_{\mathrm{a}}}+\frac{1}{4\pi }\sum _{r_{ij}\ne 0}\left(\frac{3{\mathbf{m}}_j\cdot {\mathbf{r}}_{ij}{\mathbf{r}}_{ij}}{r_{ij}^5}-\frac{{\mathbf{m}}_j}{r_{ij}^3}\right)\right),r_{ij}=\left|{\mathbf{r}}_{ij}\right|,$$

$$\frac{\partial }{\partial t}\mathbf{M}\left(t,\mathbf{x}\right)=-{\beta }_{\mathrm{L}}\mathbf{M}\times \mathbf{H}-{\alpha }_{\mathrm{L}}\mathbf{M}\times \left(\mathbf{M}\times \mathbf{H}\right),\left|\mathbf{M}\right|=1,$$

$$\mathbf{H}\left(t,\mathbf{x},\mathbf{M}\right)=\frac{1}{\mu }{\mathbf{A}}_{\mathrm{e}}:\nabla \nabla \mathbf{M}+\frac{1}{\mu }{\mathbf{K}}_{\mathrm{a}}\cdot \mathbf{M}+{\mathbf{H}}_{\mathrm{e}}+{\mathbf{H}}_{\mathrm{c}},$$

$${\mathbf{A}}_{\mathrm{e}}=\frac{1}{2}\sum _{j\ne i}{J}_{ij}{\mathbf{r}}_{ij}{\mathbf{r}}_{ij},$$

$${\mathbf{H}}_{\mathrm{c}}=-c_{\mathrm{L}}\nabla U,c_{\mathrm{L}}=\frac{{\mu }_0\mu }{V_{\mathrm{a}}},$$

$$U\left(\mathbf{x}\right):=\{\begin{array}{ll}{U}_{\mathrm{in}}\left(\mathbf{x}\right),& \mathbf{x}\in {\Omega }_{\mathrm{in}},\\ U_{\mathrm{out}}\left(\mathbf{x}\right),& \mathbf{x}\in {\Omega }_{\mathrm{out}},\end{array}$$

$$\left\{\begin{array}[]{lll}\Delta U_{\mathrm{in}}\left(\boldsymbol{x}\right)=\nabla\cdot\boldsymbol{M}\left(\boldsymbol{x}\right),&\text{in }\Omega_{\mathrm{in}},\\ \Delta U_{\mathrm{out}}\left(\boldsymbol{x}\right)=0,&\text{in }\Omega_{\mathrm{out}},\\ U_{\mathrm{in}}=U_{\mathrm{out}},&\text{on }\partial\Omega_{\mathrm{in}},\\ \partial_{n}U_{\mathrm{in}}-\partial_{n}U_{\mathrm{out}}=\boldsymbol{M}\cdot\boldsymbol{n},&\text{on }\partial\Omega_{\mathrm{in}},\\ \left|rU_{\mathrm{out}}\right|\to 0,\quad\left|r^{2}\nabla U_{\mathrm{out}}\right|\to 0,&\text{at }r\to\infty,&\text{where }r=\left|\boldsymbol{x}\right|.\end{array}\right.$$

$$\frac{\partial }{\partial t}\mathbf{M}=-\mathbf{F}\left(t,\mathbf{x},\mathbf{M}\right)-\frac{{\alpha }_{\mathrm{L}}}{{\beta }_{\mathrm{L}}}\mathbf{M}\times \mathbf{F}\left(t,\mathbf{x},\mathbf{M}\right)$$

$${\mathbf{H}}_{\mathrm{c}}\left(\mathbf{r}\right)=\frac{c_{\mathrm{L}}}{4\pi }{\int }_{{\Omega }_{\mathrm{in}}}\left(\frac{3\mathbf{M}\left({\mathbf{r}}^{\prime }\right)\cdot \left(\mathbf{r}-{\mathbf{r}}^{\prime }\right)\left(\mathbf{r}-{\mathbf{r}}^{\prime }\right)}{{\left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|}^5}-\frac{\mathbf{M}\left({\mathbf{r}}^{\prime }\right)}{{\left|\mathbf{r}-{\mathbf{r}}^{\prime }\right|}^3}\right)\mathrm{d}{V}^{\prime }.$$

$$\Delta U_{\mathrm{in}}(\mathbf{x})={\partial }_{x_1}{M}_1(x_1),\Delta U_{\mathrm{out}}(\mathbf{x})=0.$$

$$U_{\mathrm{in}}(x_1)={\int }_0^{x_1}{M}_1(s_1)\mathrm{d}{s}_1+c_1{x}_1+c_0,$$

$$U_{\mathrm{out}}(x_1)=b_0+b_1{x}_1.$$

$${\partial }_{x_1}{U}_{\mathrm{in}}(R)={\partial }_{x_1}{U}_{\mathrm{out}}(R)+M_1(R)=M_1(R),$$

$${\partial }_{x_1}{U}_{\mathrm{in}}(0)={\partial }_{x_1}{U}_{\mathrm{out}}(0)+M_1(0)=M_1(0).$$

$$U_{\mathrm{in}}(x_1)={\int }_0^{x_1}{M}_1(s_1)\mathrm{d}{s}_1+c_0.$$

$${\mathbf{H}}_{\mathrm{c}}(\mathbf{x})=-c_{\mathrm{L}}\nabla U_{\mathrm{in}}(\mathbf{x})=-c_{\mathrm{L}}{M}_1(x_1){\mathbf{e}}_1,\forall \mathbf{x}\in {\Omega }_{\mathrm{in}}.$$

$${\Omega }_R=[0,1]\times [-R/2,R/2{]}^2.$$

$${\Omega }_{R,a}:=\left\{\left(x_1^i=ia,x_2^j=-R/2+ja,x_3^k=-R/2+ka\right)\right\},$$

$$i=0,\dots ,N_1,N_{1}a=1,j,k=0,\dots ,N_2,N_{2}a=R.$$

$${\mathbf{m}}_i=\frac{1}{\sqrt{2}}\left(\cos (2\pi x_1^i){\mathbf{e}}_1+\sin (2\pi x_1^i){\mathbf{e}}_2+{\mathbf{e}}_3\right),\forall i.$$

$$\underset{a\to 0}{\lim }{\mathbf{H}}_{\mathrm{a},i}={\mathbf{H}}_{\mathrm{c}}({\mathbf{x}}_i).$$

$${\mathbf{m}}_i={\mathbf{e}}_1.$$

$$E_k:=\left|{\mathbf{H}}_{\mathrm{a},i}^{a_{k+1}}-{\mathbf{H}}_{\mathrm{a},i}^{a_k}\right|,k=0,1,\dots ,4.$$

$$2^{s_k}=\frac{E_k}{E_{k-1}}.$$

$${\mathbf{m}}_i=\frac{1}{\sqrt{2}}\left(\cos (2\pi x_1^i){\mathbf{e}}_1+\sin (2\pi x_1^i){\mathbf{e}}_2+{\mathbf{e}}_3\right),\forall i.$$

$$E_k:=\left|{\mathbf{H}}_{\mathrm{a},i}^{R_{k+1}}-{\mathbf{H}}_{\mathrm{a},i}^{R_k}\right|,k=0,1,2,3$$

$$2^{-s_k}=\frac{E_k}{E_{k-1}}.$$

$${\int }_{\mathbb{R}}K(t)t^r\mathrm{d}t=\{\begin{array}{ll}1,& r=0,\\ 0,& 0<r\le p.\end{array}$$

$$\left|\left(K_{\tau }\ast f^{\epsilon }\right)(0)-\bar{f}\right|\le C{\left(\frac{\epsilon }{\tau }\right)}^{q+2},\tau >\epsilon ,$$

Full text

1. 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 nonlinear ordinary differential equations (ODEs). At the macroscopic scale, nonlinear 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 [25]. The latter is referred to as upscaling approach in this paper.

In the partitioned-domain approaches, the entire physical domain is split in regions represented by the atomistic and the continuum models with an explicit interface between them. 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 macro model 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 macro models, where the atomistic simulations use the macro data as initial or boundary conditions, while the macro model uses the information coming from local computations of the atomistic model.

The domain partitioning approach is intended for cases when an interaction of a magnetic structure, e.g. a domain wall or a skyrmion, with an isolated heterogeneity, e.g. a crystallographic defect, is studied. In this case, the region of interest in 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 well-defined, 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 at 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 a high-frequency external fields. In what follows, a short overview of applications of the partitioned-domain and upscaling approaches in relation with 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 [16, 19, 3, 11, 23, 24]. For the overview and the comparison of various methods the reader is referred to exhaustive review articles [21, 20] discussing partitioned-domain methods in general and [18] discussing 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 [23], the problem of high-frequency wave reflections from the atomistic-continuum interface has been addressed by introducing additional numerical damping, while in [24], 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 was given. The extension of the method to problems at elevated temperatures was 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 accurately captures 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 long-range dipole-dipole interactions between 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. The aim of this paper is to present an efficient strategy to include the long-range interactions in the partitioned-domain strategy based on [23, 24] and the upscaling formalism developed in [7, 6].

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.

2. Mathematical models at different scales

2.1. Atomistic spin dynamics

At the atomistic scale, the mathematical model is the atomistic Landau-Lifshitz-Gilbert equation [8, 14], which is given by

[TABLE]

where $\gamma$ is the gyromagnetic ratio, $\lambda$ is the phenomenological damping constant, $\boldsymbol{m}_{i}$ is the direction of spin magnetic moment, $\mu$ is the length of spin magnetic moment, $J_{ij}$ are constants of Heisenberg exchange interaction between atoms $i$ and $j$, $\boldsymbol{K}_{\mathrm{a}}$ is the anisotropy tensor and $\boldsymbol{H}_{\mathrm{e}}$ is the external field. Thermal excitations are taken into account by adding a stationary stochastic field with the following statistical properties:

[TABLE]

where $\rho$ and $\nu$ are the Cartesian coordinates of $\boldsymbol{h}_{i}$, $k_{\mathrm{B}}$ is the Boltzmann constant and $T$ is temperature. Finally, the term $\boldsymbol{H}_{\mathrm{a},i}$ is a demagnetising field, which originates from long-range dipole-dipole interactions between spin magnetic moments, and is given by [2]

[TABLE]

where $\boldsymbol{r}_{ij}$ is the vector connecting atoms $i$ and $j$, and $V_{\mathrm{a}}$ is the volume occupied by a single atom111For crystallographic lattice with cubic stacking, $V_{\mathrm{a}}=a^{3}$, where $a$ is the distance between two neighbouring atoms.. Parameters $\lambda$, $\mu$, $J_{ij}$ and $\boldsymbol{K}_{\mathrm{a}}$ can be computed from electronic structure calculations [13] and are considered to be constant for a given material, and $\mu_{0}$ is the permeability of free space.

2.2. Continuum models for magnetisation dynamics

In this section, two different continuum models are presented. In the first part, a well-known continuum model from the classical micromagnetic theory is provided. It is followed by an alternative continuum model based on the upscaling approach.

2.2.1. Continuum model from classical micromagnetism

At the continuum scale, the magnetisation dynamics is modelled by the following nonlinear partial differential equation [2, 9]:

[TABLE]

where $\boldsymbol{M}$ is the normalised magnetisation field and $\beta_{\mathrm{L}}$ and $\alpha_{\mathrm{L}}$ are the same coefficients as used in (2.1). At zero temperature, exchange tensorial222Here, the standard tensor notation is used, where the tensor product of two vectors is denoted as $\boldsymbol{a}\boldsymbol{b}$, which results in a second-order tensor. The double inner product of two second-order tensors is denoted as $\boldsymbol{A}:\boldsymbol{B}=\sum_{i}\sum_{j}A_{ij}B_{ji}$, which results in a scalar. parameter $\boldsymbol{A}_{\mathrm{e}}$ can be obtained directly from the atomistic parameters:

[TABLE]

where $\boldsymbol{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 $\boldsymbol{A}_{\mathrm{e}}$ is assumed to be spatially constant. Since the anisotropy term is local, the same anisotropy tensor $\boldsymbol{K}_{\mathrm{a}}$ is used in the continuum and the atomistic equations.

The demagnetising field is denoted as a vector field $\boldsymbol{H}_{\mathrm{c}}(\boldsymbol{x})$ for all $\boldsymbol{x}\in\Omega_{\mathrm{in}}$, where $\Omega_{\mathrm{in}}\subset\mathbb{R}^{3}$ defines the interior region of a given magnetic material. The outer region is defined as $\Omega_{\mathrm{out}}=\mathbb{R}^{3}/\Omega_{\mathrm{in}}$. This field is given by

[TABLE]

where $U$ is the solution of the following PDE [2]:

[TABLE]

Here, the coefficient $c_{\mathrm{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. $\left|\boldsymbol{M}\right|=1$.

The continuum magnetisation field $\boldsymbol{M}$ is equal to the normalised ensemble average of the volume average of the atomic spin magnetic moments $\boldsymbol{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 [24].

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.

2.2.2. A continuum model based on upscaling

The basic idea behind upscaling approaches is to start by assuming a macro model, in which certain quantities are unknown and must be obtained from a given microscopic model. The form of the macro model usually requires some knowledge about the physical laws that govern the evolution of macroscopic variables. A macro model in the form of

[TABLE]

has been proposed and analysed in [7]. In this macro model, the term $\boldsymbol{F}$ is an unknown quantity, which is then upscaled using the local microscopic equation. On the other hand, the macroscopic response, $\boldsymbol{M}$, is obtained by using a suitable time discretisation. While designing such upscaling strategies, one important issue is the synchronisation of the micro problems using the coarse-scale variables, which is achieved by assigning suitable initial and boundary conditions for the micro problems. Note that the long-range field does not appear explicitly in the macro model (2.13), but the macro model 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.

2.3. Quantification of the errors associated with approximations of the demagnetising field

It is well-known that from the purely mathematical point of view, 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_{\mathrm{c}}$, which gives

[TABLE]

The atomistic expression (2.6) can be seen as a discretisation of the integral expression for $\boldsymbol{H}_{\mathrm{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 $\boldsymbol{H}_{\mathrm{c}}$ and the atomistic field $\boldsymbol{H}_{\mathrm{a},i}$, and b) the geometric error, which is due to neglecting the far-field particles in the computation of $\boldsymbol{H}_{\mathrm{a},i}$. For the sake of comparison, a brief derivation of a 1D solution for the magnetic potential is provided here.

2.3.1. Analytical solution for the magnetic potential in 1D

Consider the domain $\Omega_{\mathrm{in}}=[0,R]\times\mathbb{R}^{2}$, which is bounded in the $x_{1}$ direction and infinite in $x_{2}$ and $x_{3}$ directions. Assume also a magnetisation vector field $\boldsymbol{M}=\boldsymbol{M}(x_{1})$, which is a function of the first coordinate only. Then it follows that

[TABLE]

This is a one-dimensional problem, i.e. $U_{\mathrm{in}}$ and $U_{\mathrm{out}}$ depend only on $x_{1}$, since the magnetisation $\boldsymbol{M}$ depends only on $x_{1}$. Hence the Laplace operator $\Delta$ reduces to the second derivative in $x_{1}$, and the following solution is directly obtained:

[TABLE]

For $U_{\mathrm{out}}$ to be bounded as $x_{1}\to\infty$, it is necessary that $\partial_{x_{1}}U_{\mathrm{out}}(x_{1})=b_{1}=0$. This implies that $U_{\mathrm{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

[TABLE]

Using these equalities in the equation (2.15) yields $c_{1}=0$. Hence, $U_{\mathrm{in}}$ becomes

[TABLE]

This leads to

[TABLE]

It must be noted that in the derivation above, $U_{\mathrm{out}}$ does not necessarily decay to zero, but remains bounded for all $x_{1}$. However, the derivative $\partial_{x_{1}}U_{\mathrm{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 truly three dimensional problem, the solution itself also has to decay at infinity.

Analytical solutions are also available in a few specific three-dimensional domains with a uniform magnetisation over the domain [2].

2.3.2. A qualitative comparison between $\boldsymbol{H}_{\mathrm{c}}$ and $\boldsymbol{H}_{\mathrm{a},i}$

In what follows, the aim is to demonstrate that the terms $\boldsymbol{H}_{\mathrm{a},i}$ and $\boldsymbol{H}_{\mathrm{c}}$, which are given by (2.6) and (2.10), respectively, match qualitatively. The domain $\Omega_{\mathrm{in}}=[0,1]\times\mathbb{R}^{2}$ is taken and magnetisation $\boldsymbol{M}=\boldsymbol{M}(x_{1})$ is assumed to be dependent on the first coordinate only. From (2.17), the continuum solution is given by $\boldsymbol{H}_{\mathrm{c}}(\boldsymbol{x})=-c_{\mathrm{L}}M_{1}(x_{1})\boldsymbol{e}_{1}$. For the purpose of this example, $c_{\mathrm{L}}=1$ is assumed. To compute the atomistic field $\boldsymbol{H}_{\mathrm{a},i}$ via equation (2.6), first, a truncation of the domain $\Omega_{\mathrm{in}}$ is required. Therefore, the following hyper-rectangle is defined:

[TABLE]

Next, the atomistic lattice is defined as a uniform discretisation of $\Omega_{R}$, where the step size of the discretisation is chosen to be equal to the interatomic distance $a$. Namely,

[TABLE]

where

[TABLE]

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

[TABLE]

Figure 2.1 confirms a qualitative match between the atomistic expression (2.6) for $\boldsymbol{H}_{\mathrm{a,i}}$ and the continuum equation (2.10) for $\boldsymbol{H}_{\mathrm{c}}$. It is evident, from equation (2.6), that the computed atomistic field $\boldsymbol{H}_{\mathrm{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 Figure 2.1. The next subsection focuses on studying 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.

2.3.3. The approximation error

Given a magnetic body $\Omega\subset\mathbb{R}^{3}$ filled with a number of atoms, the continuum field $\boldsymbol{H}_{\mathrm{c}}(\boldsymbol{x}_{i})$, $\boldsymbol{x}_{i}\in\Omega$ is obtained as the limit of the atomistic field $\boldsymbol{H}_{\mathrm{a},i}$ when $a\to 0$, i.e.

[TABLE]

The aim here is to analyse the rate of convergence by considering the volume $[0,1]^{2}\times[-1,1]$ in $\mathbb{R}^{3}$ and computing the atomistic field $\boldsymbol{H}_{\mathrm{a},i}$ for decreasing values of $a$. The magnetic moments are assumed to be uniform everywhere and pointing in the $x_{1}$ direction,

[TABLE]

In particular, a sequence of atomic distances given by $a_{k}=2^{k}a_{\mathrm{min}}$, where $k=0,1,\ldots,5$ and $a_{\mathrm{min}}=0.00625$ is used, and the following differences are recorded.

[TABLE]

The convergence rate can then be obtained by computing $\{s_{k}\}_{k=1}^{4}$, where

[TABLE]

The values of $E_{k}$ and $s_{k}$ are shown in Table 2.1, and a first order convergence rate with respect to the interatomic spacing for the approximation error is observed for this specific example.

2.3.4. The geometric error

When computing the atomistic field $\boldsymbol{H}_{\mathrm{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 $\boldsymbol{H}_{\mathrm{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}\times[-R,R]$ is uniformly discretised using the interatomic distance $a=0.1$. For the computations, it is assumed that

[TABLE]

The convergence is studied for $R_{k}=2^{k},k=0,1,\ldots,4$, and the errors

[TABLE]

are recorded. To find the rate $s$ of the geometric error, which is assumed to be $O(R^{-s})$, the following ratios are computed.

[TABLE]

The results summarised in Table 2.2, show a second order convergence rate for the geometric error.

3. Multiscale modelling

3.1. Partitioned-domain approach

3.1.1. 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 between these two regions, as illustrated in figure 3.1.

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 [25]. 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 [24]; however, it has a more complex implementation than FDM. In the case of FDM-discretised continuum, the coupling requires boundary conditions for the continuum region, which can be obtained from additionally constructed pad 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 atomistic-continuum interface [22, 24]. 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 [25].

3.1.2. 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 [24], 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 long-range dipole-dipole interactions in multiscale models of magnetic materials is to model these interactions using an exclusively continuum approach. To the best knowledge of the authors, this idea has not been suggested so far in the context of the multiscale approach in application to the magnetisation dynamics. 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 figure 3.1. 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 $\boldsymbol{H}_{\mathrm{a},i}$ is equated to the continuum demagnetising field $\boldsymbol{H}_{\mathrm{c}}$, while the magnetisation $\boldsymbol{M}$ in (2.12) is, in turn, equated to the atomistic solution $\boldsymbol{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\to 0$.

There will be differences, however, depending on whether the approach is energy-based or force-based and whether FEM or FDM is used for the continuum region. In the energy-based approach and FEM-discretised continuum with the finite-element mesh refined down to the atomistic lattice, figure 3.1a, 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 a system will be energy-conserving.

In the force-based approach and FDM-discretised continuum with the structured mesh, figure 3.1b, the auxiliary mesh will be also structured and will 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 $\boldsymbol{M}$ in (2.12) is obtained by volume-averaging of the atomistic spin magnetic moments, while the atomistic demagnetising field $\boldsymbol{H}_{\mathrm{a},i}$ is obtained by interpolation of the continuum demagnetising field $\boldsymbol{H}_{\mathrm{c}}$.

It is also possible to have the force-based approach and FEM-discretised continuum. In this case, since the conservation of the total energy becomes irrelevant, the interpolation can be used, which means that mesh used for the solution of the equation for the demagnetising field can be arbitrary. The final combination of the energy-based approach and 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.

3.1.3. The scheme used in the numerical examples

In the numerical examples of this paper, the continuum region is discretised using the 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, padding atoms are constructed, see figure 3.1. The solution at padding atoms is obtained by bilinear interpolation of the continuum solution with subsequent normalisation. The normalisation is introduced to preserve the length of spin magnetic moments.

To be able 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 $\Delta x$ centred at the the node, where $\Delta 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 $\Delta 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, the magnetisation $\boldsymbol{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 atomistic-continuum interface [23, 24]. This damping acts as a low-pass filter for the waves travelling from the atomistic region to the continuum, 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 [23].

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 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.

3.2. 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 equation (2.1) also includes a white noise term. The noise term results in fluctuations of the magnetic moment vectors $\boldsymbol{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 equation (2.1) with $\boldsymbol{m}_{i}$, it is easy to see that the length of each individual moment $\boldsymbol{m}_{i}(t)$ is equal to $\left|\boldsymbol{m}_{i}(t)\right|=1$, $\forall t$. However, due to the thermal fluctuations, the statistical averages acquire reduced lengths, i.e. $\left|E[\boldsymbol{m}_{i}](t)\right|<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 subsections 3.2.2 and 3.2.3 of this paper, an extension of the algorithms from [7, 6] is presented, after introducing the mathematical tools and notations in subsection 3.2.1. In particular, it is demonstrated that the micro problems 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.

3.2.1. Averaging kernels

In this subsection, the basic mathematical tools and notations for the HMM algorithms in subsections 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 to the size of a few interatomic distances, i.e. $\eta=ma$, $m\in\mathbb{Z}^{+}$. The temporal averaging, however, takes place on a domain of size $\tau=O(\varepsilon)$, where $\varepsilon$ is a time scale, at which the microscopic dynamics undergoes some variations. For averaging, the space $\mathbb{K}^{p,q}$ of averaging kernels (weight functions) is introduced. The space $\mathbb{K}^{p,q}$ consists of functions $K$, which have compact support in $[-1/2,1/2]$, and

• •

$K$ is symmetric, i.e., $K(t)=K(-t)$

• •

$K^{(q+1)}(t)\in BV(\mathbb{R})$, where $BV$ is the space of functions with bounded variations in $\mathbb{R}$

• •

$K$ has $p$ vanishing moments, i.e.

[TABLE]

Applying a kernel $K\in\mathbb{K}^{p,q}$ to an $\varepsilon$-periodic function $f^{\varepsilon}(t)=f(t/\varepsilon)$, where $f$ is $1$-periodic, with an average defined as $\bar{f}:=\int_{0}^{1}f(s)\;\mathrm{d}s$ results in the following arbitrarily high convergence rates, see e.g. [5, 4],

[TABLE]

where $K_{\tau}(\cdot):=\frac{1}{\tau}K(\cdot/\tau)$ is a scaled kernel, $q$ is the smoothness parameter associated with $K$, and

[TABLE]

Note that a constant kernel belongs to the space $\mathbb{K}^{1,-1}$, i.e. $q=-1$, and, therefore, the corresponding error becomes $O(\varepsilon/\tau)$, in view of the estimate (3.1). In general, smoother the kernel $K$ (higher $q$), higher the convergence rates become. Moreover, if $\tau$ is an exact integer multiple of $\varepsilon$, 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 [5, 4, 12].

These local averaging kernels will be used in the description of the HMM algorithms below.

3.2.2. 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+\ell)L+1$ number of magnetic moments that are located on a set of discrete points in 1D, i.e. $\{x_{i}=ia\}_{i=0}^{(r+\ell)L}$, where $a$ represents the interatomic distance, $r\in\mathbb{Z}^{+}$ and $\ell\in\mathbb{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 (2.1) with $\boldsymbol{H}_{i}=\left(\sum_{j}J_{ij}\boldsymbol{m}_{j}+\boldsymbol{H}^{\varepsilon}_{\mathrm{e},i}(t)+\boldsymbol{H}_{\mathrm{a},i}\right)$, where $\boldsymbol{H}^{\varepsilon}_{\mathrm{e},i}(t)$ is a high-frequency external field oscillating with the wavelength $\varepsilon$ in time. The index $i$ in $\boldsymbol{H}_{\mathrm{e},i}^{\varepsilon}$ is to allow for spatially non-uniform external fields. For the exchange coefficient a nearest-neighbour interaction is assumed, i.e.

[TABLE]

The macroscopic variable is defined as the local average of $(2r+1)$ microscopic magnetic moments. To define the macro variable, assume that the coarse grid is given by $\{X_{I}=I(r+\ell)a\}_{I=0}^{L}$ with $L\ll N$, implying much fewer degrees of freedom in comparison to a full atomistic model. The macroscopic magnetisation $\boldsymbol{M}_{I}$ at a point $X_{I}$ is defined as

[TABLE]

where $\eta=(2r+1)a$ is the size of the local spatial averaging domain and $\tau>\varepsilon$ is that of a temporal averaging. From formula (3.2), it is evident that between two consecutive macroscopic points, a total number of $\ell$ magnetic moments is skipped while averaging.

Macro model. The macro model takes the form

[TABLE]

where $\tilde{I}=\{I-1,I,I+1\}$.

Micro model. To close the macro problem, $\boldsymbol{F}_{I}(t^{*},\boldsymbol{M}_{\tilde{I}})$ must be computed. To do this, a set of coupled non-linear ODEs for $\boldsymbol{m}_{Ir^{\prime}+j}(t)$, where $r^{\prime}=r+\ell$, $t\in\mathcal{I}^{\pm}_{\tau}$, $\mathcal{I}^{+}_{\tau}=t^{*}+[0,\tau/2]$, $\mathcal{I}^{-}_{\tau}=t^{*}+[-\tau/2,0]$, with $\tau/2>\varepsilon$, and $j=-r+1,\ldots,r-1$ is solved, i.e.

[TABLE]

where $\hat{\boldsymbol{M}}=\pi_{2}\boldsymbol{M}/|\pi_{2}\boldsymbol{M}|$ is obtained by a normalised second order polynomial interpolation of the macroscopic solutions $\boldsymbol{M}_{\tilde{I}}$ for $\tilde{I}=I-1,I,I+1$.

Upscaling. The last step is to upscale the quantity $\boldsymbol{F}_{I}(t^{*},\boldsymbol{M}_{\tilde{I}})$ in (3.3) by

[TABLE]

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 micro problem (3.4). Moreover, the boundary atoms and the initial data of the micro problem are forced to be equal to the coarse-scale variables, to synchronise the microscopic model with macro variables. In (3.4), the quantity $\boldsymbol{H}_{\mathrm{c}}$ is computed by solving the equation (2.12) on the macroscopic grid. Moreover, in the computation of $\boldsymbol{H}_{\mathrm{c}}$, the macro solutions $\boldsymbol{M}$ are used in the right hand side of (2.12). It is worth mentioning that the damping term is not included in the micro model, since it is modelled at the macroscopic level, 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 $\boldsymbol{H}_{\mathrm{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 [15] 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 $\Delta t$, which is much larger than a time step $\delta t$ used for a discretisation of the micro model (3.4).

3.2.3. 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

[TABLE]

be adopted. Note the differences between $\boldsymbol{H}_{\mathrm{det},i}^{\varepsilon}$ and $\boldsymbol{H}_{\mathrm{sto},i}$. The term $\boldsymbol{H}_{\mathrm{det},i}^{\varepsilon}$ is deterministic but oscillatory, while $\boldsymbol{H}_{\mathrm{sto},i}$ is stochastic and includes the filtered external field $\left(K_{\tau}\ast\boldsymbol{H}_{\mathrm{e},i}^{\varepsilon}\right)(t^{\ast})$. The superscript $\varepsilon$ in the term $\boldsymbol{H}^{\varepsilon}_{\mathrm{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.

Macro model. With a slight deviation to the algorithm at zero temperature, the macro model at nonzero temperatures is given by

[TABLE]

where $\tilde{I}:=\left\{I-1,I,I+1\right\}$, $\boldsymbol{F}_{I}$ is the missing data in the model and $s_{I}(t)\boldsymbol{M}_{I}(t)$ is the ultimate macro solution, which models the coarse-scale dynamics. In particular, $\boldsymbol{M}_{I}(t)$ has unit length, up to an upscaling error, and represents the direction and $s_{I}(t)\leq 1$, which is computed below, accounts for the reduction in the magnetisation length.

Micro model. To compute $\boldsymbol{F}_{I}(t^{*},s_{\tilde{I}}\boldsymbol{M}_{\tilde{I}})$, first the micro problem

[TABLE]

needs to be solved. Here $t\in\mathcal{I}^{\pm}_{\tau}$, $\mathcal{I}^{+}_{\tau}:=t^{*}+[0,\tau/2]$, $\mathcal{I}^{-}_{\tau}:=t^{*}+[-\tau/2,0]$, and $r^{\prime}=r+\ell$. The index $j$ is in the range $j=-r+1,\ldots,r-1$. Moreover, the final microscopic time $\tau/2$ satisfies $\tau/2>\varepsilon$ and $\hat{\boldsymbol{M}}(x)=\pi_{2}(s_{\tilde{I}}\boldsymbol{M}_{\tilde{I}})/|\pi_{2}(s_{\tilde{I}}\boldsymbol{M}_{\tilde{I}})|(x)$ denotes the normalised second order polynomial interpolation of the macroscopic solutions. Finally, the term $\boldsymbol{H}_{\mathrm{c}}$ is computed similarly as in the zero temperature HMM algorithm in the previous subsection.

Upscaling. The quantity $\boldsymbol{F}_{I}(t^{*},s_{\tilde{I}}\boldsymbol{M}_{\tilde{I}})$ is computed by

[TABLE]

where $\eta=(2r+1)a$.

Computation of $s_{I}(t)$. In the final step, the quantity $s_{I}(t^{*})$ is computed by solving the following stochastic LLG equation for $j=-r+1,\ldots,r-1$ and $t\in[t^{*},t^{*}+\tau_{\mathrm{f}}]$

[TABLE]

where $\tau_{\mathrm{f}}>\tau_{\mathrm{r}}$, and $\tau_{\mathrm{r}}$ is the time it takes to reach the thermal equilibrium, and $\hat{\boldsymbol{M}}$ is defined similarly as in the micro problem. Then, with $\eta=(2r+1)a$, the following is computed

[TABLE]

Note that in the zero and the non-zero temperature algorithms, the long-range continuum field appears in the micro problem. Moreover, for the non-zero temperature HMM, it is necessary to include the long-range continuum field in the length scaling procedure, equation (3.2.3), 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.

4. Computational examples

4.1. 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 [24], such a problem was considered and the performance and the advantages of the partitioned-domain multiscale technique were demonstrated. However, in [24], 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 micro-crack of the sample or to impurity atoms.

The domain wall in a material with the 2D $\left(111\right)$ 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.

4.1.1. Computational setup and model parameters

Since all quantities are considered to be dimensionless, $\beta_{\mathrm{L}}=1$ was used. The damping was selected to be $\alpha_{\mathrm{L}}=0.1$. Atoms were selected to be arranged according to 2D $\left(111\right)$ fcc stacking. Lattice spacing was taken to be $a=1/64$. The exchange coefficients were taken to be $J_{ij}/\mu=(2/3)a^{-2}$, which gives the continuum exchange tensor $\boldsymbol{A}_{\mathrm{e}}=A_{\mathrm{e}}\boldsymbol{I}$ with $A_{\mathrm{e}}/\mu=1$, where $\boldsymbol{I}$ is the identity tensor. The anisotropy term was not taken into account either, $\boldsymbol{K}_{\mathrm{a}}=0$. The coefficient that defines the magnitude of the demagnetisation field was taken to be $c_{\mathrm{L}}=2\pi^{2}$, which gives the approximate width of the domain wall $w_{\mathrm{D}}=1$. The external field was applied in the $x$-direction, $\boldsymbol{H}_{\mathrm{e}}=H_{x}\boldsymbol{e}_{x}$, where $H_{x}=5$ was assumed.

The width and the height of the computational region were taken to be $x_{\mathrm{L}}=4$ and $y_{\mathrm{L}}=1$, respectively. Neumann boundary conditions were used at $x=0$, $x=x_{\mathrm{L}}$, $y=0$, and $y=y_{\mathrm{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 $\Delta x=4a$. The atomistic subregion was located in the centre of the computational region. The width and the hight of the atomistic subregion were taken to be $x_{\mathrm{A}}=y_{\mathrm{A}}=48a=0.75$. The atomistic subregion contained a hexagonal void with a side of $na$, where $n\in\left\{3,4,5,6\right\}$, the centre of which was located at $x=2$ and $y=0.5$. Time step $\Delta 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 [24], 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_{\mathrm{D}}=625$, $s_{\mathrm{A}}=3\Delta 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:

[TABLE]

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$.

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, $\boldsymbol{M}$, is defined at the grid points, an auxiliary quantity $\zeta=\operatorname{arccos}\left(M_{z}\right)-\pi/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 $\zeta=0$.

4.1.2. Results

In figure 4.1, 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 of the domain wall accelerates and overtakes the lower part ($t=1.4$ and $t=1.6$). 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 to 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 $\boldsymbol{e}_{z}$ or opposite to $\boldsymbol{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 is 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 figure 4.1 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 figure 4.2. 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 $\boldsymbol{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, which is vertical. This corresponds to the moment when the upper part of the domain wall moves faster than the lower part.

To understand the dependence of the kinetics of the domain wall on the size of the void, the position of the domain wall centre was calculated for each time step. The dependencies of the positions of the domain wall centre at $y=0.25$ and at $y=0.75$ on time are shown in figure 4.3. The void is located at $x=2$ and it is clearly seen that when the domain wall passes the void, the velocity drops. The decrease of the velocity non-linearly depends on the size of the void. For the void of size $3a$ the decrease of the velocity is relatively small, while for $5a$ it is already significant. For the void of the size of $6a$, the position of the domain wall stays below $2.5$, i.e. the domain wall does not move past the void and stays at $x\approx 2$.

In figure 4.3, it is also seen that there is a difference between the positions at $y=0.25$ and at $y=0.75$. Increases and subsequent decreases of the positions, where the lower part of the domain wall overtakes and then falls behind the upper part of the domain wall, are the oscillations created by the void. The magnitude of these oscillations also depends on the size of the void and increases with the increase of the size (in the considered example, the highest oscillations in the domain wall position are observed for the void of size $6a$).

4.2. 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_{\mathrm{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, $\boldsymbol{H}_{\mathrm{c}}=-c_{\mathrm{L}}M_{x}\boldsymbol{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), $C\boldsymbol{M}$ can be added to effective field $\boldsymbol{H}$, where $C$ is an arbitrary constant, without influencing the solution of the LLG equation. This is used in the example below.

4.2.1. Computational setup and model parameters

The same $\beta_{\mathrm{L}}$ and $\alpha_{\mathrm{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}/\mu=a^{-2}$, which gives the continuum exchange parameter $A_{\mathrm{e}}/\mu=1$. Biaxial anisotropy was used,

[TABLE]

The introduction of parameter $c_{\mathrm{L}}$ into the anisotropy is explained below. The coefficient that defines the magnitude of the demagnetisation field, $c_{\mathrm{L}}$ was varied from $0.4$ to $12.8$. These parameters give the approximate width of the domain wall $w_{\mathrm{D}}=1$. The external field was applied in $z$-direction, $\boldsymbol{H}_{\mathrm{e}}=H_{z}\boldsymbol{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_{\mathrm{L}}$, because

[TABLE]

and addition of $c_{\mathrm{L}}\boldsymbol{M}$ to $\boldsymbol{H}$ does not affect the solution of the LLG equation. 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 $\boldsymbol{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_{\mathrm{L}}$.

The length of the computational region was taken to be $x_{\mathrm{L}}=6$. Neumann boundary conditions were used at $x=0$, $x=x_{\mathrm{L}}$ for the LLG equation. The 1D analytical solution of the magnetostatic equation, $\boldsymbol{H}_{\mathrm{c}}=-c_{\mathrm{L}}M_{x}\boldsymbol{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 $\Delta 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_{\mathrm{A}}=128a$. Time step $\Delta t=10^{-2}$ was used. Same $g_{\mathrm{D}}$, $s_{\mathrm{A}}$ and the width of the damping band as in example above are used.

The following initial conditions for the domain wall were used:

[TABLE]

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_{\mathrm{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_{\mathrm{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.

4.2.2. Results

In figure 4.4, the dependence of the multiscale modelling error on $c_{\mathrm{L}}$ is shown. This is the error due to numerical multiscale treatment of the demagnetisation field, compared to the reference case, where the analytical expression for the demagnetisation field is used. The error is proportional to $c_{\mathrm{L}}$.

4.3. HMM examples

4.3.1. 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 lattice consisting of points $\{ia\}_{i=0}^{100}$, where $a=0.01$ is the atomic distance and the initial configuration is given by

[TABLE]

Moreover, the high frequency external field

[TABLE]

with $f^{\varepsilon}(t)=1+\cos(0.43t)+\cos^{2}(2\pi t/\varepsilon)$ and $\varepsilon=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 $\boldsymbol{H}_{\mathrm{c}}$ in the microscopic model (3.4), is given exactly by

[TABLE]

in this one-dimensional setting, cf. the derivations in Section 2.3.1. The parameters $\beta_{\mathrm{L}}=\alpha_{\mathrm{L}}=\mu=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 $\Delta t=0.01$ on the macro-scale and $\delta t=\varepsilon/10$ on the micro-scales. Moreover, the temporal microscopic box is $\tau=5\varepsilon$. 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 4.5 demonstrates the evolution of the macroscopic dynamics, and compares it to the solution of the full atomistic model. In the atomistic simulation it is assumed that $\boldsymbol{H}_{\mathrm{a},i}=-m_{x}\boldsymbol{e}_{x}$, which is justified by the convergence of the atomistic long-range field to the continuum counterpart as $a\to 0$, see Section 2.3.3. As time increases, the magnetisation vectors are pointing to the $\boldsymbol{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.

4.3.2. 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_{\mathrm{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_{\mathrm{L}})$, see Figure 4.6.

4.3.3. A chain of magnetic particles at nonzero temperature

In this section, the aim is to show an example, where the existence of the long-range interaction has an effect over the macroscopic magnetisation length at elevated temperatures. For a numerical study, $N=101$ atomistic particles located on a one-dimensional lattice consisting of points $\{ia\}_{i=0}^{100}$, where $a=0.01$ is the atomic distance, are considered. The initial configuration is uniform and given by

[TABLE]

equipped with periodic boundary conditions. The anisotropy is assumed to be zero, and the external field is pointing in the $x$ direction, i.e.

[TABLE]

The atomistic particles are subjected to a thermal noise with a standard deviation of $D=0.01$, and the parameters $\beta_{\mathrm{L}}=\mu=1$, $\alpha_{\mathrm{L}}=10$, with the nearest neighbour interactions for the exchange coefficient, similar to previous examples, are chosen for the simulation. In Figure 4.7, 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 Figure 4.8, 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, $|\langle\boldsymbol{m}\rangle|$, based on the closure argument $E[\boldsymbol{m}_{i}\cdot\boldsymbol{m}_{j}]=E[\boldsymbol{m}_{i}]\cdot E[\boldsymbol{m}_{j}]$, and the restrictive assumption that $E[\boldsymbol{m}_{i}]=E[\boldsymbol{m}_{j}]$ for all $i,j$. The mean-field formula, for this specific setup, reads as

[TABLE]

see e.g. [6, 2] 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 long-range 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 atomic moments are correlated through the short-range exchange and the long-range dipole-dipole interactions.

5. Conclusions

This paper demonstrates a way of including the long-range dipole-dipole interactions between the atomistic spin magnetic moments into the 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 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 in both partitioned-domain and upscaling strategies, the modelling error is $O(c_{\mathrm{L}})$, where $c_{\mathrm{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 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 in the energy-based methods and 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 the continuum mesh within 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 void-affected kinetics of the domain wall was modelled using the partitioned-domain approach. Only a small region around the void was modelled at the atomistic scale, while the rest of the 2D structure was modelled using the continuum model with a coarser resolution. The domain wall structure itself was 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.

As for the upscaling strategy, the main novelty and advantage with the proposed algorithm is that the long-range atomistic communication between the magnetic moments, located in different microscopic boxes, is avoided. Nevertheless, the macroscopic effect of the long-range interactions is captured accurately in a multiscale formalism. This leads to a significant computational gain not only in comparison to a naive computation of the dipole-dipole interactions (which scales quadratically with respect to the number of particles) but also when compared to more efficient multiscale methods such as the fast multipole method [19, 17], which is a linear scaling algorithm. The fact that the long-range field is included in the microscopic model, using the continuum field $\boldsymbol{H}_{\mathrm{c}}$ allows for obtaining sublinear 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.

To sum up, an efficient multiscale modelling of the demagnetisation field for two well-established multiscale formalisms, the partitioned domain and the upscaling approaches, has been proposed. The ideas presented in this paper allow for a complete multiscale simulation of ferromagnetic materials, as they can be used when designing efficient modelling strategies taking into account both short- and long-range interactions between the spin magnetic moments and a finite temperature.