Coupling atomistic and continuum modelling of magnetism

Accepted Manuscript Coupling atomistic and continuum modelling of magnetism M. Poluektov, O. Eriksson, G. Kreiss

PII: DOI: Reference:

S0045-7825(17)30246-3 https://doi.org/10.1016/j.cma.2017.10.010 CMA 11637

To appear in:

Comput. Methods Appl. Mech. Engrg.

Received date : 5 February 2017 Revised date : 24 August 2017 Accepted date : 5 October 2017 Please cite this article as: M. Poluektov, O. Eriksson, G. Kreiss, Coupling atomistic and continuum modelling of magnetism, Comput. Methods Appl. Mech. Engrg. (2017), https://doi.org/10.1016/j.cma.2017.10.010 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Coupling atomistic and continuum modelling of magnetism M. Poluektova,∗, O. Erikssonb , G. Kreissa a Department

of Information Technology, Uppsala University, Box 337, SE-751 05 Uppsala, Sweden b Department of Physics and Astronomy, Uppsala University, Box 516, SE-751 05 Uppsala, Sweden

Abstract In this article, a new energy-based atomistic-continuum partitioned-domain multiscale method for magnetisation dynamics is proposed. The main feature of the method is the minimised mismatch between the local continuum description and the non-local atomistic description (in which non-nearest neighbour interatomic interactions are present). The error at the atomistic-continuum interface is minimised by constructing an intermediate region consisting of transition atoms, which interact differently with the neighbouring atomistic and continuum regions. When the mesh of the continuum region is refined to the atomistic scale at the atomistic-continuum interface and an energy-conserving time-stepping scheme is used for the entire computational domain, the method conserves the total energy of the system. The second feature of the method is the introduction of a damping band at the atomistic-continuum interface, which absorbs high-frequency waves that are otherwise reflected from the interface and thereby contribute to the error inside the atomistic region. Several examples of domain wall motion and spin wave propagation in one- and two-dimensional structures are used to illustrate the applicability of the method and to investigate its limitations. Keywords: magnetisation dynamics, multiscale modelling, ∗ Corresponding

author Email addresses: [email protected] (M. Poluektov), [email protected] (O. Eriksson), [email protected] (G. Kreiss)

Preprint submitted to Computer Methods in Applied Mechanics and EngineeringOctober 12, 2017

atomistic-continuum coupling, domain walls, ghost forces

1. Introduction Physical properties of materials are often highly dependent on the underlying microstructure and functionalities are often affected by the presence of defects and/or heterogeneities at the atomistic length scale. Examples of this 5

are the pinning of domain walls to impurity atoms, something which gives rise to the so-called Barkhausen effect [1], and magnetic nano-structures for enhanced coercive fields and remanence [2]. Multiscale atomistic-continuum models allow investigating the influence of atomistic-scale phenomena on macroscopic quantities without using a fine-scale resolution over the entire computational

10

domain. Over the past decades, several ways of modelling the dynamic behaviour of magnetisation in magnetic materials were developed. There are theories which describe magnetisation dynamics of atomistic systems [3, 4, 5, 6, 7], where magnetic moments of individual atoms are considered. Although being compu-

15

tationally expensive at relatively large scales (106 –109 atoms), this approach provides a reliable description of the material behaviour and, for example, was used to describe ultrafast switching in a magnetic random-access memory device [8]. On length scales of micrometres, magnetisation dynamics is usually described using micromagnetics [9, 10], where the magnetisation is considered

20

to be a continuum vector field, which is described by the set of specific partial differential equations (PDEs). However, this approach relies heavily on the assumption of smooth magnetisation distribution and cannot accurately describe regions around defects, such as nano-voids, nano-cracks or grain boundaries. Thus, when continuum models are applied at the macroscopic scale, they con-

25

tain effective parameters, in which the microstructural information is indirectly embedded. A more elaborate discussion of applicability of different modelling paradigms can be found in [7] and is out of the scope of the current work. There are multiple ways of performing multiscale simulations, one of which

2

is domain partitioning, where an explicit interface between various material de30

scriptions is implied. Domain partitioning is a form of two-way coupling between models/scales, as opposed to one-way coupling, in which the effective properties are obtained at the microscopic scale and are subsequently used at the macroscopic scale. Methods that use an atomistic-continuum domain partitioning were initially developed to describe the mechanical behaviour of crystalline

35

solids. Since atomistic-continuum techniques for mechanical and magnetic descriptions have many similarities, a brief summary of coupling methods for mechanics is given below. There is a wide variety of methods, which are applicable in the case of quasistatic behaviour. Only few examples are listed here and for the exhaustive overview the reader is referred to the excellent review article [11].

40

The development of partitioned-domain methods in application to mechanical systems started from the FEAt method [12]. At the moment, one of the most well-known methods is the quasicontinuum (QC) method [13, 14]. Another notable example is the CADD method [15] for modelling plasticity of metals in the quasistatic regime (later extended to dynamics). In this method, the movement

45

of dislocations between the atomistic and the continuum regions was included. Most partitioned-domain methods can be separated into two categories: energy-based (such as the classical QC method) and force-based (such as FEAt and CADD) methods [16]. The energy-based class of methods relies on the ability to construct the total energy of the system, from which the forces are

50

derived. However, due to the mismatch between the non-local atomistic interactions and the local continuum description, so-called ghost-force error (or simply ghost forces) often appear at the interface between the regions, which leads to a distorted configuration. The force-based methods are designed such that perfect crystals have a correct equilibrium state across the atomistic-continuum

55

interface; however, these methods lack a well-defined total energy functional, which can be problematic when material dynamics is modelled. To deal with the ghost forces, the concept of quasi-non-local (QNL) atoms, which interact with a different number of neighbouring atoms depending on direction, was introduced in [17]. However, this approach was initially appli3

60

cable only to certain crystal types, a flat atomistic-continuum interface and a second nearest neighbour interatomic interaction. In [18, 19, 20], the QNL-type coupling was significantly advanced. In [20], the generalisation of this coupling, denoted as GRAC, which completely eliminates the ghost forces under certain conditions, was presented. There are also other efficient ghost force reduction

65

techniques, for example [21], and the exhaustive review can be found in [22]. In the case of dynamic behaviour, the goal of modelling can be capturing the equilibrium properties of a system at a finite temperature. An example of the multiscale method for such problem is the finite-temperature QC method [23]. However, sometimes non-equilibrium (fast) processes should be modelled. In

70

this case, non-equilibrium disturbances should propagate within the atomistic region, but should not reflect from the atomistic-continuum interface back to the atomistic region and thereby introduce a numerical error to the solution. One of the non-equilibrium multiscale approaches is the finite-temperature CADD method [24, 25], where a fully dynamic atomistic region was coupled to a qua-

75

sistatic continuum region. In this method, the problem of wave reflections was addressed by applying a thermostat to the atoms near the atomistic-continuum interface. There are several initial attempts to construct multiscale techniques for simulating the magnetic properties of materials. In [26], an atomistic-continuum

80

multiscale approach was used to model a domain wall in the FePt/FeRh bilayers. Within the suggested framework, the continuum discretisation was based on the regular grid and a sharp interface to the atomistic region was used (i.e. there was no transition zone between disparate regions). In [27, 28], some aspects of constructing continuum micromagnetic description and estimating its parameters

85

based on atomistic spin dynamics at finite temperatures were discussed. In [29], a way of coarse-graining from the atomistic to the micromagnetic description was discussed, where a gradual transition of the discretisation between the regions was considered. However, the proposed method was limited to the static problem of finding equilibrium states in a computational domain partitioned

90

into the atomistic and the continuum regions. Recently, a way of construct4

ing a transition between two descriptions for dynamic problems was published [30]. This approach suggested an overlap of the atomistic and the continuum regions, where the exchange energies of both models are mixed. However, the error at the interface was not analysed for the method. The model was used 95

to describe the dynamics of Bloch points, which were surrounded by the atomistic description, while the rest of the material was modelled as a continuum. In [31], a multiscale approach with a sharp atomistic-continuum interface was used to simulate domain walls, spin waves and magnetic skyrmions. However, the method lacked a local-non-local transition region as well as the numerical

100

damping at the interface (e.g. such as a thermostat in the CADD method). In [32], a multiscale method based on the HMM (heterogeneous multiscale methods) framework for coupling disparate time scales in the micromagnetic equations has been proposed and analysed. In this framework, the models corresponding to different scales are integrated at different temporal resolutions

105

and effective values of the required quantities are then passed from the lower to the upper scale. The numerical example of a single spin under an oscillatory external field was analysed in this study. In [33], the multiscale modelling of magnetisation dynamics where both scales were represented by the continuum description, however with different spatial discretisations, was studied. Inter-

110

facing descriptions with different spatial discretisations by domain partitioning usually leads to computational difficulties, in particular, high frequency wave reflections at the interface due to a change in discretisation. Therefore, the concept of a wave-absorbing band at the interface was introduced. The aim of this article is to introduce a new energy-based atomistic-continuum

115

multiscale method for modelling magnetisation dynamics. In contrast to other works on multiscale magnetisation dynamics, in this article, specific emphasis is made on the error due to atomistic-continuum interface and a method, which minimises this error, is proposed. To deal with the mismatch between the non-local atomistic interactions and the local continuum description, a special

120

transition region is suggested. This transition region consists of atoms with the modified interatomic interaction, which depends on the direction of the vector 5

that connects atoms. Moreover, transition atoms, which are located next to the continuum region, interact only with the nearest neighbours. This approach is a variation of the approach proposed in [17], where the concept of QNL atoms 125

was introduced. In this article, the QNL-type concept is for the first time applied to spin systems. The proposed method is similar to [19, 20], where the atomistic-continuum coupling method for mechanical problems was developed (specific similarities and differences are discussed in section 2.5.3 of this paper). The multiscale method proposed in this article is intended for dynamic simu-

130

lations, where it is beneficial to construct a methodology based on the total energy functional. In this case, if energy-conserving time-stepping method is used (and damping is absent), the total energy of the system is conserved. Being energy-based, the proposed method suffers from the ghost-force artefacts. However, it is shown that due to QNL-like transition region, the magnitude

135

of these artefacts is reduced. More precisely, it is shown that for the case of uniformly varying magnetisation (e.g. spin spiral with a uniform twist), the torque error in the multiscale system is of the same order of magnitude as the mismatch error between the continuum and the atomistic descriptions. Since the description of magnetisation dynamics does not contain forces, but rather

140

torques, the ghost-force error is referred to as the ghost-torque error (or simply ghost torques) in this article. Moreover, in this article, the concept of the damping band [33], which was previously developed for the continuum micromagnetic description, is applied to the atomistic description. In this article, the description of the method is given at first, followed by

145

a few analytical examples and a number of numerical examples. The purpose of the analytical examples is to show that in the case of a uniformly varying magnetisation, which is referred to as the exchange patch test, the ghost torques and the atomistic-continuum mismatch error are of the same order of magnitude. In the numerical examples, the material parameters are varied to find the

150

limitations of the method and, therefore, they are deliberately selected not to represent an actual material. The numerical examples are intended to provide numerical error estimates and information regarding the dependence of the er6

ror on dimensionless quantities, such as the number of atoms per domain wall. The second purpose of the numerical examples is to show that for properly se155

lected parameters of the transition region, the solution given by a multiscale simulation does not depend on the size of the atomistic region, i.e. the solution inside the atomistic region is as if this region was a part of the fully atomistic model rather than embedded into the coarser continuum description. For this purpose, examples where the effect originates from the outside of the atom-

160

istic region (a domain wall that is initiated in the continuum region and passes through the atomistic region) and from the inside of the atomistic region (an initial inhomogeneity, which creates circular spin waves) are studied. Although, only 1D and 2D numerical examples are considered in this article, the proposed method is also valid in a 3D case for arbitrary geometry of

165

the atomistic-continuum interface. This follows from the form of the equations defining the interatomic interaction of the transition atoms, in which the dimensionality of atomic position vectors is not explicitly used. There is only one requirement for the construction of the transition region — its thickness should not be smaller than the interatomic interaction radius of the real atoms.

170

2. Methods This section summarises two different modelling approaches. Both approaches are deterministic and describe materials at temperatures close to 0 K. Moreover, the long-range dipole-dipole interaction is not taken into account and an efficient multiscale way of treating this interaction type should be a subject of

175

a future research. The transition between the material descriptions in the case of sufficiently smooth magnetisation distribution is established by coarse-graining interatomic interactions and thereby determining parameters used in the continuum description, see appendix Appendix A. The presentation of the method is split into

180

several parts. In sections 2.2 and 2.3, the atomistic and the continuum models are summarised. After this, in section 2.4, the local coarse-graining continuum

7

error is analysed. Afterwards, the atomistic-continuum coupling method is presented. In sections 2.5 and 2.6, the transition region between the continuum and the atomistic regions is described. In section 2.7, an absorbing layer at the 185

interface is introduced. Finally, some aspects of the numerical time-stepping scheme are addressed in section 2.8. 2.1. Notation Vectors in R3 are denoted with an arrow, e.g. ~a. The scalar and the vector products between vectors ~a and ~b are denoted as ~a · ~b and ~a × ~b, respectively. Second-order tensors are denoted as bold capital letters, e.g. A. The tensor product between vectors ~a and ~b is denoted as ~a~b and the result of this product is a second-order tensor. The coordinate representation of the tensor product is given by Aij = ai bj , where ak and bk are components of vectors ~a and ~b, respectively, and Aij are components of tensor A. The transpose of a tensor is denoted as B = AT , which is Bij = Aji in the coordinate representation. The coordinate representations of the scalar and the vector products of a tensor by a vector from the right, d~ = A · ~c and C = A × ~c, respectively, are given by di =

3 X

Aij cj ,

Cin =

j=1

3 X 3 X

ξnjk Aij ck ,

j=1 k=1

where ξnjk is the Levi-Civita symbol. The scalar product of a tensor by a vector from the left is ~e = ~c · A = AT · ~c. The double inner product of two tensors is denoted as f = A : D, the result of which is a scalar that in the coordinate representation is given by f=

3 X 3 X

Aij Dji .

i=1 j=1

Higher-order tensors, which are formed by the tensor product of more than two vectors, are denoted with preceding superscript showing the order, e.g. 190

3

E is

a third-order tensor. For convenience, vector/tensor indices x, y, z are used instead of 1, 2, 3 in some parts of the text.

8

2.2. Atomistic spin dynamics In the atomistic approach, the dynamic behaviour of spin magnetic moments of individual atoms is described by the atomistic Landau-Lifshitz-Gilbert equation [34, 5]:   d ~ i − αL m ~i , m ~ i = −βL m ~i×H ~i× m ~i×H dt βL =

γ , µ (1 + λ2 )

αL =

|m ~ i | = 1,

γλ , µ (1 + λ2 )

(1)

(2)

where γ is the gyromagnetic ratio, λ is the phenomenological Gilbert damping constant, m ~ i is the direction of spin magnetic moment and µ is the length of spin magnetic moment1 . Equation (1) results in a damped precession of magnetic ~ i . This field is formed by the external moments around an effective field, H magnetic field and the surrounding atoms, including relatively distant atoms, and depends on the material itself. Here the following form of the effective field is considered: 

~i =  H

X j



  X ~ ij × m ~ e, Jij m ~ j +  D ~ j + Ka · m ~ i + µH

(3)

j

where Jij are constants of the Heisenberg exchange interaction2 between atoms ~ ij describe the Dzyaloshinskii-Moriya (DM) interaction, K a i and j, vectors D 195

