A multi-scale approach for magnetisation dynamics: Unraveling exotic magnetic states of matter

Crystallographic lattice defects strongly influence dynamical properties of magnetic materials at both microscopic and macroscopic length scales. A multi-scale approach to magnetisation dynamics, which is presented in this paper, accurately captures such effects. The method is illustrated using examples of systems with localized, non-trivial topological properties, e.g. in the form of skyrmions and chiral domain walls that interact with lattice dislocations. Technical aspects of the methodology involve multi-scale magnetisation dynamics that connects atomistic and continuum descriptions. The technique is capable of solving the Landau-Lifshitz-Gilbert equations efficiently in two regions of a magnetic material --- the mesoscopic and the atomistic regions, which are coupled in a seamless way. It is demonstrated that this methodology allows simulating realistically-sized magnetic skyrmions interacting with material defects and novel physical effects, uncovered using this theoretical methodology, are described.


INTRODUCTION
It is of paramount importance in material science to simulate the behaviour of matter both at small and large length scales. New and exciting macroscopic phenomena in materials are usually driven by the perturbation of physical properties at local and very small length scales. For example, dislocations or defects of atomic positions influence plastic deformation of materials and in magnetic media, they are a source of domain wall pinning and the Barkhausen effect. In many cases, local properties cannot be disentangled from the global properties of the material, in which they are embedded, and, consequently, both spatial scales are required to be treated on equal footing. This is a highly non-trivial task, which this article addresses. Technological applications, e.g. the promising area of spintronics [1] and the evolving field of quantum computing [2], bring relevance to the field of magnetism. Scientific advancements in this area are provided by several new experimental techniques, e.g. pump-probe measurements of magnetism [3]. For interpretation of experiments, simulations are highly relevant and allows understanding or even predicting new and interesting phenomena. In essentially all fields of physics, simulation tools are used in parallel to experimental observations.
Magnetic phenomena can be studied at a variety of scales, from macroscopic to atomic. Different models are traditionally used at each level in order to consider specific physical phenomena that are significant for a given time and length scale. Two commonly used paradigms are the micromagnetic model and the atomic description that is based on the Heisenberg Hamiltonian. Micromagnetic theory [4] describes phenomena typically at micrometre length scale. This methodology is well established and has been used for several decades. It is known to provide results of a relatively high degree of agreement with experimental data, at least for the length scale and time range it is designed for. Fundamental to this model is the assumption that the magnetisation varies smoothly in space and micro-structural information is embedded into effective parameters. These conditions impose a limit on the applicability of this model to small magnetic structures or intricate atomistic arrangements. On the other hand, the Heisenberg spin model provides a discrete, atomic-scale description of magnetism. This model is capable of representing magnetisation dynamics on an atomic length scale, naturally including effects of lattice irregularities [5]. However, the computational cost of atomistic spin dynamics simulations makes it unfeasible for large computational domains. Therefore, it has difficulties describing the properties of samples at sizes that are used in some experimental measurements.
Micromagnetic models allow distributing computational points in a controlled manner. These methods have been reported in the literature and are based on finite differences, e.g. in Ref. [6,7], or finite elements, e.g. in Ref. [8]. The atomistic model does not permit redistributing computational points, but it can be used in combination with micromagnetic models leading to various multi-scale schemes. There are several versions of such coupling schemes developed in application to magnetic materials [9][10][11][12] and the reader is referred to a review article published in Ref. [13] for an exhaustive overview.
The aim of this paper is to demonstrate the applicability of a multi-scale method to analysing complex magnetisation dynamics in systems where large and small length scales play role, e.g. skyrmion dynamics and its influence from defects on the atomic scale. Limitations of this methodology are established by considering a set of examples that are challenging for any similar numerical method that aims at simulating magnetisation dynamics of materials. To perform the investigation, an implementation of a multi-scale method has been made in the Up-pASD software [14,15], taking on-board techniques from Refs. [16,17]. The micromagnetic model is shown here to have sufficient accuracy for regions where the magnetisation varies smoothly, while atomic-level simulations can reproduce the much more complex behaviour at shorter length scales. Framing a small atomistic region with a thick micromagnetic layer in a way proposed in this paper, permits estimating how perturbations propagate in the material, while avoiding reflections that pollute the region of interest and keeping the computational effort and accuracy mainly within the atomistic part. As is demonstrated below, the multi-scale method outlined here is able to simulate dynamics of topological excitations, e.g. skyrmions with realistic (experimental) size, and their interaction with atomistic defects, such as dislocations.