~ e is the external field. Summations are taken for is the anisotropy tensor and H all j, such that j 6= i. In general, the interatomic exchange interaction has a range of several lattice units (5–10), i.e. it is significant not only for the nearest neighbours, although there are materials, such as magnetic oxides, with rapidly decaying exchange

200

interaction where non-nearest neighbours can be neglected. 1 In

older literature other definition of the constants in the Landau-Lifshitz-Gilbert equation

is used: βL = γ/µ and αL = γλ/µ. 2 This definition of J contains squared length of spin magnetic moments and is taken from ij [5]; however, other definitions exist [34].

9

2.3. Continuum-based description of magnetisation dynamics At the continuum scale, the magnetisation dynamics is modelled by the continuum version of the Landau-Lifshitz-Gilbert (LLG) equation [10, 35]:   ∂ ~ − αL m ~ , |m| m ~ = −βL m ~ ×H ~ × m ~ ×H ~ = 1, (4) ∂t

where m ~ is the normalised magnetisation field and βL and αL are the same as in (1). In the case of the continuum approach, the following effective field is considered: ~ = Ae : ∇∇m ~ e, H ~ + ∇ · (D e × m) ~ + Ka · m ~ + µH

(5)

where tensor Ae contains the exchange interaction constants, tensor D e contains the Dzyaloshinskii-Moriya interaction constants3 , K a is the anisotropy ~ e is the external field. These tensors are directly obtained from the tensor and H atomistic parameters: Ae =

1X Jij ~sij ~sij , 2 j

(6)

j6=i

De =

X

~ ij , ~sij D

(7)

j j6=i

where ~sij is the vector connecting atoms i and j. Since the anisotropy term is local, the same anisotropy tensor K a is used in the continuum and the atomistic equations. Equations (6) and (7) imply that Ae and D e do not depend 205

on i (are spatially homogeneous), although an extension to a spatially inhomogeneous case is straightforward4 . The derivation of relations (6) and (7) is 3 Less

general formulations of the DM interaction, in which the DM term has the form of

De ∇ × m, ~ where De is a scalar, have been used in the past. 4 Spatially inhomogeneous case corresponds to the case when continuum parameters A and e D e are dependent on the spatial coordinate ~ x. In this case, the efficiency of the continuum model should be maintained as long as these spatial changes of the parameters are gradual and, more importantly, the continuum mesh is small enough to resolve these changes. Also, differential operator in the second term of equation (5) acts only on m, ~ i.e. this form of the equation relies on D e being a constant tensor; in the spatially inhomogeneous case, the coordinate representation from appendix Appendix B has to be used.

10

provided in appendix Appendix A. Since tensorial notation is not common in micromagnetics, an alternative coordinate representation of equations (5), (6) and (7) is provided in appendix Appendix B. To simplify the notation, j 6= i is 210

omitted in all summations over j in the rest of the article. The differential operators in the continuum equations (4) and (5) can be spatially discretised in various ways, e.g. using the finite-element (FE) or the finite-difference (FD) methods. The specific form of the discretisation is not important for the problem of the atomistic-continuum coupling and both FE

215

and FD discretisations are discussed in this paper. The FE discretisation has an advantage of a possibility of creating a mesh that is refined to the atomistic lattice spacing at the interface. Therefore, in the main description of the atomistic-continuum coupling method, section 2.5, it is assumed that the continuum is discretised using the FE method. Specific changes to the scheme in

220

the case of the FD discretisation of the continuum equations are discussed in section 2.8. 2.4. Error in the continuum-based representation The continuum approach provides an approximation of the atomistic approach and hence contains an intrinsic error, which has two contributions. The

225

first contribution is the coarse-graining error, which is the error due to the change of the description (from the atomistic to the continuum). The second contribution is the discretisation error. When equations (1) and (4) are compared, it can be seen that they have the ~ i and same exact structure, while different expressions for the effective fields (H

230

~ in the atomistic and the continuum cases, respectively) are used. Therefore, H the first step towards the quantification of the error in the solution (m) ~ is the ~ This error is considered estimation of the local error in the effective field (H). in the current section and is quantified for a given solution. Although it can be ~ leads to a small error in m, expected that a small error in H ~ the analysis of the

235

solution error is non-trivial and is left out of the scope of the current work.

11

2.4.1. Coarse-graining error, exchange The continuum solution is the vector field, which is assumed to be m ~ (~x) ∈
3 C R . The atomistic solution is the discrete set of vectors m ~ i defined at 6

lattice positions ~xi . The atomistic and the continuum effective fields depend on 240

m ~ i and m, ~ respectively. Therefore, to compare the atomistic and the continuum ~ i and H, ~ respectively), a continuum magnetisation field m effective fields (H ~ that coincides with the atomistic spin magnetic moments m ~ k at lattice positions has to be considered, i.e. m ~ (~xk ) = m ~ k , where ~xk are atomic positions. This is analogous to the truncation error analysis. At first, the exchange parts of the atomistic and the continuum effective ~ EA and H ~ EC , respectively, and which are the first fields, which are denoted as H i terms of equations (3) and (5), respectively, are compared. A Taylor expansion of m ~ j in (3) is performed: 1 1 2 3 m ~j =m ~ (~xj ) = m ~ (~xi )+~sij ·∇m ~ (~xi )+ (~sij · ∇) m ~ (~xi )+ (~sij · ∇) m ~ (~xi )+. . . , 2 6 (8) ~ EA is rewritten: where ~sij = ~xj − ~xi . Using (8) and (3), the expression for H i ~ iEA H

      X X X 1 =m ~  Jij  +  Jij ~sij  · ∇m ~ +  Jij ~sij ~sij  : ∇∇m ~ + 2 j j j     X 1 1 X 3 4  Jij (~sij · ∇)  m ~ + Jij (~sij · ∇)  m ~ + . . . , (9) + 6 24 j j

where m ~ and derivatives of m ~ are taken at point ~xi . The first term in (9) can be removed, since, when (9) is substituted into (1), it is eliminated due to the vector product with m ~ i . In the case of symmetric atomic lattices, for each pair of atoms i and j, atom k exists, such that ~sik = −~sij and Jik = Jij , which gives X

Jij ~sij = ~0,

(10)

Jij ~sij ~sij ~sij = 3O,

(11)

j

X j

12

where 3O is the third-order zero tensor. This eliminates the second and the fourth terms in (9). It should be noted that for non-symmetric lattices (e.g. with defects), relations (10) and (11) are not fulfilled, which affects the error estimate. The third term in (9) corresponds to the exchange effective field in the continuum case, i.e. it is exactly the same as in (5). Therefore,   ~ EA ~ EC 1 X 4 ~ ≤ Jij (~sij · ∇) m ~ + O a6 Jˆ ∇6u m Hi − H = 24 j    X  1 4 Jij sij 4 + O a6 Jˆ ∇6u m ≤ ∇u m ~ ~ = O a4 Jˆ ∇4u m ~ , (12) 24 j p

where ∇pu = (~u · ∇) is the directional derivative of order p along the unit vector

~ takes the maximum value, ~u, the direction of which is selected such that |∇pu m|

a is the interatomic spacing, Jˆ = max |Jij | and sij = |~sij |. Estimate (12) indicates that the difference between the atomistic and the continuum exchange
effective fields is O a4 , when the atomistic solution (m ~ i ) equals the continuum solution (m) ~ at lattice points and when m ~ and Jij are fixed. It can also be useful to rewrite this error estimate in terms of components of tensor Ae , which
are O a2 according to (6). Thus, the error estimate is given by
~ EA ~ EC ~ , Hi − H = O a2 A∗ ∇4u m

(13)

A∗ = max {|Axx | , |Axy | , |Ayy |} .

245

Here Axx , Axy , Ayy are components of tensor Ae . From error estimate (13), it can also be seen that the error in the effective field is the order of the fourth derivative of the solution. 2.4.2. Coarse-graining error, DM A similar approach is applicable to the comparison of the DM parts of the ~ DA and atomistic and the continuum effective fields, which are denoted as H i ~ DC , respectively, and which are the second terms of equations (3) and (5), H respectively. A Taylor expansion of m ~ j given by (8) can be used to obtain

13

~ DA , which is the following: H i ~ iDA H

      X X X 1 ~ ij ×m+ ~ ij × (~sij · ∇) m+ ~ ij × (~sij · ∇)2  m+ = D ~  D ~ D ~ 2 j j j   1 X ~ 3 Dij × (~sij · ∇)  m ~ + . . . . (14) + 6 j

In the case of symmetric atomic lattices, for each pair of atoms i and j, atom k ~ ik = −D ~ ij , which gives exists, such that ~sik = −~sij and D X

~ ij = ~0, D

(15)

~ ij = 3O. ~sij ~sij D

(16)

j

X j

This eliminates the first and the third terms in the expansion. The second term in the expansion corresponds to the DM effective field in the continuum case. Therefore 1 X   ~ DA ~ DC ~ ij × (~sij · ∇)3 m ˆ ∇5u m ~ ≤ D ~ + O a5 D Hi − H = 6 j     X
1 ~ ˆ ∇5u m ˆ ∇3u m ≤ ∇3u m ~ ~ = O a3 D ~ = O a2 D∗ ∇3u m ~ , Dij sij 3 +O a5 D 6 j

(17)

~ ij and ˆ = max D where D

D∗ = max {|Dxx | , |Dxy | , |Dyy |} .

2.4.3. Discretisation errors The discretisation error results from the discrete representation of the differential operators. Here the 2D case of the FD discretisation is analysed for simplicity, although such error analysis can be performed for any type of dis~ EC . cretisation. The discrete form of the exchange effective field is denoted as H h

14

For the second-order FD representation, it can be shown that the discretisation error becomes the following:

250

   
∂4 ∂4 ∂4 ∂4 ~ EC ~ EC 1 2 4 + Ayy 4 m ~ +O h ≤ + H − Hh = h Axx 4 + 4Axy 12 ∂x ∂x3 ∂y ∂x∂y 3 ∂y

1 ≤ h2 A∗ ∇4u m ~ + O h4 = O h2 A∗ ∇4u m ~ , (18) 2

where h is the spatial discretisation step.

~ DC . As for the The discrete form of the DM effective field is denoted as H h exchange, for the second-order FD representation, it can be shown that  3  X 1 3
∂ ~ DC ~ DC 4 ~ ij × ∂ m~ − Hh = h2 D ~ e + m~ ~ e · ~ s + O h H x y ij ≤ 3 3 ∂x ∂y j 6   X

1 ~ ˆ ∇3 m = O h2 D∗ ∇3u m ~ ~ . ≤ h2 ∇3u m Dij sij + O h4 = O h2 aD u~ 3 j

(19)

Here D e 6= 0 was used. Obviously, the total (local) error between the discretised continuum and the atomistic descriptions is less than or equal to the sum of the coarse-graining error and the discretisation error, i.e. ~ κA ~ κC ~ κA ~ κC ~ κC ~ κC Hi − Hh ≤ Hi − H + H − Hh ,

κ ∈ {E, D}.

This can be useful for the selection of the discretisation parameters of the multiscale system, such that the total number of degrees of freedom (atoms and continuum nodes) is minimised for a given local error. 255

2.5. Atomistic-continuum transition As discussed in the introduction, there are multiple ways of constructing an interface between the continuum and the atomistic descriptions. The two key requirements for the method, which is proposed in this article, are energy conservation (if energy-conserving time-stepping method is used and damping

260

is absent, αL = 0) and small magnitude of the error in the effective field (i.e. small ghost torques). In this section, the construction of the interface between 15

the discretised continuum and the atomistic descriptions is summarised, while the motivation behind specific choices is provided in sections 2.5.1 and 2.5.2. The connection to the atomistic-continuum modelling of the mechanical be265

haviour is explored in section 2.5.3. The discussion of time-stepping and energy conservation is delegated to section 2.8. A transition region between the continuum and the atomistic regions is introduced, see figure 1. This region consists of atoms with a behaviour, which is described by a modified system of equations governing the magnetic spin

270

motion. The dynamic behaviour of these atoms is still described by equations ~ ij , are (1) and (3); however, the interatomic interaction coefficients, Jij and D modified. Although the structure of the equations is the same as for the atomistic region, atoms in this region should be understood as generic computational points rather than actual atoms. Therefore, in this article, these atoms are re-

275

ferred to as transition atoms, while within the actual atomistic region, atoms are referred to as real atoms. The transition region is split into two parts. The transition atoms that are located in the subregion, which is adjacent to the atomistic-continuum interface, interact only with their nearest neighbours and are refereed to as the fully coarse-

280

grained atoms. The rest of the transition region consists of atoms, which interact with a different number of atoms depending on the direction of ~sij , and which are referred to as the partly coarse-grained atoms 5 . For all atoms, the transition and the real atoms, the symmetry of the interaction is enforced, i.e. Jij = Jji ~ ij = −D ~ ji , see section 2.5.1 for motivation. and D

285

Within the real atomistic region, the interaction coefficients are supposed to be known, e.g. obtained from the first principles theory [36], and, for conA venience, a superscript “A” is added to them (e.g. Jij ) in further discussion.

These coefficients are used to determine the tensor parameters that are used 5 The

partly coarse-grained atoms can also be defined as atoms, which interact with larger

number of atoms than fully coarse-grained atoms, and which interact with fewer number of atoms than real atoms.

16

Figure 1: The schematic representation of the atomistic-continuum interface. The transition region is split into two parts: the fully coarse-grained atoms, which interact only with their nearest neighbours, and the partly coarse-grained atoms. The interatomic interaction ranges are indicated with the dashed lines for the different types of atoms. The interface nodes are degrees of freedom of the continuum model and they interact with the nearest transition atoms as if these nodes are the fully coarse-grained atoms. The partly coarse-grained atoms have a non-symmetric interatomic interaction, i.e. their interaction range is direction-dependent, which is schematically illustrated as a shifted ellipse.

17

in the continuum model according to equations (6) and (7). The interaction 290

coefficients of the fully coarse-grained atoms are denoted with a superscript “F” F (e.g. Jij ). The determination of these coefficients is discussed in section 2.5.1.

The determination of the interaction coefficients of the partly coarse-grained P atoms, which are denoted with a superscript “P” (e.g. Jij ) is discussed in sec-

tion 2.6. The motivation for introducing partly coarse-grained atoms (section 295

2.5.1), their relation to the local torque error (section 2.5.2), and the actual numerical scheme for the determination of the interaction coefficients of the partly coarse-grained atoms (section 2.6) are deliberately separated. 2.5.1. Construction of the transition region The proposed method is an energy-based method. Therefore, there should exist a total energy functional of the atomistic-continuum system (e.g. see equations (A.1) and (A.12) for the atomistic exchange and the atomistic DM energies, respectively), from which the effective fields for all atoms and nodes are derived by taking the partial derivatives of the total energy by m ~ at the atoms/nodes. To ensure that this is possible, the symmetry of the interaction is enforced for all atoms and nodes. For the atomistic region, these symmetry requirements are formally written in the following way: κ κ Jij = Jji ,

κ κ ~ ij ~ ji D = −D ,

κ ∈ {F, P, A}.

(20)

Obviously, if both atoms, i and j, are real atoms, (20) is fulfilled automatically. The interatomic interaction between two different types of atoms complies with the rules/constraints that are enforced for both types of atoms, i.e. when A either i or j corresponds to a real atom, the exchange coefficient is Jij , while

when either i or j corresponds to a fully coarse-grained atom, the exchange F coefficient is Jij . This formally is written as P A Jij = Jij = Jij ,

P A ~ ij = D ~ ij ~ ij D =D ,

if i ∈ ΩA and j ∈ ΩP ,

(21)

P F Jij = Jij = Jij ,

P F ~ ij = D ~ ij ~ ij D =D ,

if i ∈ ΩF and j ∈ ΩP ,

(22)

18

300

where ΩA , ΩF , and ΩP are the sets of the real, the fully coarse-grained, and the partly coarse-grained atoms, respectively. To ensure that the method is energy-based, the continuum nodes should symmetrically interact with the atoms. This means that if a real atom, which has large interaction range, is located near the atomistic-continuum interface

305

(but inside the atomistic region), and if there is a continuum node, which is located inside the continuum region, and which is within the interaction range of the real atom, such continuum node should interact with the real atom in the same way as the real atom interacts with the continuum node. Technically this is possible, but this breaks the local nature of the continuum description. To

310

keep the continuum description local (strictly according to equations (4) and (5)), the fully coarse-grained atoms are introduced, as shown in figure 1. The boundary nodes are the degrees of freedom of the continuum model and they interact with the nearest transition atoms (which are the fully coarse-grained atoms) as if these nodes are the fully coarse-grained atoms. The interaction coefficients of the fully coarse-grained atoms can be determined from the constraint that these coefficients should lead to the continuum tensor parameters. Thus, these coefficients must satisfy 1X F J ~sij ~sij = Ae , 2 j ij

X

~ F = De , ~sij D ij

(23)

j

which correspond to (6) and (7), respectively. Moreover, the requirement that only the nearest neighbours interact leads to the following constraints:

315

F Jij = 0,

~ F = ~0, D ij

if |~sij | > a and i ∈ ΩF ,

(24)

F Jij 6= 0,

F ~ ij D 6= ~0,

if |~sij | ≤ a and i ∈ ΩF ,

(25)

where a is the lattice spacing and ~sij is the vector connecting atoms i and j. It should be noted that equation (23) does not fully determine the interaction coefficients of the fully coarse-grained atoms, as the number of equations is smaller than the number of unknowns; however, it is always possible to introduce 19

additional constraints, e.g. lattice symmetry, which make the determination 320

possible. For example, in the case of fcc stacking of the atoms and isotropic A exchange interaction (Jij depend only on the distance between atoms i and j), F all Jij corresponding to nearest neighbours can be taken equal (this was done

for the numerical examples of sections 4.3 and 4.4). The partly coarse-grained atoms are introduced because of the symmetry 325

of the interaction as well. If a real atom, which has large interaction range, is placed next to the fully coarse-grained atoms, the symmetry is broken: the real atom can “see” the fully coarse-grained atoms over several lattice spacings, while the fully coarse-grained atoms can “see” only the nearest neighbours (by definition).

330

Due to the introduced constraints, there are two geometric requirements for the atomistic-continuum transition region. The first requirement, as discussed above, is that at least one layer of the fully coarse-grained atoms should isolate the continuum region completely. The second requirement is that the width of the transition region should be larger than the interatomic interaction range of

335

the real atoms, rI . This follows from equations (20), (24) and (25). The total energy of the atomistic-continuum system consists of the “noninteractive” part, i.e. the anisotropy and the external field terms, and the “interactive” part, i.e. the exchange and the DM terms. The former part is trivial, as it does not involve any interactions between atoms/nodes, while

340

the latter part is relevant for the atomistic-continuum coupling method. The interatomic interaction energy is the sum of energies of the interatomic links for the entire atomistic region (including the transition atoms), see equations (A.1) and (A.12). The continuum region was not modified in this description, only the continuum boundary nodes were connected to the nearest fully coarse-grained

345

atoms. Therefore, the “interactive” part of the total energy of the atomisticcontinuum system is the sum of the “interactive” part of the total energy of the continuum region and the “interactive” part of the total energy of the atomistic region, which also includes the energies of the links between the closest to the continuum region layer of the fully coarse-grained atoms and the continuum 20

350

boundary nodes. 2.5.2. Ghost torques As discussed in the introduction, the major problem for all energy-based atomistic-continuum methods is the presence of the ghost forces. In the case of micromagnetic systems, the error in the force corresponds to the error in the

355

effective field. Torque that acts on spin magnetic moments is formed by the ~ Hence the error is referred to as the ghost torque. vector product of m ~ and H. To quantify this error, an error analysis as in sections 2.4.1 and 2.4.2 has to be performed for the partly coarse-grained atoms. The key difference here is that the interatomic interaction coefficients are modified and, for the partly coarse-

360

~ ij become J P and D ~ P , respectively, while for the fully grained atoms, Jij and D ij ij ~ ij become J A and D ~ A , respectively. coarse-grained atoms Jij and D ij ij At first, the exchange terms are considered. The effective field of the partly ~ EAP . In addition to this, H ~ EA is introcoarse-grained atoms is denoted as H i i duced, which is the effective field of atom i as if it is completely surrounded by the real atoms. For clarity of the explanation, it should be noted that these quantities are expressed as ~ iEAP = H

X

P Jij m ~ j,

j

~ EA = H i

X

A Jij m ~ j.

j

To find the local error in the effective field of the partly coarse-grained atoms in ~ EAP and H ~ EA should be compared for a given atomistic the multiscale system, H i i solution m ~ i. As mentioned above, in the partly coarse-grained region, the interaction coefficients and the interaction range are modified. Referring to figure 1, each partly coarse-grained atom interacts with fewer atoms on the left than on the P P right. Therefore, although the interaction symmetry is preserved (Jij = Jji ), the

lattice symmetry is broken: there might not exist atom k, such that ~sik = −~sij

P P and Jik = Jij . Therefore, all terms from equation (9), except the first term,

should be used for the error analysis. Using equations (9) and (6), the following

21

result can be obtained: ~ iEAP − H ~ iEA = H

X j



 X 1 P Jij ~sij · ∇m ~ + J P ~sij ~sij − Ae  : ∇∇m ~ + 2 j ij   1X P 3 + Jij (~sij · ∇) m ~ + O a4 Jˆ ∇4u m ~ . (26) 6 j

In the proposed multiscale scheme, it is suggested to determine the interaction coefficients of the partly coarse-grained atoms based on the minimisation of the error in the effective field. The lowest-order term in (26) can be removed by imposing the following constraint: X

P Jij ~sij = ~0,

j

∀i ∈ ΩP ,

(27)

where ΩP is the set of the partly coarse-grained atoms. Using this constraint, it is possible to obtain an estimate for the difference between the atomistic and the continuum representations of the effective field:    X  1 P ~ GE = H ~ EA − H ~ EAP = Ae − Jij ~sij ~sij  : ∇∇m ~ + O a3 Jˆ ∇3u m ~ , H i i i 2 j

(28)

where, as in section 2.4.1,

∇pu

is the directional derivative along the unit vector

~u, the direction of which is selected such that |∇pu m| ~ takes the maximum value. The ghost torque, which acts on atom i, is formed at the right hand side of the LLG equation when the effective field is vector multiplied by m ~ from the left. Therefore, the ghost torque is ~ iGE , T~iGE = m ~ ×H 365

(29)

where m ~ and its derivatives are taken at lattice position i. This is an error that arises due to the representation of the transition region as it is proposed here (using the partly coarse-grained atoms) rather than as a real atomistic region. P It can be quantified once coefficients Jij are calculated (see also discussion in

section 4.2.2).

22

For a large number of the partly coarse-grained atoms, the number of equaP ). Therefore, there tions in (27) is smaller than the number of unknowns6 (Jij

still is some freedom in the determination of the unknowns. For a further minimisation of the ghost torques, the tensorial prefactor PE i = Ae −

1X P J ~sij ~sij , 2 j ij

(30)

which appears in front of ∇∇m ~ in (28), can be minimised. It is convenient to

PF introduce the following notation: Jij is the interatomic interaction coefficient

between partly coarse-grained atoms i and j, which is calculated as if atoms i PA and j are the fully coarse-grained atoms; in a similar way, Jij is the interatomic

interaction coefficient as if atoms i and j are the real atoms. These coefficients are related to Ae by equations (23) and (6), respectively7 : 1 X PA 1 X PF Jij ~sij ~sij = J ~sij ~sij = Ae . 2 j 2 j ij Hence, the prefactor P E i can be rewritten as   1X ζ1 ζ2 E PF PA P Pi = J + J − Jij ~sij ~sij , 2 j ζ1 + ζ2 ij ζ1 + ζ2 ij 6 This

(31)

(32)

can be easily seen in 1D. If the real atoms have interaction range of na, then n partly

coarse-grained atoms have to be used. Therefore, the number of constrains (27) is n. However, it can be shown that the rank of system (27) is n − 1, which results in n − 1 independent constraints. The number of unknowns is the number of interactions of each partly coarsegrained atom with all other partly coarse-grained atoms, taking into account the interaction 2 = n (n − 1) /2. If n > 2, system (27) is underdetermined. The case of symmetry, which is Cn

n = 2 is a special case for which the system can be solved exactly. 7 Here it is implied that the sum in (31) includes as many terms as for the real atoms. However, the partly coarse-grained atoms that are located close to the atomistic-continuum interface have fewer atomistic neighbours within the interaction range of the real atoms (same is valid for the fully coarse-grained atoms). This can be resolved by introducing virtual atomistic sites, which can be located inside the continuum region, as the pad atoms are; however, these virtual sites interact with neither the partly nor the fully coarse-grained atoms, κ = 0, i ∈ Ω , κ ∈ {F, P} and j is the index of the virtual site. As is shown later, these i.e. Jij κ

sites do not participate in any calculations.

23

370

where ζ1 and ζ2 are arbitrary positive numbers, meaning of which is explained in section 2.6. There is some ambiguity in tensor minimisation, as it is possible to minimise its components or a particular invariant. Moreover, there is a number of prefactors, since there is a number of the partly coarse-grained atoms. The exact scheme of minimising the prefactors is presented in section 2.6. In the case of the DM interaction, a similar analysis leads to the following ~ P: constraint for D ij

X

P ~ ij D = ~0,

∀i ∈ ΩP .

j

(33)

The ghost torque for the DM part of the effective field is given by    X  ~ GD = H ~ DA − H ~ DAP = ∇ · D e − ~ P × m ˆ ∇2 m . (34) H ~sij D ~ + O a2 D i i i ij u~ j

~ GD , T~iGD = m ~ ×H i

(35)

~ P can be determined by minimising the following tensorial prefactor: Hence D ij PD i = De − 375

X

P ~ ij ~sij D .

(36)

j

It can be seen that when the material parameters at the continuum scale (Ae and D e ) are fixed, while lattice spacing (a) varies, the ghost torques (29) and ~ ij become O(1/a2 ) (35) are O(1) with respect to a, since in this case, Jij and D and O(1/a), respectively. Thus, constraints (27) and (33) eliminate the ghost torque terms that are O(1/a).

380

2.5.3. Relation to atomistic-continuum modelling of deformation The energy-based quasicontinuum (QC) method [13, 14] is one of the most well-known atomistic-continuum coupling methods for mechanostatic problems. In the case of the QC method, it is also possible to separate the coarse-graining error, which is the local force error inside the continuum region due to the change

385

of the description from the atomistic to the continuum, and the local force error at the interface (the ghost forces). The QC method uses the Cauchy-Born rule for the continuum model, which was analased in [37], see also [22] and references 24

therein. This model has the coarse-graining error O(a2 ), where a is the interatomic spacing. The classical micromagnetic continuum description also has the 390

coarse-graining error that is O(a2 ), see (13) and (17). However, in contrast to the Cauchy-Born rule, the micromagnetic model is not fully consistent with the atomistic description in the case of the patch test, see section 3.2. In [38], it was discussed that the original formulation of the QC method suffers from the local force error at the interface that is O(1/a), while QC-

395

QNL schemes, in which the ghost forces are eliminated under certain conditions and for the uniform deformation, e.g. [17], have the local force error O(1) with respect to a. Using a 1D test case, in [38], it was shown that the lowest achievable local force error for any energy-based QC-type method is O(1), if the interface has a bounded number of degrees of freedom.

400

Results of [38] are directly applicable to magnetostatic problems, however, only in a 1D spatial case (e.g. chain of atoms) and when m ~ is constrained to a plane8 . In this case, mechanostatic and magnetostatic problems are mathematically equivalent, while the force in the mechanostatic formulation corresponds to the torque in the magnetostatic formulation, which is shown in section 3.2. As

405

it is possible to create such magnetisation field in 3D, the lowest achievable local torque error (the ghost torques) for any energy-based atomistic-continuum coupling method for magnetostatics is also O(1) with respect to the lattice spacing. This is achieved by the method proposed in this article. Also, in [39], it has been shown that although the local force error at the

410

interface for QC-QNL scheme is O(1), the solution error is O(a1+1/p ) in the W 1,p norm. In this article, the solution error is analysed numerically in section 4.2 for the dynamic case. As mentioned in the introduction, the proposed atomistic-continuum transition method is similar to the geometric reconstruction method for mechanostatic

415

problems, which was first introduced in [18], and which was further extended to non-flat geometries in [19, 20], where it was denoted as GRAC. This method sug8 If

m ~ is not constrained to a plane, the specific form of the energy functional is different.

25

gests to modify the interatomic interaction for atoms that are located close to the atomistic-continuum interface. As these atoms lack a number of neighbouring atoms, for each atom in the interface neighbourhood, a geometric reconstruc420

tion technique is used to obtain the positions of the missing neighbours from the positions of the existing neighbours by extrapolation. Following this, the same interatomic interaction potential function as in the atomistic region (which is an arbitrary three times continuously differentiable function), is applied to these atoms. The reconstruction technique uses a number of free (extrapolation) pa-

425

rameters. In this article, it is also suggested to introduce a region consisting of special atoms (the transition atoms) next to the interface, for which the interatomic interaction is modified. However, instead of performing reconstruction of missing neighbours, in this article, it is suggested to modify the interatomic interaction function between the existing neighbours directly. As in the case of

430

spin statics/dynamics the interaction types are known (the exchange and the DM interactions), it is suggested to use the same functional dependencies, but with different coefficients, which become free parameters. Thus, both schemes modify the interatomic interaction of the atoms that are closer to the interface than the interaction range of the real atoms, both schemes use the solutions at

435

the lattice points of the existing neighbours of the transition atoms to perform this modification, and both schemes have a number of free parameters which have to be determined in a preprocessing step. In the GRAC scheme, two constraints for the reconstruction parameters are introduced: for the uniform deformation case, the local energy consistency

440

(the transition atoms must have the same energy as the real atoms) and the local force consistency (the forces that act on the transition atoms must be zero) are enforced, e.g. see equations (3.3) and (3.9) in [19]. In this article, constraints (27) and (33) are enforced, which are the analogies of the local force consistency; however, in this article, the ghost torques are not removed

445

completely for a uniformly varying magnetisation, the error is only reduced (see discussion in section 3). Constraints (27) and (33) enforce the lowest-order terms of the ghost torques of the exchange and the DM interactions, respectively, to 26

be zero. In both methods, the number of constraints is lower than the number of unknowns. Therefore, both methods suggest to solve a constrained least 450

square problem to find the parameter values, e.g. see equation (2.13) in [20] and see section 2.6 of this article. Since in the GRAC scheme, the ghost forces are removed completely only for the uniform deformation case, i.e. there is an error in the force for a non-uniform deformation, the minimisation problem for the reconstruction coefficients is designed to minimise the prefactor of the

455

lowest-order term of the force error (although with some simplifications). In this article, the minimisation problem (40) is constructed to minimise the prefactor entering the remaining lowest-order terms in the torque errors. 2.6. Interaction coefficients of transition atoms P ~ P satisfying equations (20), (21), (22), (27), and In general, any Jij and D ij

460

(33) can be used, as the number of constrains is smaller than the number of unP ~ P is underdetermined, knowns. Since the problem of finding unknown Jij and D ij

in section 2.5.2, it was suggested that one of the possible ways of constructing P ~ P is the minimisation of tensorial prefactors (30) and (36), which enter Jij and D ij

the lowest-order terms in the ghost torques, given the above constraints. Only 465

P the determination of Jij is discussed below, since a similar technique can be

~ P. applied to unknown D ij Within this scheme, it is proposed to consider the trace of tensor P E i of each partly coarse-grained atom for the minimisation problem. Using (32), the following estimate can be obtained:  1 X  ζ ζ2 1 PF PA P = J + J − J tr (~ s ~ s ) ≤ tr P E ij ij i ij ij ij 2 j ζ1 + ζ2 ζ1 + ζ2 ! 1 X PF PA ζ1 + Jij ζ2 1 X P Jij 2 zij a2 , (37) = |~sij | = Jij − 2 2 j ζ1 + ζ2 j

where the following quantities are introduced:
P ∗ zij = Jij − Jij wij ,

∗ Jij =

PF PA Jij ζ1 + Jij ζ2 , ζ1 + ζ2

27

2

wij =

|~sij | . a2

(38)

There is a number of the partly coarse-grained atoms, hence a number of tensorial prefactors P E i ; however, the minimisation problem requires a single scalar quantity. Each individual component of the sum in (37), which is zij a2 , 470

P corresponds to a contribution of a particular Jij to the torque error of atom i.

Therefore, it is suggested to take the L2 norm of vector-column containing all zij as the scalar quantity for minimisation. This quantity corresponds to the global ghost torque prefactor error due to the interaction coefficients between the partly

475

coarse-grained atoms being different to the interaction coefficients of the real P 2
1/2 P zij , atoms. Thus, to determine unknown Jij , it is suggested to minimise where the sum is over all pairs of ij for which interaction is unknown, with constraints (20), (21), (22), (27). F It is useful to introduce rij as the distance between point (~xi + ~xj ) /2, which

is the middle of the link between atoms i and j, and the closest fully coarseA grained atom. In a similar way, rij is a distance between the middle of the link

between atoms i and j and the closest real atom. Since parameters ζ1 and ζ2 can be arbitrary positive numbers, by selecting A ζ1 = rij ,

F ζ2 = rij ,

it is possible to obtain a smooth transition between the real atomistic interaction ∗ coefficients and the fully coarse-grained interaction coefficients. In this case, Jij 480

A F is close to Jij and Jij when both atoms (i and j) are next to the real and the

fully coarse-grained regions, respectively. It is convenient to rewrite the minimisation problem in a matrix-vector form. In this case, constraints (20), (21), (22), (27) can be rewritten as M J = B,

(39)

P where J is a column containing all unknown Jij . In order to write (39), at

first, constraints (20), (21), (22) are used to reduce the number of unknowns. Afterwards, elements of matrix M and column B are obtained from equation 485

A F (27). Hence, column B contains Jij and Jij multiplied by sijp , which are the

components of ~sij , while matrix M contains various sijp . Since the problem 28

is underdetermined, the number of rows in the matrix M is smaller than the number of columns. ∗ Columns Z, J ∗ and diagonal matrix W , which contain zij , Jij , and wij ,

respectively (the same indexing as in equation (39) is implied), are introduced. Finally, the unknown J is obtained by minimising kZk given constraints (39). This problem has a unique solution, which can be explicitly written as T T J = W −1 MW MW M W


−1

(B − M J ∗ ) + J ∗ ,

(40)

where MW = M W −1 . Of course, solution J can be obtained using one of the 490

efficient iterative methods instead of evaluating (40) directly. 2.7. Wave-absorbing band at the interface When dynamic behaviour of a material is considered, the major problem for the multiscaling technique concerns reflections of high frequency waves from the interface separating two scales. In relation to magnetic systems, this problem was partly discussed in [40]. In [33], a wave-absorbing band for magnetisation dynamics was introduced for the case of continuum domains with different spatial discretisations. In this article, a wave-absorbing band is used in the atomistic region at the proximity of the interface. Since there is already a modified region (transition atoms) within the atomistic region, it is convenient to introduce the additional damping within this region, although, obviously, an arbitrary number of atomic layers can be used for the damping band. Thus, the equation describing the dynamic behaviour of transition atoms becomes the following [33]:  
d ~ i − αL m ~ i − γiD m m ~ i = −βL m ~i×H ~i× m ~i×H ~i× m ~i×m ~A i , dt

(41)

where γiD is the strength of the damping for atom i and m ~A i is the normalised average of the solution over a certain region. The last term in the equation (41) “forces” m ~ i and m ~A i to become parallel acting effectively as a low-pass filter. The damping strength increases gradually with atomic positions approaching the interface. Moreover, the damping strength is proportional to the time 29

495

derivative of m ~ i . The following expression for the damping strength is used: 2 s  d riR D m γi = gD (42) ~ i , P R dt ri + ri − a where gD is a constant determining the damping strength, riR is the distance to

the closest atom outside of the damping band (real atom), riP is the distance to the closest interface atom/node or pad atom and a is the interatomic spacing. This relation is a modified expression for γ D from [33], where a 1D model was considered. In this article, the distance riR is scaled. The reader is referred to 500

[33] for the explanation of the specific structure of equation (42). Within the damping band, spin magnetic moments of atoms are damped to an average value, m ~A i , which is calculated over a certain region (the averaging window surrounding atom i), which has the width and the height of sA . Within this region, an additional quantity m ~ A∗ is introduced: i m ~ A∗ cxx x2 + ~cyy y 2 + ~cxy xy + ~cx x + ~cy y + ~c0 . i =~

(43)

Coefficients ~c are fitted to all m ~ j within the averaging window using least squares. Finally, m ~A i is obtained by normalisation: m ~ A∗ i m ~A i = A∗ . m ~i

(44)

As mentioned above, such damping band is characterised by two parameters: gD and sA . These parameters depend on the difference between the continuum spatial step and the atomistic lattice spacing and can be fitted numerically by simulating spin wave propagation in a 1D atomistic-continuum system, for de505

tails see [33]. Moreover, in [33], it was shown numerically that there are optimal values of these parameters, which minimise the error from wave reflections. In particular, small gD or small sA lead to inefficient damping, thus there are wave reflections from the atomistic-continuum boundary, while large gD or large sA lead to a very strong damping, thus wave reflections from the boundary between

510

the damping layer and the real atomistic region dominate. In addition to this, there is one more geometrical parameter — the width of the damping band. Although in the beginning of this section, it was suggested to make the damping 30

band as wide as the transition region, the width of the damping band can be arbitrary and the wave reflection error decreases with the increase of the width 515

of the damping band. 2.8. Time-stepping scheme In the case of dynamic problems, the energy-conserving property of the numerical method is useful. The total energy of the system is preserved when energy-conserving time-stepping numerical method, such as the mid-point method

520

[41], is used. The mid-point class of methods conserves the Hamiltonian structure of the LLG equation, the Lyapunov structure of the LLG equation, length of magnetic moments and is unconditionally stable [35]. This numerical method requires the symmetry conditions enforced for the interactions between all atoms/nodes, see [41] for details. Since such symmetry conditions are prescribed for the mul-

525

tiscale method proposed in this paper, see section 2.5.1, the application of the mid-point method for the time-integration of the dynamic equations results in the conservation of the total energy. The proposed multiscale method is energy-conserving only when the continuum mesh is refined to the atomistic scale at the interface9 , as illustrated

530

in figure 1. This can be achieved by representing the continuum using finiteelements10 . If the continuum mesh is not refined at the interface, which is the case when the finite-difference continuum discretisation is used, then pad atoms11 should be introduced, which breaks the symmetry required for the energy conservation of the mid-point scheme (there can exist atom k which “sees” 9 This

means that arbitrary geometry of the atomistic-continuum interface can be used as

well as arbitrary continuum mesh, as long as the continuum boundary nodes exactly coincide with the atomistic lattice sites and the distance between neighbouring boundary nodes is not more than a. 10 For the finite-element representation of the continuum region, the energy-conserving (in the case when damping is absent) implicit numerical method by Bartels and Prohl [42] can be used. 11 Pad atoms are atoms, magnetic moment of which is entirely determined by interpolation between neighbouring continuum nodes.

31

535

continuum node q, while node q does not “see” atom k). Obviously, in 1D, the refinement of the mesh is not necessary. Moreover, the mesh refinement is not necessary if the numerical integration scheme is not energy-conserving in the first place (e.g. the backward Euler method); however, the use of the transition atoms still reduces the error due to the coupling of local (continuum) and

540

non-local (atomistic) descriptions, as was shown in section 2.5.2. The atomistic and the continuum regions can either be solved together implicitly (with the same time step), which is possible when the continuum mesh is refined to the atomistic scale, or decoupled, which allows using different time steps for the regions and is more convenient in terms of implementation, but in-

545

troduces an error in the total energy of the system, i.e. the method stops being energy-conserving. For this reason, in the cases when the energy conservation is essential, a decoupled scheme is not advised. However, in the case of 2D, it is sometimes convenient to use a continuum mesh that is not refined to the atomistic scale, which already introduces an error in the total energy. In this case,

550

the non-conservation of the total energy by the decoupled time-stepping scheme does not seem to be a significant disadvantage. Therefore, a possible way of constructing a decoupled scheme is discussed below. The scheme is similar to the one discussed in [43]. For simplicity and clarity of the explanation, a 1D case is discussed.

555

The mid-point method [41], which is used in this article, is implicit and, in the case of an atomistic-continuum interface, it requires values of m ~ I (tn+1 ) at the interface points to advance from time step tn to tn+1 . These values depend on the state of the atomistic region at tn+1 . Therefore, an extended continuum grid is assumed to fill the entire computational region, including the

560

atomistic region; however, nodes that are located in the atomistic region are treated differently. These nodes are referred to as the auxiliary nodes, see figure 2. The time evolution now takes the following steps: (a) At a given time step tn , the solution at these auxiliary nodes is obtained by volume-averaging the atomistic solution with a subsequent normalisation.

32

Figure 2: The schematic representation of the numerical time-stepping scheme in 1D. Atomic positions are marked as open circles, continuum nodes are marked as open squares, and auxiliary nodes are shown as crosses. Nodes, where the solution is provided by the boundary conditions and the interpolation, are shown with solid squares and solid circles, respectively. The sequence of computational steps is indicated with numbers.

565

(b) Using the extended mesh and the boundary conditions for the entire region, two time steps are performed, and the solutions at tn+1 and tn+2 are obtained. (c) Using values of m ~ I (tn−1 ), m ~ I (tn ), m ~ I (tn+1 ) and m ~ I (tn+2 ) at the interface 570

points, interpolated values m ~ I (τ ), τ ∈ (tn , tn + ∆t) are obtained. The values of magnetic moments at the interface are interpolated in time using circle splines [44]. In this case, the lengths of magnetic moments are naturally preserved. (d) These interpolated values serve as boundary conditions for the atomistic

575

region, which is advanced with its own time step up to tn+1 . (e) At tn+1 , all values at the continuum auxiliary nodes are overwritten by a new volume-averaging, see step (a). The continuum solution at tn+2 that was obtained in step (b) is disregarded. This scheme can be straightforwardly generalised to the case when the con-

580

tinuum time step is not divisible by the atomistic time step. Moreover, the continuum auxiliary node values in step (a) can also be improved by an additional time-averaging. In the case of 2D, during step (c), interpolation of the 33

pad atoms is performed as well. The values are interpolated between several nodes and this can be performed by polynomial interpolation with a subsequent 585

normalisation.

3. Analytical examples: 3D and 1D patch tests (statics) The expression for the ghost torques of the proposed atomistic-continuum coupling method for arbitrary m ~ was obtained in section 2.5.2. The aim of this section is to analyse this error further for the specific test case of a uniformly 590

changing magnetisation, as this case is analogous to the test problem that is often utilised in testing of atomistic-continuum coupling methods for mechanical systems. Moreover, in this section, an equivalence with the mechanical description in 1D case is established and an important consequence of the selection of the continuum model is highlighted.

595

The test that is often utilised to access the accuracy of atomistic-continuum methods in application to mechanics is the uniform deformation problem [22], which is sometimes called the patch test. In this case, the error is expressed as ghost forces, which appear at the atomistic-continuum interface of a uniformlydeformed configuration. If the energy minimisation problem is solved then these

600

forces lead to a modified non-uniform strain, hence an error in the solution appears. As mentioned in the introduction, there are atomistic-continuum methods in application to mechanics that are free of these artefacts in the case of a uniform deformation. The static problem of finding magnetisation distribution in a magnetic ma-

605

terial by minimisation of the total energy is analogous to the mechanical problem of finding deformed equilibrium configuration. Moreover, in a 1D spatial case with spin magnetic moments constrained to the plane, a mathematical equivalence can even be established. Therefore, the representative test for an atomistic-continuum method for magnetic materials is the analogy of the uni-

610

form deformation test, which is the state with the constant angle between neighbouring atomic spins.

34

In this section, only the exchange interaction is considered. The analogy of a patch test for a magnetostatic problem with the DM interaction is somewhat ambiguous. Therefore, the test, which is considered in this section, is referred 615

to as the exchange patch test. As presented in sections 2.2 and 2.3, the solutions of the spin dynamics and the micromagnetic models are constrained to a spherical manifold, i.e. the solution m ~ (~x) is a vector field such that |m| ~ = 1. This means that even for the case of a uniformly varying magnetisation as a function of angle (exchange patch test,

620

e.g. 1D spin spiral with a uniform twist), the solution m ~ depends non-linearly on ~x. This has two consequences. The first consequence is that if a similar non-linear solution is considered in a mechanical atomistic-continuum problem, methods that are free of the ghost forces (in the case of uniform deformation), such as the methods based on the geometric reconstruction [20], will result in an

625

error in the force. Although, the GRAC scheme [20] could probably be modified such that the reconstruction of the solution is performed on a spherical manifold (reconstruct the angles between vectors, rather than vectors themselves). This might eliminate the ghost torque problem for the exchange patch test case and is an interesting problem for a future research.

630

The second important consequence of the non-linearity of the exchange patch test solution is that the continuum description should be also non-linear and preserve the non-linearities of the atomistic description. However, there is an intrinsic mismatch between the atomistic and the continuum models — different torques are required to create the same uniform twist in the atomistic and the

635

continuum domains. This is demonstrated in section 3.2.1, where spin spiral solutions are considered. In solid mechanics, the mismatch between non-linear atomistic and continuum models, for the case of the uniform deformation, was removed by constructing a Cauchy-Born continuum representation [16]. Therefore, a similar construction of continuum model for magnetism is suggested in

640

section 3.3.

35

3.1. 3D exchange patch test As mentioned above, a uniform deformation in mechanics corresponds to a magnetostatic configuration, for which angles between any equidistant spin magnetic moments along any given direction are constant. This can be expressed as m ~ (~x + ~c) × m ~ (~x) + m ~ (~x − ~c) × m ~ (~x) = ~0,

∀~x, ~c ∈ R3 .

(45)

If m ~ is expressed using spherical coordinates (i.e. using polar and azimuthal angles), equation (45) is fulfilled, for example, when the azimuthal angle is constant, while the polar angle depends linearly on the spatial coordinates. The ghost torque, which is given by equation (29), can be rewritten as T~iGE =

3 3 X X p=1 q=1

E Ppq m ~ ×

  ∂2m ~ ~ , + O a3 Jˆ ∇3u m ∂xq ∂xp

(46)

E are the components of the prefactor tensor P E where Ppq i , while xp are the

components of ~x. The vector products inside the sum simplify to  ∂2 1 m× ~ m ~ =m ~ (~x)× lim m ~ (~x + h~ep + h~eq ) − m ~ (~x − h~ep + h~eq ) − h→0 4h2 ∂xq ∂xp   1 m ~ (~x) × m ~ (~x + h~c+ ) + −m ~ (~x + h~ep − h~eq ) + m ~ (~x − h~ep − h~eq ) = lim h→0 4h2  +m ~ (~x) × m ~ (~x − h~c+ ) − m ~ (~x) × m ~ (~x + h~c− ) − m ~ (~x) × m ~ (~x − h~c− ) = ~0, (47) where ~c+ = ~ep + ~eq and ~c− = ~ep − ~eq . The vector products are cancelled out by

645

equation (45). Thus, for the exchange patch test case, the ghost torque is  
T~iGE = O a3 Jˆ ∇3u m ~ = O aA∗ ∇3u m ~ . (48) This means that for fixed m ~ and Ae , the ghost torque for the exchange patch

test is O(a). It should be noted that this is only the local error in torque, which is acting on spin magnetic moment of atom i, given a patch test solution for m. ~ 3.2. 1D problems The purpose of this section is to show that there is an intrinsic mismatch 650

between the atomistic and the continuum descriptions of magnetism at the 36

boundary of the computational domain even for a uniform twist. This mismatch is calculated for a 1D case, although, obviously, the mismatch persists in 2D and 3D cases. For this purpose, two half-infinite problems are considered. After this, a 1D atomistic-continuum transition problem is analysed for the comparison. For a 1D chain of atoms interacting via the Heisenberg exchange, atomic spin directions can be expressed as m ~ i = ~ex sin φi + ~ey cos φi ,

(49)

since, without loss of generality, one degree of freedom per spin can be used for this problem. The total energy of the atomistic system is written in the following form (see equation (A.1) for the exchange energy expression): EA = −

i+N 1X X J|j−i| cos (φi − φj ) , 2 i

J0 = 0.

(50)

j=i−N

Here, double-indexed Jij was replaced by single-indexed Jn , which is an exchange constant between particle i and particles i ± n. Particle i interacts with 2N particles. Spin statics problem is written in the following way (right hand side of equation (1) is zero): ~ i = ~0, m ~i×H

A

~ i = − ∂E , H ∂m ~i

∀i.

(51)

Using (49), equations (51) simplify to TiA = 0, 655

TiA = −

i+N X ∂E A J|j−i| sin (φi − φj ) , =− ∂φi

(52)

j=i−N

where TiA is a scalar value of a torque12 , which is acting on spin magnetic moment m ~ i . Here the equivalence with the mechanical elastostatic 1D problem can be seen, in which forces are equated to zero and are given by Fi = −∂E/∂ui , where ui is the displacement of particles. For the purpose of this example, the continuum description is discretised using a second-order finite-difference discretisation13 with the step size ∆x. In the 12 In 13 It

~ i is parallel to ~ez , while T A is the projection of m ~ i on −~ez . this example, m ~ i ×H ~ i ×H i can also be shown that a higher-order discretisation does not solve the problem high-

lighted in section 3.2.1.

37

case of static 1D problem with one degree of freedom, magnetostatic equations (right hand side of (4) equals zero) can be rewritten in a scalar from, as was done in the atomistic case. The scalar torque at node k can be obtained: TkC = −Ae

1 (sin (φk − φk−1 ) + sin (φk − φk+1 )) , ∆x2

Ae = a2

N X

Jn n2 ,

(53)

(54)

n=1

where a is the interatomic spacing. Equation (54) follows from equation (6). 660

Torques TkC and TkA are exactly the same for the case when ∆x = a and N = 1, which is the case of the non-zero interatomic interaction between the nearest neighbours only. 3.2.1. Separate half-infinite problems Two half-infinite 1D geometries are considered: the discretised continuum and the atomistic problems, see figure 3ab. A spin spiral with a uniform twist, for which φi = εai, is considered. To avoid an additional discretisation error, ∆x = a is taken. To equilibrate continuum domain solution, a compensating torque, which is analogous to a reaction force in a mechanical system, should be added to node 0: T CR = T0CR = −Ae

N X 1 sin (εa) = − Jn n2 sin (εa) . ∆x2 n=1

(55)

When the atomistic chain is considered, the compensating torques should be added to N initial atoms: TnAR = −

−n−1 X

j=−N

J|j| sin (−jεa) ,

n = 0, . . . , N − 1.

(56)

Therefore, the total added torque is T

AR

=

N −1 X n=0

TnAR

=−

N X

Jn n sin (nεa) .

(57)

n=1

The mechanical equivalence explains why atomic forces/torques are summed in 665

(57). 38

Suppose an atomistic-continuum multiscale problem is constructed by combining these two half-infinite geometries. In this case, the mechanical equivalence helps to see that if a system, which is constructed out of two separate parts, is to be equilibrated, these parts should be deformed by the same force/torque. However, it can be seen that when the half-infinite atomistic and continuum problems are considered separately, a slightly different torque is required to achieve the identical solution:   X   N N 3 3 3 3 3 CR X

ε a n ε a AR 2 5 5 5 5 T Jn n εa − −T = +O a ε +O a ε − Jn n nεa − = 6 6 n=1 n=1 N ε3 a3 X  
Jn n2 n2 − 1 + O a5 ε5 Jˆ . (58) = 6 n=2

This means that it is not possible to combine these two separate descriptions within the partitioned-domain energy-based approach, such that the identical

constant gradient of φ is maintained in both regions. The imbalance in internal forces/torques will lead to an error in the solution. The error in the torque 670

is of the third order with respect to the gradient, ε, and will be referred to as the boundary mismatch error. For the fixed continuum parameters and the fixed continuum solution, the boundary mismatch error becomes O(a), since
Jn = O 1/a2 . In this case, the aim of the coupling method is not to increase the error, which is already present in the multiscale system due to this mismatch.

675

Since the boundary mismatch error does not disappear in 3D, and the ghost torque is of the same order of magnitude with respect to the lattice spacing for the 3D patch test, it can be inferred that the proposed multiscale method reduces the error in torques down to the intrinsic boundary mismatch error between the atomistic and the continuum descriptions for the case of uniformly

680

changing m, ~ see (45), and the exchange interatomic interaction. However, in the case of arbitrary magnetisation field, the lower-order terms appear, and the torque error is O (1), see equations (28) and (29). The boundary mismatch error (58) is different to the coarse-graining error, which is estimated in section 2.4.1. First of all, the coarse-graining error is

685

defined as the local truncation error in the effective field at the continuum nodes 39

Figure 3: The schematic representation of the half-infinite 1D discretised continuum geometry (a), the half-infinite 1D atomistic geometry (b) and the infinite 1D multiscale geometry (c), which consists of the continuum and the atomistic regions. In (a) and (b), the compensating torques that are applied to nodes/atoms are schematically represented by the arrows with different lengths (this is only a schematic representation, as for this problem, the compensating torques are the vectors that are oriented out of plane of the image). In (c), the distinction is made between the fully coarse-grained atom (brown), the partly coarse-grained atoms (blue) and the real atoms (red).

which are inside the computational domain, where the symmetry is present. The boundary mismatch error is at the edge of the computational domain, where the symmetry is broken. Secondly, the coarse-graining error was defined as a local

690

~ Here, the compensating torques are summed. For the half-infinite error in H.   chain, the local error is O aεJˆ , which is easily seen, for example, for atom (N − 1), since (10) is not fulfilled at the edge of the domain (due to the absence of the symmetry). 3.2.2. The multiscale problem An infinite 1D problem is considered, in which the domain is split into the

695

continuum and the atomistic regions, see figure 3c. The solution with a uniform gradient is considered: φi = εai (spin spiral with a uniform twist). The atomistic-continuum interface is constructed such that there is one fully coarsegrained atom and N partly coarse-grained atoms. According to equation (23), the fully coarse-grained atom interacts with nodes 1 and (−1) with coefficients

700

F F = J0(−1) = Ae /a2 . J01

In this case, the torque which is acting on all nodes up to (−2) is zero (by the symmetry of the discretised continuum description), Tk = 0, k ≤ −2. Interface 40

node (−1) interacts with atom 0 and with node (−2) with the same coefficients, therefore, T(−1) = 0. Fully coarse-grained atom 0 interacts only with interface 705

node (−1) and partly coarse-grained atom 1 and the interaction coefficients are equal, therefore, T0 = 0. For atoms starting from N + 1 (real atoms), the torque is also zero (by the symmetry of the atomistic description), Tk = 0, k ≥ (N + 1). The torque which acts on transition atoms i = 1, . . . , N is the following: Ti = − =− =ε

i+N X j=0

i+N X j=0

i+N X j=0

P Jij sin (φi − φj ) = −

i+N X j=0

P Jij sin (εa (i − j)) =

3

P Jij


(i − j) ε3 a3 + O a5 ε5 εa (i − j) − 6

P Jij a (j − i) +

!

=

i+N   ε3 a3 i+N   X ε3 a3 X P 3 3 P Jij (i − j) + O a5 ε5 Jˆ = Jij (i − j) + O a5 ε5 Jˆ . 6 j=0 6 j=0

(59) The sum which enters O (ε) term is zero due to constraint (27), which is an important feature of the method. It can be seen that the third-order derivative 710

of m ~ in (48) results in ε3 (due to m ~ being expressed via trigonometric functions of φi , which linearly depend on ε). When equations (58) and (59) are compared, it can be seen that the error in the 1D multiscale problem is indeed of the same order of magnitude as the boundary mismatch error with respect to the gradient, ε, and the interatomic

715

spacing, a. 3.3. Cauchy-Born continuum model for exchange interaction As shown in section 3.2.1, the continuum model has a mismatch error at the boundary when compared to the atomistic model. In the case of modelling of the mechanical behaviour of crystals, such problem was solved by introducing a Cauchy-Born continuum model, which is fully compatible with atomistic models for the case of constant deformation gradients. Similar model can be constructed for the case of the exchange interaction between spins. It can be derived from the principle that when equation (45) is fulfilled, the atomistic and the continuum

41

models should have exactly the same values of the exchange parts of the effective fields, which are given by the first terms of equations (3) and (5), respectively. It can be shown (see appendix Appendix C), that this leads to the continuum exchange effective field in the form of ~ EC = H

X j

Jj ~sj · ∇m ~

sin |~sj · ∇m| ~ . |~sj · ∇m| ~

(60)

This model is local, as the effective field at any spatial point depends only on the gradient of m ~ at the same point. However, it has built-in information on the crystallographic lattice and the interatomic interaction coefficients. Here Jj is 720

the interaction coefficient of an atom with the neighbour, to which the atom is connected by a vector ~sj . These coefficients can also be understood as Jij and ~sij , which were introduced previously, but for which index i is omitted based on the assumption that all atoms in the lattice are of the same type. This model can be useful in future for constructing exchange patch test com-

725

patible atomistic-continuum multiscale methods. It is not, however, useful for the approach proposed in this article, since the construction of the fully coarsegrained atoms relies on the form of the continuum model, and in this article, they were constructed for the classical model, which is presented in section 2.3. This means that the fully coarse-grained atoms also contain the boundary mis-

730

match error — different torques are required to create the same uniform twist in the chain consisting of only real atoms and in the chain consisting of only fully coarse-grained atoms. The patch test compatible atomistic-continuum method should have a different representation of the transition region that relies on (60).

4. Numerical examples 735

In this article, the emphasis is made on the numerical technique, therefore, all quantities, including x and t, are considered to be dimensionless, which is common in the field of numerical analysis [35]. The purpose of this section is to investigate the accuracy and the limitations of the proposed method numerically, as functions of relative quantities, such as the number of atoms per domain wall. 42

740

The time-stepping was performed using the mid-point method [41] in both atomistic and continuum models. In the continuum numerical model, the secondorder central finite-difference method was used to discretise the equations spatially. Several cases are considered in this section. In section 4.2, a number of 1D

745

fully atomistic simulations are compared to 1D fully continuum simulations as well as to simulations, where the computational region was split into the atomistic and the continuum regions. In the latter case, the total number of degrees of freedom was kept the same as in the atomistic simulations, i.e. ∆x = a. Such comparison gives an error between the atomistic and the atomistic-continuum

750

models resulting only from the introduction of a different description and a transition region. In addition to this, when the atomistic and the continuum simulations are compared, errors due to the coarse-graining and the discretisation are quantified. In section 2, a concept of the fully coarse-grained atoms, which interact

755

only with the nearest neighbours, was introduced. It can be shown that in the 1D case, the discrete set of equations describing such fully coarse-grained atomistic region is exactly the same as the system obtained by the discretisation of the continuum LLG equation using the second-order finite-difference method with a spatial step ∆x = a. This was used in 1D simulations, section 4.2,

760

where a transition between the real atoms and the continuum region was tested, and where the continuum region was represented by the fully coarse-grained atoms. In the 1D atomistic-continuum examples, the fully coupled time-stepping scheme was used, which conserved the total energy of the system, as neither physical damping (αL = 0) nor numerical damping (gD = 0) was present.

765

In atomistic-continuum modelling of mechanical systems, one of the standard tests involves placing gradients (constant or non-constant) over the atomisticcontinuum interface. In these examples, all energy-based atomistic-continuum methods for mechanics perform either equally (advanced methods for the case of a constant gradient) or worse than just the continuum model discretised

770

with ∆x = a, since the atomistic-continuum interface breaks the symmetry of 43

the interaction. Therefore, it is important to establish the magnitude of the increase of the atomistic-continuum error with respect to the continuum model error (or coarse-graining error) for the case of ∆x = a). The main purpose of the 1D examples of this article is to estimate this and thereby try to establish 775