METHODS
The type of atomistic-continuum coupling considered in this paper is the partitioned-domain method. In this scheme, the computational domain is split into regions that belong to different models, with an explicit interface between the regions. There are two conceptually distinct atomistic-continuum coupling methodologies -the energy-based and the force-based coupling [18]. There are two major challenges for partitioneddomain atomistic-continuum coupling methods: handling non-local atomistic interactions, which is relevant for the energy-based methods applied in both statics and dynamics, and averaging high-frequency excitations within the atomistic solution that must not reflect from the boundary between the atomistic and the coarselydiscretised continuum regions. It must be emphasised that the latter is not an issue of coupling methods and results from disparate discretisations of the regions. Atomistic-continuum coupling methods originate from the field of mechanics and date back to 1990s and, therefore, the highlighted issues were largely addressed in application to modelling deformation, e.g. Ref. [19] solves the local/non-local coupling issue in the energy-based methods in mechanics and Ref. [20] proposes a way of dealing with high-frequency atomistic excitations. In the field of magnetism, the local/non-local coupling has been addressed in Ref. [16], while a technique for averaging high-frequency spin motion at the interface has been proposed in Ref. [17].
The technical aspects of the method proposed and used here involves the merging of an atomistic and a continuum description of the magnetisation dynamics into an established software -UppASD [15]. Most of the technical aspects used in this work have been introduced previously [5,14,16,17] and, therefore, relevant technical aspects of this work are presented in Supplementary Note S1.

RESULTS
It is well-known that it is not possible to construct a completely error-free atomistic-continuum coupling, as discussed in Ref. [18]. Previous studies [16,17] focused on analysing the errors introduced by the coupling approach both analytically and computationally. The aim of the examples of this paper is to show that it is possible to construct the coupling such that atomistic and continuum regions behave in a similar way and the errors introduced by the atomistic-continuum interface are sufficiently small and do not influence the microscopic and the macroscopic behaviours. In the general case, the atomistic-continuum coupling errors depend on the continuum discretisation and on the relation of magnetisation gradient to atomistic lattice spacing [16].
In the computational examples and figures, Cartesian coordinates x, y, z correspond to the horizontal, the vertical and the out-of-plane axes, respectively. In figures, the red and the blue colours indicate magnetisation pointing in-plane and out-of-plane, respectively.

A. Domain wall motion
A basic domain wall motion, with a Hamiltonian that involves Heisenberg exchange and uniaxial anisotropy, is considered. The simulations consider two regions of opposing magnetisation, each having orientation parallel to the lowest energy direction, specified by the easy axis [21]. The precise material parameters for the simulations of domain wall motion are detailed in Supplementary Note S3a). In this example, the domain wall width is approximately 50 interatomic distances. The simulation shown in Fig. 1a) shows that the domain wall, which starts in the continuum region, is not modified or distorted when it reaches the atomistic region -white/yellow box in the centre of Fig. 1a). This illustrates that the handshake between the two regions is seamless and that an object, which is not expected to change its shape, correctly maintains its original configuration when moving across interfaces between the regions. The same conclusion is reached when considering the domain wall that is leaving the atomistic region and entering the continuum region. The dynamics of the domain wall is found in Supplementary Video 1.

B. Domain wall motion and the Barkhausen effect
The cause of the Barkhausen effect has been discussed previously in Ref. [22]. Materials usually contain irregularities in the crystallographic lattice, such as dislocations, stacking faults, impurity atoms, vacancies and lattice distortions due to micro-stresses. A domain wall typ- Figure 1. a) Domain wall crossing atomistic region, which is indicated by the white/yellow square. b)-f) Snapshots at successive steps of the time evolution of the domain wall crossing an atomistic region containing a round hole, which is indicated by the black colour. The atomistic region is indicated by the black frame.
ically needs a local increase in energy around the defect to cross one of these irregularities and, therefore, becomes pinned -parts of the domain wall near the defect lag behind the rest of the domain wall. The unpinning of the domain wall causes a detectable change in the magnetic flux, which was picked up in Barkhausen's coils.
The numerical examples investigated here consider defects in the form of microscopically small, regular shaped holes in the hexagonal crystal lattice. Snapshots for different times of the domain wall motion are shown in Fig. 1b)-f) and are seen to represent the wall moving over the defect. Regions of the domain wall that come close to the hole move quickly towards the hole and get stuck around it, as can be seen in Fig. 1c). The domain wall is stuck for an appreciable time, even if the wall in the rest of the material keeps moving, as can be seen in Fig. 1e). Eventually, in just a relatively short period of time (≈ 200−500 ps) a substantial amount of the material changes its magnetisation. In an experiment, such fast changes of the magnetisation would lead to big changes in the magnetic flux that is connected to the Barkhausen effect (see Supplementary Video 2). Although micromag-netic models have successfully discussed the Barkhausen effect [23] and generate results similar to those of Fig. 1, the approach adopted here has greater descriptive possibilities due to treatment of the defect region in an atomistic way, with unique properties of each atom around the defect and with parameters of the spin Hamiltonian that are evaluated from the first principles theory.