the limitations of the method. In section 4.3, a 2D example of a domain wall moving in a material with a void is shown. The region around the void was modelled as an atomistic region, while the surrounding area was modelled as a continuum. Finally, in section 4.4, a similar 2D geometry, however without a void, is considered, where the initial

780

state of the magnetisation had a small-scale inhomogeneity in the atomistic region, which corresponded to a local perturbation of the magnetisation due to an external excitation. In 2D continuum, the mesh size was 4 times larger than the interatomic spacing; therefore, a damping band at the atomistic-continuum interface was used. Moreover, in these 2D examples, a non-refined mesh was

785

used as well as the decoupled time-stepping algorithm described in section 2.8. Therefore, the total energy of the system in these examples is not preserved. The purpose of these examples is to test the behaviour of the transition atoms and of the damping band rigorously in a more realistic scenario. 4.1. Material parameters Since all quantities are considered to be dimensionless, βL = 1 was used. The damping was absent (αL = 0) in order to investigate the numerical errors of the proposed method for the most clear-cut example of a non-dissipative system. The DM interaction was not considered, since the coarse-graining technique and the atomistic-continuum transition for the DM interaction is effectively the same as for the exchange interaction, and, in addition, the DM interaction is a less significant interaction compared to the exchange. The interatomic exchange coefficients were based on the following function [45]:
Jij = f sij a−1 Cf ,

 π f (x) = sin πx − exp (− (x − 1) p) x−3 , 2

44

(61)

790

where sij are distances between atoms i and j, a is the interatomic spacing, p is the parameter controlling the rate of the decay and constant Cf was calculated such that the coarse-grained exchange interaction results in Ae = Ae (~ex~ex + ~ey ~ey ), where Ae = 1. Although the parameters of the interatomic interaction do not reflect an actual material, they still contain two

795

main features: the multiple neighbour interaction (not only nearest) and the oscillating decaying nature of the coefficients. Expression (61) does not introduce a loss of generality, it is entirely possible to replace the interaction coefficients defined by this equation with materials specific ones, obtained from the first principles electronic structure theory [36]. In the real atomistic region, a cer-

800

tain cut-off radius of the interaction, rI , is used, which has to be taken into account during the construction of the transition region, the width of which should be larger than rI . The example, which rigorously tests the atomistic-continuum transition, is the domain wall movement, where the magnetisation gradient can be easily controlled by the material parameters. Here the anisotropy tensor K a = Ky ~ey ~ey + Kz ~ez ~ez , where Ky = 2π 2 + 1 and Kz = 1, was selected. The external ~ e = Hx~ex , where Hx = −1. The width of field was applied in x-direction, H q −1 the domain wall can be approximated by wD = π 2Ae (Ky − Kz ) and, with

the parameters defined above, wD = 1. The following initial conditions for the