C. Topological magnetic states
In condensed matter physics, magnetic skyrmions are spatially localized excitations that preserve their shape. These skyrmions are homotopically distinct from the ferromagnetic state, meaning that there exists no continuous transformation that transforms the magnetic texture of a skyrmion to the ferromagnetic state. This topological argument has been discussed critically in the literature, particularly with respect to stability and the possibility of using these objects as carriers of information [24,25]. Skyrmions have been observed experimentally, both in static and dynamic conditions [26].
Magnetic skyrmions come in two versions -hedgehog skyrmions have their moments pointing in the radial direction, while vortex skyrmions have their moments circulating around the centre. A particular type of skyrmion depends on the direction of D ij vector with respect to vector r ij connecting atomic positions i and j. The material parameters used in the subsequent simulations, with the exception of skyrmion-dislocation interactions, are listed in Supplementary Note S3b). These parameters correspond to a one-atom thick layer of iron resting atop a substrate of iridium. It has already been demonstrated experimentally that skyrmions exit in this system [27].

Spin spiral and Skyrmion lattice
Considering only the Heisenberg exchange, energy is lowered for an atom and its neighbours when their atomic moments become parallel. In contrast, the energy associated to the Dzyaloshinskii-Moriya interaction can be minimized when the magnetic moments of neighbouring atoms are non-collinear and the minimum energy is found when moments form an angle of 90 degrees. In materials where both terms contribute significantly, a competition between interactions may introduce many different low-energy states. One family of states of these magnets is known as the spin spiral. For these systems, locally, the variation of spins becomes relatively small in one direction and relatively large in an orthogonal direction. There is a tendency for directions of minimal variation to align, creating long lines of aligned magnetisation. In this example, the multi-scale approach is used to describe this class of materials.
The time evolution of the magnetic moments achieving the spin spiral state is shown in Fig. 2a)-c) and in Supplementary Video 3. Initially, the system is thermal- ized by using an annealing process -the temperature is decreased gradually starting from 100 K until it reaches 0.1 K. After the thermalisation process, the system obtained a ferromagnetic state with all moments pointing along the z-direction (data not shown). Afterwards, the magnetic moments are allowed to evolve in time at T = 0.1 K, until the system converges into a spin spiral state, as shown in Fig. 2a)-c). The stripes of opposing magnetisation are easy to identify by the z-component of the moments. If an external field of a sufficient magnitude is applied along the z-axis to the spin-spiral state (5 T in the example here), skyrmions may appear in the material, forming what is known as a skyrmion lattice, see Supplementary Video 4 and Fig. 2d)-f). In simulations of both the skyrmion lattice and the spin spiral configuration, the border between the micromagnetic and atomistic region is completely transparent. This implies that the interface between two regions does not introduce any spurious errors and the time evolution of the magnetisation behaves according to expectations in both regions -the size and the shape of skyrmions is the same in both regions.
In the above simulations, the introduction of micromagnetic regions brings an improvement of computational efficiency, compared to purely atomistic simulations. To characterise the coarseness of the system, an additional variable c is introduced, which represents the number of atoms that fit between the continuum mesh nodes. Note that this parameter only specifies the coarseness of the simulation nodes of the continuum region. Hence, the value c = 1 specifies the case when the atomistic and continuum regions have equal density of simulation atoms and simulation nodes, a case which represents essentially the computational effort of a purely atomistic simulation. In Fig. 3, the simulation time (cpu clock time measured in seconds after the simulation reaches 500 time steps) is shown for a 2D ferromagnet with rectangular shape, in which both Heisenberg exchange and DM interaction are significant, as a function of the number of simulation elements (atoms/nodes), for different values of c, i.e. for different coarseness. The material parameters of this simulation are given in Supplementary Note S3c. The geometry considered in the simulations had rectangular shaped regions of continuum and atomstic regions. In the simulation considered here, the initial magnetisation points along the z-direction and since there is no external magnetic field, the magnetisation evolves in time towards a spin-spiral state. The simulation time scales linearly with the size of the system for any choice of coarseness, as shown in Fig. 3. It may also be seen that the simulation time decreases significantly with increasing values of c. Naturally this comes at a cost, in terms of accuracy, and a suitable choice of c must always be chosen with care, depending on shape of the simulation cell and the dynamical responce one is interested in. A closer inspection of the computational efficiency of the multiscale approach (data not shown) shows that the efficiency scales exponentially with c, clearly demonstrating the major advantage of the multi-scale modeling approach suggested here, compared to a stright atomistic description.

Single-skyrmion creation via microwave fields
In this section, the creation process of skyrmions via an external magnetic field is illustrated, in a magnetic medium in which both the Heisenberg exchange and the DM interaction are important. Simulation details and materials parameters are given in Supplementary Note S3d). Initially, the magnetic moments point along the z-direction and are forced to align parallel by an exter- nal field of 6 T. Single skyrmions are then stabilized by a local torque. Details are given in Supplementary Note S5. The simulated time evolution of these skyrmions is shown in Fig. 4. The results from the micromagnetic model are shown in the topmost region, while the atomistic approach, on a hexagonal crystal lattice, is presented in the bottom region. The external field is applied to the entire material throughout the simulations, as it is found to keep the majority of magnetisation aligned along the initial direction, which bounds the size of the skyrmion [26].
The skyrmions shown in Fig. 4 represent a metastable state. In absence of the external field, the system can lower its energy by expanding the radius of the skyrmion, while due to the external field, atomic moments favour alignment along the z-direction. The competition between the interactions results in a skyrmion with a size following a damped oscillatory mode. This behaviour is known as the so-called "breathing mode". In Fig. 4 the inner magnetic texture of the skyrmions is shown, both for the atomistic and the micromagnetic regions. Due to the selected direction of the DM vector, the topological excitations in both regions are vortex skyrmions as shown in Fig. 4e) and it can be seen that they have similar shape and size [28]. The breathing mode is responsible for the alternating radius of the skyrmions depicted for various time steps in the simulations as shown in Fig. 4b)-d).
Since the resolution of the pulsed local field is constrained by the resolution of the underlying domain, slight differences in its effect are expected and are revealed at the first few steps of the simulation. This illustrates that a continuum description, even when tuned to have materials parameters selected to be as close as possible to represent the behaviour of an atomistic region, does not capture the correct time evolution of the more accurate atomistic model. However, as the magnetic configuration evolves towards a stable stationary state, skyrmions of the continuum and atomistic descriptions become increasingly similar as shown in Fig. 4d). The time evolution of both skyrmions is shown in Supplementary Video 5.