domain wall were used: m ~ = ~ex sin ν + ~ey cos θ cos ν + ~ez sin θ cos ν, !! √ Hx π 2 (x − x0 ) π ν = arcsin , θ = arcsin tanh + , Ky wD 2

(62)

where x0 is the position of the centre of the domain wall. In the numerical examples, the material parameters as well as the numerical 805

parameters were varied. The material was characterised by two parameters: the ratio of the domain wall width to the interatomic spacing14 , wD /a, and the decay rate of the interatomic interaction, p. The macroscopic exchange was 14 As

for real materials, they can be categorised according to the number of atoms covered

45

Figure 4: The geometries of two 1D regions with the initial position of the domain wall, which were used in the simulations. The position of the atomistic-continuum interface is shown for the cases when the computational region was split into two subregions. The initial magnetisation of the regions is given by equation (62), which results in primarily ~ey and −~ey orientations of magnetic moments within the regions [0, 1] and [2, xL ], respectively.

fixed (Ae = 1) as well as the domain wall width (wD = 1), while Jij were changing with the variation of a and p according to equation (61). 810

4.2. Domain wall movement in 1D At first, domain wall movement in the fully atomistic 1D systems was simulated to calculate a reference error due to time stepping. This error was established using a variation of the time step for models with p = 0.5, rI /a = 8, and wD /a ∈ {4, 16, 64} (results are not shown in figures). Here the error is defined as

815

the L1 norm of the difference between solutions (components of all m ~ i ) divided by 3N , where N is a number of atoms. Time step ∆t = 10−2 was selected such that the error is below 10−3 , which corresponds to an average angular difference of 0.1◦ between the magnetisation distributions provided by two solutions. For problems that are discussed in the subsections below, one-dimensional

820

regions with lengths of xL = 6 and xL = 10 were used; initial position of the domain wall was x0 = 1.5; for the atomistic-continuum simulations, the position of the interface between the continuum and the atomistic descriptions was xI = 3, see figure 4. Neumann boundary conditions were used at x = 0 and x = xL . by a domain wall, e.g., for Ni it is between 800 and 900, for Fe – between 400 and 500, and for Gd – approximately 60, if material parameters are taken from [5] and the exchange constant is calculated for a 2D arrangement of atoms on the (111) fcc crystallographic plane.

46

825

4.2.1. Accuracy of coarse-graining As discussed above, the coarse-graining procedure of the continuum model already introduces a certain error in the solution. Since the assumption of the smoothly-varying magnetisation distribution is used, the error should increase with the increasing gradient of magnetisation. Examples, where fully atomistic

830

systems, which provide reference solutions for error estimates, were compared to fully continuum systems and where wD /a was varied from 4 to 128 (a geometry of length xL = 6 was used, solutions were compared at tr = 10, spatial step of ∆x = a was used, p = 0.5 and rI /a = 8 were used), are shown in figure 5a. It can be seen that the coarse-graining error decreases with the decreasing

835

gradient of magnetisation (i.e. increasing number of atoms per domain wall width). For large gradients, the coarse-graining error makes the simulations unreliable; however, the coarse-graining technique should not be applied in such cases anyway. It is used only when the subsequent increase of the mesh size is possible, to obtain an advantage in computational speed, and for this, a

840

sufficiently smooth magnetisation field is required. The coarse-graining error also depends on the range of the interatomic interaction and should decrease when the range of the interaction decreases. Therefore, in addition to variation of wD /a, parameter p, which characterises the decay rate of the interatomic interaction, was varied. For this comparison, xL = 10,

845

tr = 10, rI /a = 48, wD /a = 16 were used and p was varied from 2−4 (slow decaying interaction) to 21 (fast decaying interaction), see figure 5b. For a fixed magnetisation gradient, the error increases with the decrease of p, which means that for some materials under conditions which lead to fast changing magnetisation, the coarse-graining technique is not applicable due to an unacceptable

850

error. 4.2.2. Introduction of transition region A 1D region of length xL was split into two sub-regions — atoms in the region [0, xI ] were fully coarse-grained atoms (equivalent to 1D continuum discretised at the atomistic scale, i.e. when ∆x = a), followed by rI /a transition atoms, 47

Figure 5: The dependence of the error of the continuum and the atomistic-continuum simulations on the physical parameters: the number of atoms per domain wall (a), wD /a, and the decay rate of the interatomic interaction (b), p, where the fully atomistic simulations were used as a reference. The error is the L1 norm of the difference between the solutions divided by 3N , where N is a number of atoms. The dependence of errors, which are based on the L2 norm divided by 3N and on the L∞ norm, on wD /a is shown in (c). The dependence of the error, which is based on the L∞ norm, on time is shown in (d) for various wD /a.

48

855

followed by real atoms. A domain wall was initiated in the continuum region and due to the external field it moved through the transition region into the real atomistic region. Since for these simulations ∆x = a, they provide the error due to the introduction of the atomistic-continuum transition region, when compared to the reference solution, which, in this case, is the fully atomistic

860

simulation. In figure 5ab, it can be observed that just by introducing an atomistic description in the right part of the computational region and a number of transition atoms adjacent to the interface, without changing the number of degrees of freedom, the error increases compared to the uniform continuum simulation. The

865

error provided by such atomistic-continuum simulation depends on the gradient of magnetisation (expressed here as the number of atoms per domain wall, wD /a) and the decay rate of the interatomic interaction (controlled by parameter p) in the same way as the error of the uniform continuum simulation does. Obviously, in the application of such a multiscale method, the continuum

870

regions are usually much coarser than the atomistic regions, with ∆x being up to an order of magnitude larger than the interatomic spacing. Therefore, in addition to the previous results, errors of the continuum simulations with twice fewer degrees of freedom were calculated (∆x = 2a in these particular examples, since the domain wall should be properly resolved in the continuum). As seen

875

in figure 5ab, the error significantly increases in this case and is larger than the error resulting from the atomistic-continuum system. This is important and shows that when the continuum region is only twice coarser than the atomistic region, the error of the solution is already dominated by the continuum discretisation rather than the atomistic-continuum transition.

880

In this comparison, the atomistic-continuum error is larger than the coarsegraining error (of the continuum model with ∆x = a), although it is of the same order of magnitude with respect to magnetisation gradient. Since the computational region is homogeneous, the atomistic-continuum system suffers from an additional error due to coupling between the subregions. In more

885

elaborate examples, when the computational domain is non-homogeneous or 49

has non-trivial geometry, the atomistic-continuum model should provide a more reliable solution than the continuum model due to limited resolution of the latter. In section 2.5.2, it was discussed that for similar atomistic-continuum method 890

in application to mechanics (QC-QNL), it was shown that the solution error is O(a1+1/p ) in the W 1,p norm [39]. Here, the discrete solutions are compared (3 values per atom/node and N atoms/nodes), therefore, the errors corresponding to the L1 and L2 norms are calculated as the L1 and L2 vector norms of the difference between solutions (components of all m ~ i ) divided by 3N , while

895

the error corresponding to the L∞ norm is calculated as the L∞ vector norm of the difference between solutions. The division is introduced because N is changing with a, as the length of the computational domain is fixed. The reference solution is always the real atomistic solution. Due to the analogy with the mechanical problems that were analysed in [39], it can be expected that

900

the L1 , L2 , and L∞ error norms of the atomistic-continuum system have convergence rates of 2, 1.5, and 1, respectively. However, in figure 5c, it can be seen that the convergence rates for these norms are different, more precisely, these rates are provided in table 1. This can probably be explained by the fact that this example is the dynamic problem — the domain wall (the region of a

905

high gradient) moves from the continuum region into the atomistic region and the transition atoms are subjected to a high gradient only for a limited amount of time, while the most amount of time they are subjected to a small gradient (see equation (62) for the profile of the solution, which is also illustrated in topleft subfigure of figure 6). Errors in figures 5abc correspond to reference time

910

tr = 10, when the domain wall has moved through the interface and is located in the atomistic region. In figure 5d, the dependence of the L∞ norm on time is shown for different a (as wD is fixed). Although it can also be expected that the error temporarily increases at the moment when the domain wall moves through the transition region, which is around t ≈ 5 in figure 5d, this again is

915

not observed. The error of the atomistic-continuum gradually shifts from the error of the continuum system starting from t ≈ 5 almost independently of a. 50

error

L1 norm div. by 3N

L2 norm div. by 3N

L∞ norm

cont.

1.7851

2.4223

1.9608

atom.-cont.

1.8841

2.3912

1.9034

Table 1: The convergence rates for different error norms, which were calculated as slope of log10 (error) vs log10 a, where a−1 ∈ {8, 16, 32, 64, 128}.

Equation (62), which is used as the initial condition for the simulations, provides the precise profile of the moving domain wall (under constant external field) only for the infinite 1D continuum domain15 . However, in the numerical 920

examples, a bounded domain is considered. Moreover, in the case of the atomistic description, which is used as the reference for error calculations, equation (62) is not the precise solution due to the finite coarse-graining error. Simulations in figure 5 were performed with absence of any damping, therefore, as soon as simulation starts, small waves that travel within the system (and reflect

925

from the boundaries) appear. These waves originate from the boundaries and from the fact that equation (62) with fixed ν is not a local minimum of the atomistic energy. Obviously, these non-equilibrium waves are slightly different in the atomistic, the continuum and the atomistic-continuum problems. It is possible that the prefactor of the error due to these non-equilibrium waves is

930

significantly larger than the prefactor of the error due to the domain wall moving through the atomistic-continuum interface. This might explain why the dominant observable term is O(a2 ) and not O(a1+1/p ). To verify this, the analysis of the solution error in the static case is required. Finally, the effect of ghost torques can be quantified based on the interaction

935

coefficients of the examples above. This effect is the local error in the torque given by equation (29), which is acting on the partly coarse-grained atoms, 15 The

profile of the static domain wall can be found by fixing ν and finding the local

minimum of the continuum energy. Afterwards, ν can be found by enforcing equilibrium at θ = 0, π when the external field is present.

51

compared to a representation of the same region with the atomistic model. It can be quantified by the tensorial prefactor given by equation (30), which, in the case of 1D geometry, becomes a scalar. It was observed that the absolute 940

value of this prefactor decreases with the increase of the size of the transition region. In the case when 8 partly coarse-grained atoms were used, PiE was less than 0.02, while for 48 partly coarse-grained atoms, it was less than 0.0043. 4.3. Domain wall movement in 2D The second numerical example is the domain wall motion in a two-dimensional

945

system with a defect. This is a realistic test case, since experimentally domain wall motion in thin films is a frequently investigated topic [46]. The field-induced movement of the domain wall in the presence of a void in the magnetic structure was investigated. This void can correspond to a micro-crack of the sample or to impurity atoms that are produced during the thin-film growth of such magnetic

950

film. Domain wall movement was simulated in the region that is partitioned into the atomistic and the continuum subregions. The (111) plane of the fcc lattice was considered. The region contained a hexagonal void with a side of 6a. The void was surrounded by the atomistic description, which was embedded into the

955

continuum, see figure 6. Since the continuum region was discretised using the regular finite-difference mesh, which was not refined to the atomistic scale, a thin layer of pad atoms was used (since atoms next to the atomistic-continuum interface are local, they require only nearest neighbours). The width and the height of the 2D region were xL = 6 and yL = 2, re-

960

spectively; the initial position of the domain wall was x0 = 1.5. The following material parameters were used: p = 0.5, rI /a = 8, wD /a = 64. The same time step as in section 4.2 was used, ∆t = 10−2 . Neumann boundary conditions were used at x = 0 and x = xL ; periodic boundary conditions were used at y = 0 and y = yL . The continuum region was discretised using spatial step of ∆x = 4a. In

965

this example, the atomistic region (consisting of the real atoms and the transition atoms), had a size of 14∆x (width and height), which resulted in 590 real 52

Figure 6: The geometry of the 2D computational domain with the initial position of the domain wall. The region around the void and the atomistic-continuum interface is magnified. The initial profile of the domain wall, my component of magnetisation, is shown in the upper left subfigure.

atoms. Within the atomistic region, the behaviour of the transition atoms was modified with the additional numerical damping (see section 2.7 for details). Param970

eters of the damping region, the optimal values of which depend on the width of the damping region and the difference between discretisations, were fitted by simulating spin wave propagation. The detailed procedure is described in [33] for a number of numerical examples. The following values were used here: gD = 625, sA = 3∆x. The width of the damping band was selected to be 16a,

975

as the large width of the damping band ensures that the atomistic solution is not contaminated by wave reflections. It should be noted that 8 layers of the fully coarse-grained atoms and 8 layers of the partly coarse-grained atoms (as interatomic interaction range is 8a) were used for this example. The reason for this is that within the damping band, the solution is already modified (high fre-

980

quencies are eliminated) and there is no point of keeping the atomistic precision within this region. By using larger number of layers of the fully coarse-grained of atoms, computational efficiency is increased, as number of summations (over

53

neighbours) is significantly smaller for the fully coarse-grained atoms than for the real atoms. 985

When the domain wall passes through the atomistic region, its dynamics and movement are influenced by the void. The distribution of the magnetisation is shown in figure 7, where only the central part of the 2D region is shown. In this figure, it can be observed that the width of the domain wall is affected by the void and the domain wall becomes thinner in the proximity of the void. From

990

the vector field plot, it can be observed that the atomic magnetic moments have considerable positive mx component at the centre of the domain wall, which was approximately −0.05 at t = 0 (given by initial conditions). This means that the void affects the magnetic structure of the domain wall. Moreover, the simulations show that the presence of the void also temporarily slows down

995

the domain wall movement, which is a well-known effect [1]. An animation illustrating this effect may be found in the supplementary material [47]. To benchmark the multiscale framework, a second example, where the size of the atomistic region was decreased to 12∆x (which corresponds to 204 real atoms), was calculated. Other parameters were kept the same and the error in

1000

the real atomistic region was calculated. In this case, the error is the L1 norm of the difference between the magnetic moments of the real atoms of two solutions (with the sizes of atomistic regions of 12∆x and 14∆x) divided by 3N , where N is the number of real atoms. The evolution of this error with time is shown in figure 8, where it can be observed that there is a peak when the domain

1005

wall passes around the void; however, the error stays below 10−2 (the average difference between the magnetisation distributions provided by two solutions is below 0.8◦ ). After the domain wall passes through the atomistic region, the error decreases. This shows that the solution error due to the size of the atomistic region does not grow in time and for this example is of the same order

1010

of magnitude as the time-stepping error. As mentioned above, the time step was based on 1D simulations of the domain wall movement such that the error is below 10−3 .

54

Figure 7: The distribution of my component of magnetisation around the void in the computational domain at t = 4.5 (left). The regions around the void (upper right) and the atomistic-continuum interface (lower right) are magnified. The magnetic moments of the individual atoms are shown using the vector field (in the case of the continuum, vectors correspond to the rescaled nodal values), where colour represents mz component of magnetisation. For the animation corresponding to this figure see [47].

55

Figure 8: The dependence of the error of the atomistic-continuum simulations due to the size of the atomistic region on time for the problem of the domain wall movement. The error was calculated as the difference between simulations with different sizes of the atomistic regions consisting of 204 and 590 real atoms, respectively.

4.4. Excitation in the atomistic region An evolution of a small-scale inhomogeneity, which corresponds to the local 1015

perturbation of magnetisation due to an external excitation, was modelled using the proposed multiscale technique. The inhomogeneity was placed inside the real atomistic region and examples with different sizes of the atomistic region were compared. The simulation box was partitioned in a similar way as was used in section

1020

4.3, however without a void. The same atomistic stacking and material parameters were used, however without the anisotropy, Ky = Kz = 0, and the external field, Hx = 0. Time step ∆t = 10−5 was used. Two examples were studied, for which the atomistic region (consisting of the real atoms and the transition atoms), had a size of 12∆x (example 1) and 14∆x (example 2), which resulted

1025

in 295 and 681 real atoms, respectively. The initial conditions corresponding to the inhomogeneity were the following:

56

Figure 9: The vector field plot of the magnetic moments of the individual atoms (the real atomistic region only) at different times t showing the evolution of the inhomogeneity created by the initial conditions. Colour represents mz component of magnetisation.

m ~ i = ~ey cos θi + ~ez sin θi ,   (1 − u /u ) θ , u < u i p p i p θi =  0, otherwise,

(63)

where ui = |~ri − ~rp | is the distance between atom i and the centre of the inhomogeneity, which was taken to be ~rp = (0.5xL − up ) ~ex +(0.5yL − up ) ~ey . Radius

of inhomogeneity was taken to be up = 3.5a and the maximum angle θp = 45◦ . With these parameters, the inhomogeneity fits exactly in a quarter of the real 1030

atomistic region of example 1. The evolution of the magnetisation with time for example 1 is shown in figure 9, where only the magnetic moments of the real atoms are shown. It can be seen that the inhomogeneity creates circular waves, which propagate through the atomistic region. Visually it can be seen that these waves do not reflect from

1035

the surrounding region (the transition atoms with the numerical damping) and propagate as if this real atomistic region was a part of the material represented entirely by the atomistic description. To benchmark the multiscale framework, as is done in section 4.3, examples 1 57

Figure 10: The dependence of the error of the atomistic-continuum simulations due to the size of the atomistic region on time for the problem of the excitation in the atomistic region. The error was calculated as the difference between simulations with different sizes of the atomistic regions consisting of 295 and 681 real atoms, respectively.

and 2 (with different sizes of the atomistic region) were compared. The difference 1040

between the atomistic magnetic moments in the real atomistic region, referred to as an error, was calculated (the L1 norm of the difference between the magnetic moments of the real atoms divided by 3N ). The evolution of this error with time is shown in figure 10, where it can be observed that the initially identical solutions start to deviate; however, the error plateaus at the level of 10−3 (which

1045

is also the level of the error on which the selection of the time step was based).

5. Conclusions In this article, a new energy-based atomistic-continuum multiscale coupling method for magnetisation dynamics has been proposed. The method is constructed such that the mismatch between the local continuum description and 1050

the non-local atomistic description is minimised. This is done by introducing an intermediate region, which consists of the transition atoms that interact differently with the neighbouring atomistic and continuum regions. It has been suggested to treat the transition atoms as computational points rather than actual atoms, although they are stacked according to the crystallographic lattice

58

1055

of the atomistic region. The transition atoms, which are located next to the continuum region, interact only with the nearest neighbours; however, the ones that are next to the real atoms, have almost the same interaction range as the real atoms. The interaction coefficients between the transition atoms are evaluated during a preprocessing step and it is sufficient to perform this calculation

1060

only once for a given geometry of the atomistic-continuum interface. In the case when damping is absent and an energy-conserving time-stepping method, such as the mid-point method, is used for the coupled system, the total energy of the dynamic atomistic-continuum system is conserved. The sufficient requirement for the energy conservation is the refinement of the continuum mesh

1065

down to the atomistic scale at the interface, i.e. each interface atom is also a boundary node of the continuum. The method can also be used in the case when the continuum mesh is not refined to the atomistic scale. In this case, a set of links between the atoms and the continuum nodes should be constructed, such that the interaction symmetry is enforced (atom i interacts with continuum

1070

node j in the same way as node j interacts with atom i), which is a non-trivial task in general case. It is also possible to use a non-refined mesh and pad atoms for providing boundary conditions for the atomistic region. In this case, an error in the total energy is introduced, although the use of the transition region still decreases the coupling error. A way of decoupling regions and using

1075

independent time steps in the continuum and the atomistic regions has been suggested. In this article, it was shown that the proposed method has O (1) local error in torque with respect to the interatomic spacing, a, if the continuum parameters and the continuum solution are fixed. The constraints, which are applied to

1080

the transition region, ensure that O (1/a) error terms are eliminated. It was also shown that in the 1D case (and spin magnetic moments constrained to a plane), the magnetostatic and the mechanostatic problems are mathematically equivalent. In [38], it was shown that the lowest achievable local error for any energy-based atomistic-continuum coupling method for mechanical systems is

1085

O (1). Thus, the proposed method is asymptotically optimal. 59

The analytical example of an exchange patch test was considered (the static solution for which angles between any equidistant spin magnetic moments along any given direction are constant) and it was shown that the proposed method has O (a) local error in torque for the case of a 3D patch test. In addition to 1090

this, it was shown that there is an intrinsic mismatch between the atomistic and the continuum descriptions — different torques are required to create the same uniform twist in the atomistic and the continuum domains. This boundary mismatch error appears at the edge of the computational domain, where the interaction symmetry is broken, and is also O (a), if the continuum parameters

1095

and the continuum solution are fixed. Thus, in the case of the exchange patch test, the interface error is of the same order of magnitude with respect to a as the atomistic-continuum boundary mismatch error. Also, a way of resolving the boundary mismatch error (of the exchange interaction) has been suggested — a replacement of the classical continuum model for the exchange effective field by

1100

a Cauchy-Born-like model. However, in order to use this model in a multiscale framework, the behaviour of the transition atoms has to be modified as well, as the boundary mismatch error will simply move to the fully coarse-grained atoms, if the proposed method is used. The main advantage of the proposed method is its flexibility, as well as the

1105

best achievable error for the local-non-local (atomistic-continuum) coupling of the descriptions of magnetism. It is relatively easy to implement the transition region in existing spin dynamics/statics codes, as the method only requires the modification of the interatomic interaction coefficients, but not the functions. As discussed above, the method can also be applied for magnetostatic problems

1110

(in this case, only the part regarding the construction of the transition region is relevant). In the case of dynamics, it can be applied when the continuum mesh is refined to the atomistic scale and a coupled time-stepping is used (i.e. both atomistic and continuum regions are solved together with the same time step) or when the mesh is not refined and the decoupled scheme is used. In the former

1115

case, the method has a useful property of energy conservation (if damping is absent and the mid-point time-stepping method [41] is used). 60

In the case of dynamic simulations, an additional numerical damping has to be applied at the atomistic-continuum interface to avoid reflections of highfrequency waves, if there is a difference between the atomistic lattice spacing 1120

and the spatial discretisation step of the continuum region. In this article, it was proposed to add a damping term to the transition atoms, although the damping band can be wider, as the error due to reflections reduces with the increase of the width of the damping layer. The solution is attenuated towards an average solution and thus a band of damped atoms acts as a low-pass filter for the waves.

1125

A more extended discussion of such damping applied to the multiscale coupling of regions described by the continuum LLG equation can be found in [33]. The major drawback of this approach is certainly the upkeep of a relatively large number of transition atoms during the dynamic calculations, since they should surround the real atomistic region. Moreover, there are a few minor

1130

difficulties, such as necessity to evaluate the interaction coefficients between the transition atoms during a preprocessing step. These coefficients are functions of the geometry of the atomistic-continuum interface, hence, in the case of adaptive multiscaling, where this geometry changes, these coefficients have to be recalculated on the fly. In addition to this, the coefficients of the damping band,

1135

which depend on the relation between the continuum spatial discretisation step and the interatomic distance, have to be calculated/estimated prior the simulation. The precise values of these coefficients (giving the minimum error in the atomistic region) should be evaluated based on separate spin wave simulations. A number of numerical examples was constructed to investigate the applica-

1140

bility and the limitations of the method. It was shown that the coarse-graining error, as well as the interface error, decreases with the decreasing gradient of magnetisation and with the increasing decay rate of the interatomic interaction. Moreover, in the atomistic-continuum multiscale systems, the error is dominated by the discretisation error of the continuum region rather than the error

1145

introduced by the transition region. This was shown using the 1D example of a domain wall propagation from the continuum into the atomistic region. Obviously, a more accurate discretisation of the continuum PDE, e.g. using a 61

higher-order method, will reduce the error in the continuum region, hence is recommended for the multiscale systems. Using a number of 2D examples, it 1150

was shown, that there is almost no dependency of the solution on the size of the atomistic region (the error due to the size can be neglected), if the atomistic region is sufficiently large. Since neither the dimensionality of the position vectors of atoms nor the type of the crystallographic lattice is explicitly used in the derivation of the interaction coefficients within the transition region, this

1155

formalism can also be used in 3D directly with an arbitrary geometry of the atomistic-continuum interface.

Appendix A. Coarse-graining interatomic interactions The coarse-graining of the interatomic interactions is performed by comparing and matching energies of the continuum and the atomistic systems. Al1160

though the procedure of obtaining the continuum coefficients for the case of the exchange interaction is published elsewhere [9, 10], the major steps are repeated here for the reference information and the case of the Dzyaloshinskii-Moriya interaction is added. Appendix A.1. Exchange interaction The total exchange energy of the atomistic system can be written in the following form: E EA = −

1 XX Jij m ~i·m ~ j, 2 i j

|m ~ k | = 1, ∀k.

(A.1)

The effective field of atom i is given by hEA =− i 1165

1 ∂E EA , µ ∂m ~i

(A.2)

where µ is the length of the spin magnetic moments. It should be noted that HiEA = µhEA is used in equation (3), while constant µ is moved into βL and αL , i see equation (1), in order to avoid dragging µ throughout the derivations of the paper. 62

The exchange energy of atom i is εEA =− i

1X Jij m ~i·m ~ j, 2 j

(A.3)

so that the total exchange energy is the sum of the energies of the atoms. Scalar products in (A.3) can be expressed through angles φij between m ~ i and m ~ j . The continuum magnetisation field, m, ~ is constructed such that m ~k=m ~ (~xk ), where ~xk are atomic positions. Using the assumption of a smooth magnetisation field (small φij ), the following transformations are performed: m ~i·m ~ j = cos φij ≈ 1 −

φ2ij , 2

|φij | ≈ |m ~j −m ~ i | = |m ~ (~xj ) − m ~ (~xi )| ≈ |(~xj − ~xi ) · ∇m ~ (~xi )| .

(A.4)

(A.5)

P After adjusting energy (A.3) by a constant ( Jij ), the following expression is obtained: εEA = i

  1X 1X 2 T Jij |~sij · ∇m ~ (~xi )| = Jij (~sij ~sij ) : (∇m ~ (~xi )) · (∇m ~ (~xi )) , 4 j 4 j

(A.6)

where ~sij = ~xj − ~xi . The exchange energy of atom i is related to the continuum energy density as εEC = εEA i /V0 , where V0 is an elementary volume associated

with atom i. Thus, in the continuum case, the energy density depends on vector field m ~ locally and is written in the following way:   1 T εEC = Ae : (∇m) ~ · (∇m) ~ , 2V0 Ae =

1X Jij ~sij ~sij . 2 j

(A.7)

(A.8)

In the continuum case, the effective field is determined by the functional derivative hEC = − E EC =

V0 δE EC , µ δm ~

Z

Ω

63

εEC dV,

(A.9)

(A.10)

where V0 appears because m ~ is defined as the normalised volume average magnetisation and E EC is the total exchange enegry of the magnetic body. Factor 1/2 in (A.7) disappears when the functional derivative is evaluated, see chapter 8 in [10]. This leads to the exchange part of the continuum effective field: H EC = µhEC = Ae : ∇∇m, ~

(A.11)

which enters equation (5), and where tensor Ae contains the exchange interac1170

tion constants. Appendix A.2. Dzyaloshinskii-Moriya interaction Similar technique as for the exchange interaction can be applied to the DM interaction, for which the total energy is the following: E DA =

1 XX ~ Dij · (m ~i×m ~ j) , 2 i j

|m ~ k | = 1, ∀k.

(A.12)

The effective field of atom i is given by hDA =− i

1 ∂E DA . µ ∂m ~i

(A.13)

For the DM interaction, HiDA = µhDA is also used. The DM energy of atom i i is εDA = i

1X~ Dij · (m ~i×m ~ j) . 2 j

(A.14)

The continuum magnetisation field, m, ~ is constructed and the assumption of a smooth magnetisation field is used, which leads to m ~j ≈m ~ i + ~sij · ∇m ~ (~xi ) .

(A.15)

The following expression for (A.14) is obtained εDA = i

 1 X~ 1X~ Dij ·(m ~ (~xi ) × (~sij · ∇m ~ (~xi ))) = − Dij ~sij : (∇m ~ (~xi ) × m ~ (~xi )) . 2 j 2 j (A.16)

As in the case of the exchange interaction, the DM energy of atom i is related to the continuum energy density: εDC =

1 DA 1 1 ε =− D T : (∇m ~ × m) ~ =− m ~ · (∇ · (D e × m)) ~ , V0 i 2V0 e 2V0 64

(A.17)

De =

X

~ ij , ~sij D

(A.18)

j

where two equivalent forms are used (the equivalence can be shown by switching from the tensor notation to the coordinate notation). In the continuum case, the effective field is determined by the functional derivative hDC = −

E DC =

V0 δE DC , µ δm ~

Z

εDC dV.

(A.19)

(A.20)

Ω

Factor 1/2 in (A.17) also disappears when the functional derivative is evaluated. This leads to the DM part of the continuum effective field: H DC = µhDC = ∇ · (D e × m) ~ ,

(A.21)

which enters equation (5), and where tensor D e contains the DM interaction constants.

Appendix B. Coordinate representation of the effective field In the case of the continuum approach, the effective field can be rewritten using the coordinate representation: Hk =

3 X 3 X p=1 q=1

Apq

3 3 3 3 ∂ 2 mk X X X ∂mu X + Dpq ξkqu + Kkp mp +µHke , ∂xq ∂xp p=1 q=1 u=1 ∂xp p=1

Apq =

1X Jij sijp sijq , 2 j