Skyrmion-micro-stress interaction
In this example, a smaller region of different magnetic anisotropy is considered, where locally residual microstresses in a material causes an enhanced anisotropy constant [29]. The simulation shown in Fig. 5 is set up similarly to that of Ref. [30]; however, the triangular region with unique anisotropy (coloured differently in the figure) is treated by an atomistic model on a hexagonal lattice. A hedgehog skyrmion is introduced (created as described in the previous section) and driven along the track by means of a spin current induced by STT. The material parameters for this example are given in Supplementary Note S3e). The direction of the spin current used in this simulation is given by vector (0.995, 0.1, 0). Small y-component is used to compensate the Magnus force that drives the skyrmion downwards as already reported in Ref. [31], resulting in motion along a straight line. A weak spin-current of 4.5 m s −1 was found to be unable to push the skyrmion past the region of enhanced anisotropy. In this case, the skyrmion moves towards the irregularity, where it dissipates energy and angular momentum, which leads to the reduction in size, as shown in Fig. 5a)-c) and Supplementary Video 6. A strong spincurrent induced by the STT can push the skyrmion past the region of enhanced anisotropy -a value of 6 m s −1 is used for the simulation in Fig. 5d)-f) and Supplementary Video 7 and was found to be strong enough to push the skyrmion past the micro-stress region. An intermediate spin-current of 5 m s −1 sustains the skyrmion as a stationary entity, even when it meets the magnetic irregularity (data not shown).

Skyrmion-dislocation interaction
An edge dislocation is a common crystallographic defect where one plane of atoms is missing a half, which distorts nearby planes of atoms. Such defects are unavoidable in experimental samples, where they may influence the magnetic properties in a decisive way and in this subsection, the dislocation-affected skyrmion motion is illustrated. Results of this section cannot be obtained by a continuum model alone and illustrate the strength of the atomistic-continuum multi-scale approach. The process of determining the positions of atoms within the dislocation is described in the Supplementary Note S6. Material parameters are given in Supplementary Note S3f).
The sequence of images in Fig. 6 shows snapshots of the skyrmion when it moves across the computational domain. Fig. 6a)-f), Supplementary Video 8 and Fig. 6g)-l), Supplementary Video 9 correspond to the STT vectors with components (1.25, 0.175, 0) and (5, 0.5, 0) m s −1 respectively. In the latter case (i.e. relatively strong STT), the skyrmion interacts with the defect, however, it overcomes the region around the defect and subsequently leaves the atomistic region, moving to the continuum. For a weaker STT current, the skyrmion moves at the top of the row of atoms of the dislocation. Once it comes to the end of the missing row of atoms, which is indicated by the red lines in Fig. 6, it moves to the bottom of the dislocation and changes the direction -it moves backwards, against the direction of the STT. After reaching for the second time the left hand side of the dislocation defect, the direction the skyrmion reverses again. Such oscillatory behaviour of the skyrmion around the dislocation is repeated several times. Due to damping, the skyrmion travels a shorter path each iteration, until it is finally trapped by the defect as shown in Fig. 6f). This counter-intuitive behaviour, where for a certain time the skyrmion moves against the direction favoured by the STT, has not been reported in literature before and is impossible to observe within a continuum model.
In Fig. 7 and Supplementary Video 10, the time evolution of x and y components of the magnetisation of the skyrmion is shown in the neighbourhood of the dislocation. The colours indicate the direction of the major component of the magnetisation in the xy-plane -red, yellow, green and blue colours correspond to the directions of the magnetisation primarily along right, down, left and up directions, respectively. The same colour pattern is used to describe the topology of the magnetic in-teractions at the endpoints of the dislocation. Since the DM vector is in-plane and is perpendicular to the bond direction, the dislocation on the left endpoint introduces a local distortion of the DM interaction in such a way that it favours magnetisation pointing opposite to x-direction. For the right endpoint of the dislocation, the situation is reversed. As a consequence, the ground-state magnetic configuration is achieved when the atomic magnetic moments on the left and the right endpoints are tilted to the left and to the right, respectively, as shown in green and red colors in Fig. 7a). The white background colour indicates that the magnetisation is pointing along the zdirection.
As depicted in Fig. 7, the skyrmion stays on the upper side of the defect when moving to the right and in the lower side when moving to the left. Initially, the moments on the left and the right sides of the skyrmion point to the left and to the right, respectively, as can be seen in Fig. 7a). However, interaction with the dislocation introduces an asymmetry in the components of the magnetisation. When the skyrmion is approaching the dislocation, Fig. 7b), the red part of the skyrmion is facing the green endpoint of the dislocation. This creates a repulsive potential with two local magnetisations directions that are antiparallel and to lower the exchange interaction, these regions try to avoid each other. Since the STT current is pushing the skyrmion to the right in the figure, this external torque overcomes the repulsive potential and the skyrmion is moved along the upper side of the dislocation. The DM interaction coupled to the lattice defect causes a chirality and in the example shown in Fig. 7 favours a clockwise rotation of the magnetisation. This pushes the skyrmion towards the upper part of the dislocation when it is moving to the right. In Fig. 7d) the skyrmion is located at a position where it is bounded by ferromagnetic exchange interaction at both sides of the dislocation, causing a local minimum in the energy landscape. However, since the skyrmion has a finite linear momentum produced by the STT, it overcomes this bound state. As the skyrmion continues to move to the right along the dislocation it reaches the endpoint and here it stays for a short time on the right side of the dislocation. At this point the STT balance out the intrinsic exchange and DM interaction between the skyrmion and the lattice dislocation, thus, the right part of the skyrmion is pinned at the right side of the dislocation in another bound state, Fig. 7e). The inertia of the skyrmion together with the clockwise rotation of the magnetisation induced by the DM interaction now favours the skyrmion to move along the left direction under the dislocation. Notice also that the intensity of the spin current induced by the STT is not strong enough to overcome the inertia of the skyrmion during this short time. Figure 7f) shows another possible bound state, similar to the one in Fig. 7d), which again is overcome by the inertia of the skyrmion. Finally, in Fig. 7g) the skyrmion is pinned or trapped by the attractive potential between the green region of the dislocation and the green part of the skyrmion. On this way back to the final state, the skyrmion looses energy and linear momentum until it is trapped in the left side of the dislocation and, since the STT current is not sufficient to take the skyrmion out of this bounding potential, the skyrmion is pinned on the left side of the dislocation.

CONCLUDING REMARKS
In this communication, an efficient computational methodology, which enables a multi-scale description of magnetisation dynamics, and software are presented. The method proposed here is of the partitioned-domain type and is implemented in the UppASD software. Using a number of examples, it is demonstrated that the continuum and the atomistic descriptions can be interfaced in a seamless way. As shown here, this enables simulations of magnetic phenomena at atomistic length scales, which are coupled to a micromagnetic description and which allow simulating cells that are of the same length scale as experimental samples.
As illustration of the multi-scale method, several examples of chiral and topological magnetic states are used -domain walls and skyrmions that interact with local perturbations of atomistic length scale. With this technique, an atomistic description of the Barkhausen effect is obtained. In addition, skyrmion dynamics is investigated. Of the more conspicuous results presented here, one may note the intricate dynamics of a skyrmion interacting with a lattice defect -an edge dislocation. The counter-intuitive motion, which may go against an external STT, is analysed in detail and is discussed in terms of the magnetic interactions coupled to the topology of the defect and the topology of the skyrmion.
The continuum regions have certain discretisation freedom compared to the atomistic region, which leads to the major advantage of the multi-scale approach. However, as dynamics of objects, such as skyrmions and domain walls, is modelled, the mesh density of the continuum region should allow accurate resolution of the modelled objects. This would also mean that sometimes, if magneti-sation gradient is relatively large (i.e. fields change relatively rapidly), the continuum mesh must be relatively fine, which alleviates the computational advantages of using the continuum regions. One of the possible solutions in this case, is using the dynamic remeshing and/or dynamically moving atomistic/continuum regions. This, however, might require additional resources and increases the complexity of the software.
The method proposed here allows simulating magnetic phenomena in general, e.g. in the field of magnonics or for studies of racetrack memories, with topological or chiral objects. Simulations of such phenomena calls for realistically-sized simulation cells, which only a micromagnetic simulation model allows for. This level of description must be combined with an atomistic description, to treat regions where unavoidable lattice defects influence the magnetisation dynamics, and the method outlined here allows for such a description opening up for simulations of many exciting magnetic phenomena and technologies.