Dpq =

X

sijp Dijq ,

k = 1, 2, 3, (B.1)

(B.2)

(B.3)

j

1175

~ m, ~ e, where Hk , mk , Hke , sijp and Dijq are the components of vectors H, ~ H ~ ij , respectively; Kkp are the components of tensor K a ; ξkqu is the ~sij and D Levi-Civita symbol. 65

Appendix C. Cauchy-Born exchange interaction The Cauchy-Born continuum model for the exchange interaction is derived from the principle that when equation (45) is fulfilled (the state of the uniformly varying magnetisation), the atomistic and the continuum models should have exactly the same values of the exchange parts of the effective fields, which are given by the first terms of equations (3) and (5), respectively. A continuum magnetisation field, m, ~ that coincides with the atomistic spin magnetic moments, m ~ k , at the lattice positions is considered, i.e. m ~ (~xk ) = m ~ k , where ~xk are atomic positions. For given atom i and its neighbour j, vector ~sij is defined as ~sij = ~xj − ~xi . Additional vectors are defined as p~n = ~xi +

n sij ~rij , N

n ≤ N,

~rij =

~sij , sij

sij = |~sij | .

(C.1)

In this case, the angle between m ~ (~ pn ) and m ~ (~ pn+1 ) is the same for all possible n, by equation (45). Moreover, all m ~ (~ pn ) lie in the same plane. The continuum limit can be obtained by N → ∞. This results in angle θ, which is defined as θ = arccos (m ~ (~xi ) · m ~ (~ p)) ,

p~ = ~xi + z~rij ,

z ∈ [0, sij ] ,

(C.2)

and which is a linear function of z. All m ~ (~ p) belong to the same plane by equation (45). Therefore, m ~ (~ p) = m ~ (~xi ) cos θ +

~v (~xi ) sin θ, |~v (~xi )|

~v = ~rij · ∇m. ~

(C.3)

Here the fact that m ~ ⊥ ~v was used (due to |m| ~ = 1). Vector ~v is the directional derivative of m. ~ Therefore, dm ~ dθ |~v (~xi )| = = . dz z=0 dz z=0

(C.4)

Since θ (z) is linear and θ (0) = 0,

θ (z) = z |~v (~xi )| .

(C.5)

This gives m ~j =m ~ (~xi ) cos θij +

~v (~xi ) sin θij , |~v (~xi )| 66

θij = sij |~v (~xi )| .

(C.6)

The exchange part of the effective field (see equation (3)) is vector multiplied by m ~ i in the LLG equation (1), therefore, any components of m ~ j that are parallel to m ~ i can be excluded. Thus, the exchange effective field is the following: ~i = H

X

Jij m ~j =

j

X j

Jij

X ~v (~xi ) sin |~sij · ∇m ~ (~xi )| sin |sij ~v (~xi )| = Jij ~sij ·∇m ~ (~xi ) . |~v (~xi )| |~ s · ∇ m ~ (~xi )| ij j (C.7)

~ i became local, therefore, it can be used It can be seen that the expression for H 1180

for the continuum model.

Acknowledgements The authors would like to thank the reviewers of this paper for very detailed comments, which helped to improve the paper significantly. The authors would like to acknowledge the support of eSSENCE. The authors would also like to 1185

acknowledge the support from the Swedish Research Council (VR) and the KAW foundation (grants 2013.0020 and 2012.0031). The authors are grateful to dr. D. Thonig for drawing authors’ attention to ref. [31] of this paper.

References [1] C. C. H. Lo, A review of the Barkhausen effect and its applications to 1190

nondestructive testing, Materials Evaluation 62 (7) (2004) 741–748. [2] G. Scheunert, O. Heinonen, R. Hardeman, A. Lapicki, M. Gubbins, R. M. Bowman, A review of high magnetic moment thin films for microscale and nanotechnology applications, Applied Physics Reviews 3 (1) (2016) 011301. doi:10.1063/1.4941311.

1195

[3] V. P. Antropov, M. I. Katsnelson, B. N. Harmon, M. van Schilfgaarde, D. Kusnezov, Spin dynamics in magnets: Equation of motion and finite temperature effects, Physical Review B 54 (2) (1996) 1019–1035. doi:10.1103/PhysRevB.54.1019.

67

[4] B. Skubic, J. Hellsvik, L. Nordstr¨ om, O. Eriksson, A method for atomistic 1200

spin dynamics simulations: Implementation and examples, Journal of Physics: Condensed Matter 20 (31) (2008) 315203. doi:10.1088/0953-8984/20/31/315203. [5] R. F. L. Evans, W. J. Fan, P. Chureemart, T. A. Ostler, M. O. A. Ellis, R. W. Chantrell, Atomistic spin model simulations of magnetic

1205

nanomaterials, Journal of Physics: Condensed Matter 26 (10) (2014) 103202. doi:10.1088/0953-8984/26/10/103202. [6] C. Etz, L. Bergqvist, A. Bergman, A. Taroni, O. Eriksson, Atomistic spin dynamics and surface magnons, Journal of Physics: Condensed Matter 27 (24) (2015) 243202. doi:10.1088/0953-8984/27/24/243202.

1210

[7] O. Eriksson, A. Bergman, L. Bergqvist, J. Hellsvik, Atomistic spin dynamics: Foundations and applications, Oxford University Press, 2017. [8] A. Bergman, B. Skubic, J. Hellsvik, L. Nordstr¨ om, A. Delin, O. Eriksson, Ultrafast switching in a synthetic antiferromagnetic magnetic random-access memory device, Physical Review B 83 (22) (2011) 224429.

1215

doi:10.1103/PhysRevB.83.224429. [9] W. F. Brown, Micromagnetics, Interscience Publishers, 1963. [10] A. Aharoni, Introduction to the theory of ferromagnetism, Oxford University Press, 1996. [11] R. E. Miller, E. B. Tadmor, A unified framework and performance

1220

benchmark of fourteen multiscale atomistic/continuum coupling methods, Modelling and Simulation in Materials Science and Engineering 17 (5) (2009) 053001. doi:10.1088/0965-0393/17/5/053001. [12] S. Kohlhoff, P. Gumbsch, H. F. Fischmeister, Crack-propagation in b.c.c. crystals studied with a combined finite-element and atomistic model,

1225

Philosophical Magazine A: Physics of Condensed Matter Structure Defects and Mechanical Properties 64 (4) (1991) 851–878. 68

[13] E. B. Tadmor, M. Ortiz, R. Phillips, Quasicontinuum analysis of defects in solids, Philosophical Magazine A: Physics of Condensed Matter Structure Defects and Mechanical Properties 73 (6) (1996) 1529–1563. 1230

[14] V. B. Shenoy, R. Miller, E. B. Tadmor, D. Rodney, R. Phillips, M. Ortiz, An adaptive finite element approach to atomic-scale mechanics–the quasicontinuum method, Journal of the Mechanics and Physics of Solids 47 (3) (1999) 611–642. doi:10.1016/S0022-5096(98)00051-9. [15] L. E. Shilkrot, R. E. Miller, W. A. Curtin, Multiscale plasticity modeling:

1235

Coupled atomistics and discrete dislocation mechanics, Journal of the Mechanics and Physics of Solids 52 (4) (2004) 755–787. doi:10.1016/j.jmps.2003.09.023. [16] E. B. Tadmor, R. E. Miller, Modeling materials, Cambridge University Press, 2011.

1240

[17] T. Shimokawa, J. J. Mortensen, J. Schiøtz, K. W. Jacobsen, Matching conditions in the quasicontinuum method: Removal of the error introduced at the interface between the coarse-grained and fully atomistic region, Physical Review B 69 (21) (2004) 214104. doi:10.1103/PhysRevB.69.214104.

1245

[18] W. E, J. F. Lu, J. Z. Yang, Uniform accuracy of the quasicontinuum method, Physical Review B 74 (21) (2006) 214115. doi:10.1103/PhysRevB.74.214115. [19] C. Ortner, L. Zhang, Construction and sharp consistency estimates for atomistic/continuum coupling methods with general interfaces: A

1250

two-dimensional model problem, Siam Journal On Numerical Analysis 50 (6) (2012) 2940–2965. doi:10.1137/110851791. [20] C. Ortner, L. Zhang, Energy-based atomistic-to-continuum coupling without ghost forces, Computer Methods In Applied Mechanics and Engineering 279 (2014) 29–45. doi:10.1016/j.cma.2014.06.019. 69

1255

[21] A. V. Shapeev, Consistent energy-based atomistic/continuum coupling for two-body potentials in one and two dimensions, Multiscale Modeling & Simulation 9 (3) (2011) 905–932. doi:10.1137/100792421. [22] M. Luskin, C. Ortner, Atomistic-to-continuum coupling, Acta Numerica 22 (2013) 397–508. doi:10.1017/S0962492913000068.

1260

[23] L. M. Dupuy, E. B. Tadmor, R. E. Miller, R. Phillips, Finite-temperature quasicontinuum: Molecular dynamics without all the atoms, Physical Review Letters 95 (6) (2005) 060202. doi:10.1103/PhysRevLett.95.060202. [24] S. Qu, V. Shastry, W. A. Curtin, R. E. Miller, A finite-temperature

1265

dynamic coupled atomistic/discrete dislocation method, Modelling and Simulation in Materials Science and Engineering 13 (7) (2005) 1101–1118. doi:10.1088/0965-0393/13/7/007. [25] B. Shiari, R. E. Miller, D. D. Klug, Multiscale simulation of material removal processes at the nanoscale, Journal of the Mechanics and Physics

1270

of Solids 55 (11) (2007) 2384–2405. doi:10.1016/j.jmps.2007.03.018. [26] F. Garcia-Sanchez, O. Chubykalo-Fesenko, O. Mryasov, R. W. Chantrell, K. Y. Guslienko, Exchange spring structures and coercivity reduction in FePt/FeRh bilayers: A comparison of multiscale and micromagnetic calculations, Applied Physics Letters 87 (12) (2005) 122501.

1275

doi:10.1063/1.2051789. [27] N. Kazantseva, D. Hinzke, U. Nowak, R. W. Chantrell, U. Atxitia, O. Chubykalo-Fesenko, Towards multiscale modeling of magnetic materials: Simulations of FePt, Physical Review B 77 (18) (2008) 184428. doi:10.1103/PhysRevB.77.184428.

1280

[28] U. Atxitia, D. Hinzke, O. Chubykalo-Fesenko, U. Nowak, H. Kachkachi, O. N. Mryasov, R. F. Evans, R. W. Chantrell, Multiscale modeling of magnetic materials: Temperature dependence of the exchange stiffness, 70

Physical Review B 82 (13) (2010) 134440. doi:10.1103/PhysRevB.82.134440. 1285

[29] T. Jourdan, A. Marty, F. Lan¸con, Multiscale method for Heisenberg spin simulations, Physical Review B 77 (22) (2008) 224428. doi:10.1103/PhysRevB.77.224428. [30] C. Andreas, A. K´akay, R. Hertel, Multiscale and multimodel simulation of Bloch-point dynamics, Physical Review B 89 (13) (2014) 134403.

1290

doi:10.1103/PhysRevB.89.134403. [31] A. De Lucia, B. Kr¨ uger, O. A. Tretiakov, M. Kl¨ aui, Multiscale model approach for magnetization dynamics simulations, Physical Review B 94 (18) (2016) 184415. doi:10.1103/PhysRevB.94.184415. [32] D. Arjmand, S. Engblom, G. Kreiss, Temporal upscaling in

1295

micromagnetism via heterogeneous multiscale methodsArXiv:1603.04920. [33] M. Poluektov, O. Eriksson, G. Kreiss, Scale transitions in magnetisation dynamics, Communications in Computational Physics 20 (4) (2016) 969–988. doi:10.4208/cicp.120615.090516a. [34] L. Bergqvist, A. Taroni, A. Bergman, C. Etz, O. Eriksson, Atomistic spin

1300

dynamics of low-dimensional magnets, Physical Review B 87 (14) (2013) 144401. doi:10.1103/PhysRevB.87.144401. [35] I. Cimr´ak, A survey on the numerics and computations for the Landau-Lifshitz equation of micromagnetism, Archives of Computational Methods in Engineering 15 (3) (2008) 277–309.

1305

doi:10.1007/s11831-008-9021-2. [36] A. I. Liechtenstein, M. I. Katsnelson, V. P. Antropov, V. A. Gubanov, Local spin-density functional-approach to the theory of exchange interactions in ferromagnetic metals and alloys, Journal of Magnetism and Magnetic Materials 67 (1) (1987) 65–74.

1310

doi:10.1016/0304-8853(87)90721-9. 71

[37] W. E, P. B. Ming, Cauchy-Born rule and the stability of crystalline solids: Static problems, Archive for Rational Mechanics and Analysis 183 (2) (2007) 241–297. doi:10.1007/s00205-006-0031-7. [38] M. Dobson, There is no pointwise consistent quasicontinuum energy, IMA 1315

Journal of Numerical Analysis 34 (4) (2014) 1541–1553. doi:10.1093/imanum/drt051. [39] M. Dobson, M. Luskin, An optimal order error analysis of the one-dimensional quasicontinuum approximation, Siam Journal On Numerical Analysis 47 (4) (2009) 2455–2475. doi:10.1137/08073723X.

1320

[40] O. Chubykalo-Fesenko, U. Nowak, R. W. Chantrell, D. Garanin, Dynamic approach for micromagnetics close to the Curie temperature, Physical Review B 74 (9) (2006) 094436. doi:10.1103/PhysRevB.74.094436. [41] M. d’Aquino, C. Serpico, G. Miano, Geometrical integration of Landau-Lifshitz-Gilbert equation based on the mid-point rule, Journal of

1325

Computational Physics 209 (2) (2005) 730–753. doi:10.1016/j.jcp.2005.04.001. [42] S. Bartels, A. Prohl, Convergence of an implicit finite element method for the Landau-Lifshitz-Gilbert equation, Siam Journal On Numerical Analysis 44 (4) (2006) 1405–1419. doi:10.1137/050631070.

1330

[43] L. Debreu, E. Blayo, Two-way embedding algorithms: A review, Ocean Dynamics 58 (5-6) (2008) 415–428. doi:10.1007/s10236-008-0150-9. [44] C. H. S´equin, K. Lee, J. Yen, Fair, G2 - and C2 -continuous circle splines for the interpolation of sparse data points, Computer-Aided Design 37 (2) (2005) 201–211. doi:10.1016/j.cad.2004.06.009.

1335

[45] I. Turek, J. Kudrnovsk´ y, V. Drchal, P. Bruno, Exchange interactions, spin waves, and transition temperatures in itinerant magnets, Philosophical Magazine 86 (12) (2006) 1713–1752. doi:10.1080/14786430500504048.

72

[46] B. Heinrich, Z. Celinski, J. F. Cochran, A. S. Arrott, K. Myrtle, Magnetic anisotropies in single and multilayered structures, Journal of Applied 1340

Physics 70 (10) (1991) 5769–5774. doi:10.1063/1.350156. [47] M. Poluektov, Supplementary video material.

73

Highlights: 

A new energy-based atomistic-continuum partitioned-domain multiscale method for magnetisation dynamics is proposed.



The error due to the coupling of the local continuum description and the non-local atomistic description is minimised, so that the method is asymptotically optimal.



A way of constructing of an intermediate region consisting of QNL-like transition atoms is presented.



A damping band at the atomistic-continuum interface, which absorbs high-frequency waves that are otherwise reflected from the interface, is introduced.