I. DATA AVAILABILITY
The source code of the developed tool is freely available upon request to the authors (http://katalog.uu.se/empinfo/?id=N12-216).

II. ACKNOWLEDGEMENTS
We are grateful for the support from STandUpp, eSSENCE, the Swedish Research Council, the Knut and Alice Wallenberg Foundation, the Foundation for Strategic Research and the Swedish Energy Agency. We acknowledge D. Thonig for fruitful discussions.
In the atomistic approach, the dynamic behaviour of magnetic moments of individual atoms is described by the atomistic Landau-Lifshitz-Gilbert equation [1,2]: where γ is the gyromagnetic ratio, λ is the phenomenological Gilbert damping constant, m i is the direction of spin magnetic moment with |m i | = 1, µ is the length of spin magnetic moment and h i is the effective field. The spintransfer torque (STT) is represented by the vector τ i while b 1 , b 2 and the components of the vector g are STT coefficients characteristic of the material. The vector s ik is connecting atoms i and k (where x k are atomic positions) and the summation is over atoms in a neighborhood of atom i, such that k = i. The following form of the effective field is considered: where J ij are constants of the Heisenberg exchange interaction between atoms i and j, vectors D ij describe the Dzyaloshinskii-Moriya (DM) interaction, K a is the anisotropy tensor and h e is the external field. Summations are taken for all j, such that j = i, where any number of neighbours (i.e. non-nearest) can be included. At the continuum scale, the magnetisation dynamics is modelled by the continuum version of the Landau-Lifshitz-Gilbert (LLG) equation [3,4]: where m is the normalised magnetisation field (|m| = 1) and β L and α L constants have the same meaning as in Eq. (1). In the case of the continuum approach, the following effective field is considered: where tensor A e contains the exchange interaction constants while the tensor D e includes the Dzyaloshinskii-Moriya interaction constants. The anisotropy tensor is represented by K a and h e is the external field. After a quick inspection arXiv:1910.07807v1 [cond-mat.mtrl-sci] 17 Oct 2019 of Eqs. (1) and (5), it is easy to see that both equations retain a similar mathematical form. It is straightforward to prove that both equations are invariant under the following transformations: where s ij 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. The tensors A e and D e and the vector d are spatially homogeneous, i.e. these tensor/vector parameters do not depend on i.
The usual way of comparing the atomistic and the continuum models is by considering the continuum solution that coincides with the atomistic solution at all lattice points [5]. To find the difference between the models, the continuum solution is fixed and the asymptotic behaviour of the difference is found as a function of the atomistic lattice spacing. In the case of magnetism, a continuum magnetisation field m is considered in such a way that it coincides with the atomistic spin magnetic moments m k at lattice positions, i.e. m (x k ) = m k , where x k are atomic positions [6].
Most atomistic lattices of crystalline materials can be categorised as point-symmetric atomistic lattices, i.e. for each pair of atoms i and j, there exists an atom k which is of the same type as atom j, such that s ij = −s ik . The derivation of Eqs. (8) and (9) can be found in Ref. [6], where it is shown that for the point-symmetric atomistic lattices, the atomistic and the continuum models are consistent with the error estimate as a → 0, where a is the atomistic lattice spacing and the effective fields are indicated in Eqs. (4) and (7). The proof of the consistency between the atomistic an the continuum spin-transfer torque terms, Eqs. (2) and (6), respectively, can easily be derived by Taylor-expanding solution m (x i ) = m i , as was performed in Ref. [6]. This results in the following asymptotic behaviour: In multiscale partitioned-domain atomistic-continuum coupling approach, the entire computational region is split into two subregions -the atomistic and the continuum domains. There is a "sharp" atomistic-continuum interface (i.e. a curve in two-dimensional systems, a surface in three-dimensional systems) between these two regions, as illustrated in Fig. 1.
The continuum region is discretised using the finite-difference method. The atomistic-continuum interface encapsulates a layer of finite-difference nodes and is linear between the neighbouring nodes. The coupling of the regions is implemented by constructing padding atoms and padding nodes, which provide boundary conditions for the atomistic and the continuum regions, respectively. Real atoms interact with padding atoms, while the solution at padding atoms is obtained from the continuum region. Finite-difference nodes interact with padding nodes, while the solution at padding nodes is obtained from the atomistic region.
The solution at padding atoms is obtained by bilinear interpolation of the continuum solution with normalisation. The normalisation is introduced to preserve the length of spin magnetic moments. The solution at padding nodes is obtained by normalised weighted average of the atomistic solution inside the box with side ∆x centred at the padding node, where ∆x is the inter-node distance. For all atoms inside the box, the weight is assigned as the area of the intersection of the box with side a centred at the atom and the box with side ∆x centred at the padding node. The normalisation is introduced to preserve the nodal length of the vector field solution.
To reduce the ghost-torque error for systems with non-nearest-neighbour atomistic interactions, the behaviour of atoms close to the atomistic-continuum interface is modified [6]. More specifically, J ij and D ij are modified for a finite set of atoms. Two types of modified atoms can be distinguished -the fully coarse-grained atoms, which interact only with the nearest neighbours, and the partly coarse-grained atoms, which also interact with non-nearest neighbours. The modification of J ij and D ij is performed such that the fully coarse-grained atoms completely isolate the atomistic-continuum interface. More details about the exact way of modifying the atoms can be found in Ref. [6]. To reduce the high-frequency wave-reflection from the atomistic-continuum interface, additional numerical damping is added to atoms close to the atomistic-continuum interface [6,7]. This damping acts as a low-pass filter for the waves travelling from the atomistic region to the continuum, as the solution is "attenuated" to an average solution within a certain window. Due to dispersive nature of the spin waves, the damping is non-linear and depends on time derivative of the solution. The analysis of the dynamics of the damping layer and the exact form of the modification can be found in [7].
The time stepping was performed simultaneously for all degrees of freedom, i.e. all finite-difference nodes and atoms. The solution at padding atoms/nodes at a certain time step is obtained from the continuum/atomistic solution of the same time step. Mid-point [8] and Depondt [10] methods were used to solve Eqs. (1) and (5)  There exist some configurations of Heisenberg or Dzyaloshinskii-Moriya interactions that cannot be represented by the presented micromagnetism model. There are two cases known to the authors. The first example occurs in magnetic materials with large angle between adjacent magnetic atoms, e.g. in antiferromagnets or in non-collinear materials. Examples of such materials are bcc Cr, fcc Fe and most elemental forms of Mn. The problem of these materials is that spins tend to vary abruptly in space. These variations lead to larger errors in differentiation via finite difference. When the moment on consequent nodes is very close to be antiparallel, the result of the differentiation is dominated by numerical error, and the result of such simulation is completely meaningless. Furthermore, the local averages used for padding elements and damping band yield also erroneous results on antiparallel or nearly-antiparallel configurations. Variations of the finite difference method that accept antiferromagnets are available in literature.
The other class of exchange configurations that are problematic are those for which tensor A e or D e result in a null matrix. It is mathematically possible for this to happen on the Heisenberg exchange, but to the best knowledge of the authors there is no physical material that can exhibit such configuration. A known physical example of this happening for DM interaction is found in the pyrochlore antiferromagnet [9]. The planes perpendicular to the [111] direction in this material present a Kagome structure with a DM configuration such that if r ij = −r ik then D ij = D ik . This can be achieved in Kagome structures without breaking the antisymmetry of the DM interaction because there is no translational symmetry on the exchange vectors. A consequence of this configuration is that r ij D ij = −r ik D ik . Since in the structure for each neighbour to an atom there exist a neighbour in the opposite direction, it follows that i =j Having a micromagnetics domain with 0 antisymmetric exchange does not represent the effects observable on the atomistic material. As an example, skyrmions and anti-skyrmions can exist on the atomistic representation, but cannot occur on a material without DM interaction.

b) Derivative errors on the interface.
The spatial derivative of the magnetization function is considered in several parts of this model. The micromagnetics effective field terms for exchange indirectly introduce these derivatives via the transformed parameters explained in Supplementary note S4. The spin-transfer torque model directly computes this derivative via finite-differences in the continuum. When the setup is made in such a way that padding nodes fall very closely or directly over atoms, the error due to interpolation is minimized. As a demonstration, two configurations with ∆x a = 1, STT, Heisenberg and DM are displayed. On one of them, atoms and nodes are aligned, while on the other, there is an offset of a 2 . A skyrmion is introduced via microwave fields and pushed through the boundary. The structure and trajectory of the skyrmion is shown for each setup in Fig. 2. Thus, in Fig. 2a), the non-aligned atoms on the interface are shown above. At the bottom of the figure, the trajectory of the skyrmion shows a clear deformation when crossing the interface. Notice that the blue area on the right represents the continuum region. In Fig. 2b) the aligned atomistic and continuum grids show no deformation on the interface.
When derivatives are computed on nodes next to the interface, the interpolation error from the padding node propagates to the approximated derivative. This error is studied for the 1D setup and shown in Fig. 3. The elements m 2 , m 1 and m 0 are finite difference nodes, while s 2 , s 1 and s 0 represent atoms. Elements m 0 and s 0 are padding elements. The error is studied for the derivative at m 1 . An exact solution for the magnetization at a point f (x) is assumed to exist. Furthermore, the solution at non-padding elements is assumed to be exact, that is The padding node m 0 has its moment determined as which, since the moments in our model are always unit vectors, it can be simplified using the cosine law as where θ is the angle between f (x + h) and f (x + 3h). When f (x) is sufficiently smooth, the angle θ is close to 0. It is known that f (x) needs to be a smooth function to be represented by the continuum, so this assumption does not add a new constrain to the model.
The derivative at m 1 is calculated via finite difference as which can be expressed in terms of f (x) as Applying Taylor expansion on the previous expression results in This shows the error introduced by this approximation is linear with respect to h. The absolute error is generally small. It is hard to measure its effect in domain walls. However, more sensitive configurations, like spin ices or skyrmions, specially when all STT, DM and Heisenberg exchange are present, do show visible artifacts near the interface. A higher order interpolation is not generally possible since there are not necessarily enough atoms surrounding a padding node. One possible solution is to discard padding nodes and calculate a finite difference stencil coefficients for each node next to an atom in the interface. These stencils can be used instead of the suggested in Supplementary note S4 to calculate new exchange parameters for nodes in the interface.  0.0 mRy direction (0, 0, 1) After thermalization Field (0, 0, −0.2)T The uniaxial anisotropy favours magnetism along the z direction. An external magnetic field of −0.25 T along the z-direction, provides a torque that moves the wall. For the domain wall-stress interaction, the stress presents anisotropy K1 = 0.001 mRy within the defect. b) Iron-Iridium monolayer. The parameters for these set of simulations are taken from literature, corresponding to an iron monolayer on iridium [11]. The parameters are written in function of the magnetic moment magnitude, a moment magnitude of 1 µ B is used for simplicity.
For skyrmion-stress another value of DM was used, to tune the size of the skyrmion allowing for larger ∆x a . Parameter transformations derivations.
A procedure is described below to calculate micromagnetics solutions using the existing atomistic solver. The exchange parameters needed for an atomistic simulation to behave as the micromagnetics description are here derived step by step. The finite-difference stencil is defined as in Fig. 4. The goal is to place UppASD simulation elements (considered atoms in the initial design of the tool) in the locations of finite difference nodes, and obtain from the solvers implemented in UppASD software [12] the same calculations a finite difference solver would carry out. To do so only parameters need to be tweaked. To find the relationship between parameters, let us focus on a finite difference node, p. The magnetic moment at the node p is m p . To make notation simpler, labels are assigned to magnetisation vectors at neighbouring space locations. The subindices n, s, e, w, a and b refer to positive and negative y, x and z directions, respectively. Labels are chosen after North, South, East, West, Above and Below.
where A pq are the entries of the stiffness matrix A e , x i is the i th coordinate, and d is the number of dimensions of the space.Non-diagonal elements of A e are zeros and this leaves us with a simpler expression: Using this notation, discrete approximations of second order partial derivatives of m are of the form: Thus, the Heisenberg term in the Hamiltonian must be: The Hamiltonian appears in a cross-product by the local magnetization in the equation of motion. Since m p × m p = 0, the term involving the local magnetization m p can be ignored.
Having written the continuum Hamiltonian exchange term in this form, it is possible to appreciate that it is equivalent to the atomistic Hamiltonian exchange when J ij is defined as: Then, a correspondence between De ij and D j * can be found:

SUPPLEMENTARY NOTE S5
A pulsed field is used to generate locally a strong magnetic field of 600 T in the −z direction. This causes a fliping of spins of atoms in a small region of the simulation box. Such strong fields can be achieved by e.g. a Laser via the inverse Faraday effect [13], and enable the creaton of skyrmions. The pulsed field used to create the skyrmions, has a radius of 10 lattice spacing (27.15Å in total) and lasted 2 ps. After the pulse was applied, moments in the regions affected quickly flip their spin. Neighbouring magnetic moments rotate to reduce the high exchange energy introduced by the abrupt change in magnetization direction, and a skrmion state is stabilized.

SUPPLEMENTARY NOTE S6
Atoms tend to keep a constant distance from each other in the lattice. As an atom is displaced, its neighbours pull it back, trying to preserve the inter-atomic distance. The resulting behaviour is similar to that of a mesh of springs or elastic bands. An added constraint is that the atoms surrounding the dislocation must minimize the potential energy with respect to the rest of neighbouring atoms. This is a requirement since the continuum cannot adapt to the distorted atomistic lattice, and the transition atoms need to behave in the same way as the continuum. Placing the atoms is then formulated as an optimization problem, where the atoms at the boundaries of the domain are locked to lattice points and one row of atoms is removed in the interior. The target function is the elastic potential energy for the configuration of atoms. Assuming the material is isotropic, the potential energy can be written as where k is a material parameter and ∆p i is the displacement of the atom from its rest position. The positions are found using an iterative gradient descent algorithm on the energy while preserving the constrained atoms in place.
Most of the atomic distances are very close to the lattice parameter. Around 97.5% of the bonds are shorter than 1.09a, but the longest bond is of 1.36a. As some atoms were removed, no bond between atoms is shorter than 1a. Varying bond distances have an effect in the interaction strengths. For Heisenberg exchange, a polynomial is fitted to the values of a known material with cubic (bcc) lattice. The material used as reference is iron-bcc. The calculation of Dzyaloshinsky-Moriya interaction vectors is a bit more cumbersome. The Levy-Fert mechanism can be used to model the DM for some materials, but the topology of the surrounding layers is needed [14]. For the sake of demonstration, a much more simplified model was used. The DM vectors chosen are coplanar to the material, parallel to the bond and their magnitude is proportional to the distortion of the length. The magnitude of the DM vectors when the link is not distorted was chosen to be 20% of the Heisenberg exchange for first neighbours. A script was written in order to generate this configuration. The values are calculated for J ij and D ij as J ij =J 0 (−0.05988 r ij 3 + 0.95532 r ij 2 − 4.9516 r ij + 8.31778) (55) where J 0 = 1 mRy is the scale of the exchange, and R is either the identity matrix when the DM vectors lying along the bond connecting atoms, or a 90 degree rotation matrix for DM vectors lying perpendicular to the bond. The evaluation of the polynomial for J ij gives a value of 1 when r ij is exactly one lattice spacing (2.48Å).