Proposal for a micromagnetic standard problem for materials with DzyaloshinskiiMoriya interaction
Abstract
Understanding the role of the DzyaloshinskiiMoriya interaction (DMI) for the formation of helimagnetic order, as well as the emergence of skyrmions in magnetic systems that lack inversion symmetry, has found increasing interest due to the significant potential for novel spin based technologies. Candidate materials to host skyrmions include those belonging to the B20 group such as FeGe, known for stabilising Blochlike skyrmions, interfacial systems such as cobalt multilayers or Pd/Fe bilayers on top of Ir(111), known for stabilising Néellike skyrmions, and, recently, alloys with a crystallographic symmetry where antiskyrmions are stabilised. Micromagnetic simulations have become a standard approach to aid the design and optimisation of spintronic and magnetic nanodevices and are also applied to the modelling of device applications which make use of skyrmions. Several public domain micromagnetic simulation packages such as OOMMF, MuMax3 and Fidimag already offer implementations of different DMI terms. It is therefore highly desirable to propose a socalled micromagnetic standard problem that would allow one to benchmark and test the different software packages in a similar way as is done for ferromagnetic materials without DMI. Here, we provide a sequence of welldefined and increasingly complex computational problems for magnetic materials with DMI. Our test problems include 1D, 2D and 3D domains, spin wave dynamics in the presence of DMI, and validation of the analytical and numerical solutions including uniform magnetisation, edge tilting, spin waves and skyrmion formation. This set of problems can be used by developers and users of new micromagnetic simulation codes for testing and validation and hence establishing scientific credibility.
I Introduction
In computational science socalled standard problems (or benchmark or test problems) denote a class of problems that are defined in order to test the capability of a newly developed software package to produce scientifically trustworthy results. In the field of micromagnetism which, to a significant extent, relies on results produced by computer simulations, the micromagnetic modeling activity group (mag) at the National Institute of Standards and Technology has helped to define and gather a series of such standard problems for ferromagnetic materials on their website. mumag () Those five problems cover static as well as dynamic phenomena and over the years more standard problems have been proposed including the physics of spin transfer torque, Najafi2009 () spin waves Venkat2013 () and ferromagnetic resonance. Baker2017 () Thus far, however, standard problems for materials with DzyaloshinskiiMoriya interaction (DMI) have not been defined in the literature. Originally, the socalled DMI was phenomenologically described by Dzyaloshinskii Dzyaloshinskii1958 (); Dzyaloshinskii1964 () to explain the effect of weak ferromagnetism in antiferromagnets, and later it was theoretically explained by Moriya Moriya1960 () as a spin orbit coupling effect. The DMI effect is observable in magnetic materials with broken inversion symmetry and can be present either in the crystallographic structure of the material Bogdanov1989 (); Leonov2015 () or at the interface of a ferromagnet with a heavy metal. Crepieux1998 (); Fert2013 (); Wiesendanger2016 () In contrast to the favoured parallel alignment of neighbouring spins from the ferromagnetic exchange interaction, the DMI favours the perpendicular alignment of neighbouring spins. The competition between these interactions allow the observation of chiral magnetic configurations such as helices or skyrmions, where spins have a fixed sense of rotation, which is known as chirality.
Skyrmions are localised and topologically nontrivial vortexlike magnetic configurations. Although they were theoretically predicted almost thirty years ago, Bogdanov1989 () only recently have skyrmions started to attract significant attention by the scientific community because of multiple recent experimental observations of skyrmion phases in a variety of materials with different DMI mechanisms. Muehlbauer2009 (); Yu2010 (); Milde2013 (); Romming2013 (); Leonov2016 (); MoreauLuchaire2016 (); Pollard2017 (); Shibata2017 (); Nayak2017 () The magnetic profile of a skyrmion changes according to the kind of DMI present in the material. Wellknown skyrmionic textures are Néel skyrmions, Bloch skyrmions and antiskyrmions. Fert2013 (); Wiesendanger2016 (); Nayak2017 (); Hoffmann2017 () The former two are named according to the domain walllike rotation sense of the spins.
The fixed chirality of spins imposed by the DMI causes skyrmions to have different properties from structures such as magnetic bubbles Nagaosa2013 (); Kiselev2011 (); BernandMantel2017 () or vortices Butenko2009 (); Mruczkiewicz2017 (); Chen2018 (). In addition, the antisymmetric nature of the DMI has an influence on the dynamics of excitations such as spin waves, making them dependent on their propagation direction in the material.
In this paper, we define a set of micromagnetic standard problems for systems with different DMI mechanisms. This set of problems is aimed at verifying the implementation of the DMI by comparing the numerical solutions from different software with semianalytical results from published studies Rohart2013 (); Bogdanov1994 () where possible. In this context, we test these problems using three opensource micromagnetic codes, OOMMF, OOMMFreprt () MuMax3 Vansteenkiste2014 () and Fidimag. Fidimag () In Section III we introduce our analysis by defining the theoretical framework to describe ferromagnetic systems with DMI, which is used to obtain numerical and analytical solutions. Consequently we describe the problems starting by the specification of a one dimensional sample in Section IV, where the DMI has a distinctive influence on the boundary conditions. Then in Section V we test the stabilisation of skyrmionic textures in a disk geometry for different kinds of DMI. In Section VII we compute the spin wave spectrum of an interfacial system in a long stripe and show the antisymmetry produced by the DMI. Finally, in Section VI we analyse a skyrmion in a bulk system, where the propagation of the skyrmion configuration across the thickness of the sample is known to be modulated towards the surfaces because spins acquire an extra radial component. Rybakov2013 ()
Ii The DzyaloshinskiiMoriya interaction
The DzyaloshinskiiMoriya interaction Dzyaloshinskii1958 (); Dzyaloshinskii1964 (); Moriya1960 () (DMI) is a spin orbit coupling effect that arises in crystals with a broken inversion symmetry. In these materials the combination of the exchange and spin orbit interactions between electrons leads to an effective interaction between magnetic moments of the form
(1) 
where the vector depends on the induced orbital moments. In general, the DM vector will be nonzero, although it is strongly constrained via Neumann’s principle, that is, that the Hamiltonian shares (at least) the symmetry of the underlying crystal system. In this context, for two ions it is usually possible to strongly constrain the direction of through symmetry arguments, such as the ones given in Refs. Yosida1996, ; Moriya1960, . For example, if a mirror plane runs perpendicular to the vector separating two ions, passing through its midpoint, and if the separation is parallel to , then this operation sends , , and . This means that the transformation causes components proportional only to to change signs, forcing to vanish but leaving and nonzero. Yosida1996 ()
When dealing with the continuum (micromagnetic) version of the DMI, the same above considerations apply but may be generalized through the use of a phenomenological approach based on Lifshitz invariants (LIs) Landau1980 (); Bogdanov2002 (). Systems featuring LIs range from ChernSimons terms in gauge field theories, Lancaster2014 () to chiral liquid crystals, Sparavigna2009 () but also include magnetic systems hosting DMIs. In this latter context they are said to describe inhomogeneous DMIs Bogdanov1989 (); Bogdanov2002 () owing to the spatial variation of the magnetization that they describe, in the form
(2) 
where . The precise forms of the LIs are dictated by the crystal symmetry of the system and they determine the micromagnetic expression of its DMI energy. In this continuum limit, the energy written in terms of LIs encodes symmetry constraints elegantly using only a single parameter , by including the way in which the magnetisation (or spin) changes along the different spatial directions. These continuum expressions of the DMI energy are equivalent to the discrete version (equation 1).
The DMI phenomenon also occurs at surfaces Crepieux1998 () and in interfacial systems Crepieux1998 (); Bogdanov2001 (); Fert2013 () because of the breaking of symmetries. In the latter case, chiral interactions which lead to LIs in the free energy of the system, could arise from broken symmetries reflecting lattice mismatch, defects or interdiffusion between layers Bogdanov2001 (). A specific example of this Fert2013 () is interfacial DMI arising from the indirect exchange within a triangle composed of two spins and a nonmagnetic atom with strong SO coupling. Fert2013 () From the atomistic description of interfacial DMI it is possible to derive expressions in the continuum based on LIs as shown in Refs. Yang2015, ; Rohart2013, .
In general, for bulk and interfacial systems the application of the LIbased continuum theory for inhomogeneous DMI, leads to the prediction of a rich variety of noncollinear magnetic structures such as vortex configurations. Bogdanov1989 (); Bogdanov1994 (); Bogdanov2001 (); Leonov2015 () An extended discussion on the theory of the DMI and further examples are discussed in Section S1 of the Supplementary Material.
Iii Analytical model of chiral configurations
For the description of chiral structures in magnetic materials with DMI it is customary to write the magnetisation in spherical coordinates that spatially depend on cylindrical coordinates Bogdanov1994 (); Leonov2015 ()
(3) 
where are spherical angles. In the general case, and with being the cylindrical coordinates. In the case of two dimensional systems or magnetic configurations without modulation along the thickness of the sample, and are specified independent of the direction.
For a chiral ferromagnet without inversion symmetry, we are going to describe stable magnetic configurations in confined geometries Rohart2013 () when considering symmetric exchange and DMI interactions, an uniaxial anisotropy and in some cases, an applied field. According to this, the energy of the magnetic system is given by
(4)  
where is the exchange constant, is the uniaxial anisotropy constant, the applied field and the last term is the DMI energy density, which can be written as a sum of Lifshitz invariants Bogdanov1989 (); Bogdanov1994 (); Leonov2015 () (see Section II). For a material with symmetry class or , the DMI energy density is specified as
(5)  
(6) 
For a thin film with interfacial DMI or a crystal with symmetry class , located in the plane, the energy density of the DMI is
(7)  
(8) 
For a crystal with symmetry class , the DMI energy density reads Bogdanov1989 ()
(9)  
(10) 
Axially symmetric magnetic configurations that are uniform along the direction can be found by substituting the magnetisation (equation 3) into equation III, with . Accepted solutions for are obtained according to the structure of the DMI. Bogdanov1989 (); Bogdanov1994 (); Leonov2015 () For the class material, , with . In interfacial systems, the DMI has the structure of a symmetry class material, where with . For the symmetry class an accepted solution is .
Onedimensional wire  

Dimensions  
Magnetic parameters  
Iv Onedimensional case: Edge tilting
In a onedimensional magnetic system in the direction, where , and at zero field, we can simplify the expression for the energy (equation III) using and (bulk), where , or (interfacial), where . We therefore obtain the following differential equation after minimising the energy with a variational approach Rohart2013 ()
(11) 
where and . The positive sign in the boundary condition refers to the interfacial case and the negative sign to the class material. We solve equations 11 using the shooting method. For this, we refer to the alternative condition for at or derived in Ref. Rohart2013, , which is valid when the system has a large anisotropy and reads
(12) 
Depending on the chirality of the system, which can be observed from the simulations, we fix the condition and vary until finding a solution that satisfies . The upper sign refers to the interfacial case and the bottom sign to the bulk DMI case.
In Fig. 1a and b we compare results from the theory and
simulations of the one dimensional problem, for systems with interfacial
() and bulk () DMI, respectively. For every case we use
permalloylike parameters to test the problem, as specified in
Table 1. This material has associated an exchange length of
and a
helical length of . Simulations were
performed with the finite difference OOMMF, Fidimag and MuMax3 software. In our
examples we used a discretisation cell of
volume, whose dimensions are well below the exchange length. The profile of the
component and either the component of the magnetisation, for the case
of interfacial DMI, or the component for the bulk DMI case, specify the
chirality of the magnetic configuration. To obtain the correct chirality in the
simulations, the DMI energy expression must be carefully discretised, as
explained in the Appendix A. For a onedimensional
system, and since we are using common magnetic parameters, the major difference
between the profiles of systems with different type of DMI, is the
orientation of the spin rotation. Therefore, for a crystal with or
symmetry, the profile of the component resemble the profile of
the interfacial DMI case, which is according to the spin rotation favoured in
the and symmetries. Accordingly, we only show the interfacial and
bulk DMI solutions in Fig. 1a and b, respectively. In the
plots, data points from the simulations are compared with the solutions of
equations 11. In general, OOMMF and Fidimag produce
similar results that perfectly agree with the theoretical curves obtained from
the solutions of the differential equations (as shown in
Fig. 1). Specifically, in the interfacial case the average
relative error (between the semianalytical and simulation curves) for the
component is about 3.8% and for is about 0.3%. Equivalent
magnitudes are found for the bulk DMI system. In the case of MuMax3, a similar
agreement is found when imposing periodic boundary conditions along the
direction of the onedimensional system because the DMI calculation is
implemented with Neumann boundary
conditions Vansteenkiste2014 (); Leliaert2018 () rather than free
boundaries
V Twodimensional case: magnetisation profile of a skyrmion
Twodimensional disk  

Radius  
Thickness  
Magnetic parameters  
It has been shown in Refs. Carey2016, ; Beg2015, ; Rohart2013, that in a confined geometry spins at the boundary of the system slightly tilt because of the boundary condition and, due to the confinement, skyrmions can be stabilised at zero magnetic field. Experimentally observed skyrmion configurations in materials with three different types of DMIs have been reported in the literature: (i) Interfacial DMI, which favours Néel spin rotations and is equivalent to the DMI found in systems with crystal symmetry. (ii) The so called bulk DMI, which favours Bloch spin rotations and is found in systems with symmetry class or , such as FeGe. And, recently, (iii) a DMI found in systems with symmetry class where structures known as antiskyrmions can be stabilised Nayak2017 (); Hoffmann2017 () (antiskyrmions have also been found in interfacial systems but they are best described within a discrete spin formalism Hoffmann2017 (); Camosi2017 ()). These three DMI mechanisms can be described by a combination of Lifshitz invariants with a single DMI constant.
We propose a two dimensional cylindrical system of radius and thickness to test the stabilisation of skyrmions using the three aforementioned DMIs, using permalloylike magnetic parameters as in Section IV (see Table 2).
We summarise in Fig. 2 results obtained for three different skyrmion structures stabilised with the three kind of DMIs. These magnetic configurations were simulated with finite difference codes, as shown in Fig. 2a, which shows a good agreement between them, thus we only plot results from Fidimag in Fig. 2b (see Section S3 in the Supplementary Material for a more detailed comparison). The three skyrmions (see Fig. 2c) are energetically and topologically equivalent, and the magnetization profile depends only on the accepted solution for the angle when described in spherical coordinates. Therefore, the out of plane component of the spins must match for the three configurations. Solving this system analytically, we can calculate the angle for the skyrmion solution with the corresponding boundary condition Rohart2013 (), by minimising equation III. We compare the out of plane component of the spins, , with that of the simulations by extracting the data from the spins along the disk diameter, which we show in Fig. 2a. As in the onedimensional case, we observe the characteristic canting of spins at the boundary of the sample.
To distinguish the three different systems, we compute the skyrmion radius by finding the value of where , and plot the radial component of the spins (see Appendix B) located at a distance from the disk centre. Since spins are in plane at , then and the radial component (see Appendix B and equation 3) as a function of is . Therefore, , and . According to this, we see in Fig. 2b the simulated skyrmion radial profiles at for the , and symmetry class materials, which agree with the theory (shown in dashed lines and curves).
In Fig. 2c we illustrate the three different configurations. The radial component of the magnetisation is shown with a colormap and the outofplane component is shown in grayscale, where white means , thus it is possible to distinguish the region that defines the skyrmion radius, which is highlighted in black, and the slight spin canting at the disk boundary.
For this twodimensional problem, the three simulation packages produce similar results and agree well in the solutions, matching the boundary conditions from the theory. The theory predicts a skyrmion radius of . For the calculation of the skyrmion radius in the simulation results, we use a third order spline interpolation of the profile from Fig. 2a. According to this, OOMMF and Fidimag give a radius of and MuMax3 produces a radius of 22.1 nm, which is a slightly better approximation to the theoretical result.
Cylinder  

Radius  
Thickness  
Magnetic parameters  
Vi Threedimensional case: Skyrmion modulation along thickness
Skyrmions hosted in interfacial systems are in general effectively two dimensional structures since these samples are a few monolayers thick. In contrast, in bulk systems skyrmions can form long tubes by propagating their doubletwist modulation along their symmetry axis, Buhrandt2013 (); Milde2013 (); Rybakov2013 (); Carey2016 (); Charilaou2017 () which we will assume is along the sample thickness. Moreover, it has been shown by Rybakov et al. Rybakov2013 () that there is an extra spin modulation along the skyrmion axis that can be approximately described by a linear conical mode solution. This extra modulation is energetically favourable in a range of applied magnetic field and sample thickness values, where the latest is defined below the helix period . This modulation is greatest at the sample surfaces and is not present at the sample centre along the thickness. It can be identified by an extra radial component acquired by the spins, which is also maximal near the region where in every slice normal to (see Fig. 3c).
We define an isolated skyrmion in a FeGe cylinder of diameter and thickness with its long axis (the thickness) in the direction (see the system illustrated in Fig. 4). We simulate the cylinder using finite differences with cells of volume. Relaxing the system with an initial state resembling a Bloch skyrmion, Beg2015 () and an applied magnetic field of , we stabilise a skyrmion tube modulated along the thickness of the sample. A characteristic parameter of FeGe is the helical length .
Results from Fidimag simulations are shown in Fig. 3. We obtain a skyrmion tube along the sample thickness with a radius that slightly varies along the direction (we define as the radius where in a planecut). In a slice of the cylinder at , we compute a skyrmion radius of , and this value decreases towards the top and bottom surfaces down to , which is negligible in the scale of the chosen mesh discretisation. We emphasize that we are not considering the demagnetizing field in this problem, which can enhance this effect.
We analyse the magnetisation field profiles for different slices (different ) by plotting the components of the spins located across a layer diameter (or from the skyrmion centre), , up to the sample boundary, , as shown by the red dashed line in every snapshot of Fig. 3c. To quantify the radial modulation of spins, and because of the axial symmetry of skyrmions, we calculate the cylindrical components and (see Appendix B). Consistent with the results of Ref. Rybakov2013, , Fig. 3a reveals that at the middle of the sample, i.e. at the planecut located at , there is no extra radial component of the spins, which is observed in a two dimensional skyrmion. In addition, due the confined geometry the azimuthal component slightly increases in magnitude towards the sample boundary with an opposite sense of rotation than that of the skyrmion. Towards the sample surfaces, located at , spins obtain an extra radial component that increases linearly with the distance. We illustrate this effect in Fig. 3b for the bottom layer at . Interestingly, the maximum of and are slightly shifted with respect to the position at , or where , which can be seen from the data points around the dashed line of Fig. 3b. This same effect occurs at the top layer, but with the radial component pointing inwards towards the skyrmion centre, thus looks like a mirror image of that of Fig. 3b. Furthermore, we notice that the radial component towards the boundary also changes sign as does. Snapshots with a zoomed view of the sample for the bottom, middle and top layers are shown in Fig. 3c. We show spins at the skyrmion boundary, where , colored according to their radial component, and with the background colored according to the component.
The linear dependence of the radial component as a function of towards the surfaces is shown in Fig. 4, where we plot as a function of at three different positions in every layer: the centre, , close to the skyrmion radius (according to the discretisation of the mesh), , and at the sample boundary, . These spatial positions are shown as dots in the cylinder system at the top of Fig. 4, with lines denoting where the data is being extracted. From the curves of Fig. 4 we notice that the radial increment is maximal close to the skyrmion radius and is slightly smaller, and with opposite orientation, at the cylinder boundary normal to the radial direction.
Our results show that the skyrmion at the slice does not have a radial modulation and the skyrmion size remains nearly constant across the sample thickness. Hence, it would be possible to use a twodimensional model, similar to the one used in Sections III and V, to describe its profile. In performing this comparison (see Section S4 in the Supplementary Material) we noticed that the skyrmion in the cylinder system has a larger skyrmion radius than the model predicts. In Ref. Rybakov2013, an approximate solution is provided as an ansatz for the angle, which is based on a linear dependence on . Although this approximation qualitatively describes the effects observed from the simulations, a more accurate solution would be possible to obtain by taking the general case and , but it generates a nontrivial set of nonlinear equations to be minimised. Because of the consistent skyrmion size across it is likely that the dependence on in the angle only appears as a weak term or a constant, which differentiates the solution from that of the twodimensional model.
Testing this problem using the OOMMF code we obtain equivalent results with the same skyrmion radius size. In the case of MuMax3, we had to use a slightly different mesh discretisation, since the software only accepts an even number of cells, which we adjusted to get a similar sample size. Results from MuMax3 simulations produce a skyrmion with larger radii compared to Fidimag and OOMMF, with magnitudes of approximately 15.9 nm close to and 15.6 nm at the cylinder caps. Nevertheless, the tendencies of the radial profiles of the magnetisation are still close to the ones obtained with the other codes. Details of these simulations are provided in Section S5 of the Supplementary Material. Although a cylinder system is also suitable for finite element code simulations, a cuboid geometry is more natural to a finite difference discretisation. Hence, we performed a similar study using a cuboid with periodic boundary conditions. In general results on this geometry are equivalent to the cylinder but with two main differences: the periodicity removes the effects at the boundaries and the skyrmion is slightly larger in radius. These solutions are shown in Section S6 of the Supplementary Material.
Vii Dynamics: asymmetric spin wave propagation in the presence of DMI
Analyzing the dynamics of a magnetic system is a standard method to obtain information about the magnetic properties of the material, such as damping or the excitation modes of the system, among others. In particular, it is known that the spin wave spectrum of a material with DMI saturated with an external bias field, is antisymmetric along specific directions where spin waves are propagating. CortesOrtuno2013 (); Moon2013 () These directions depend on the nature of the DMI, and from the antisymmetry it is possible to quantify a frequency shift from modes with the same wave vector magnitude but opposite orientation, i.e. from waves travelling in opposite directions. This frequency shift depends linearly on the DMI magnitude of the material and hence it is a straightforward method for measuring this magnetic parameter. This has been proved in multiple experiments based on Brillouin light scattering. Belmeguenai2015 (); Di2015 (); Tacchi2017 () Accordingly, a standard problem based on spin waves offers the possibility to test the DMI influence on the spin dynamics of the system.
Twodimensional stripe  

Dimensions  
Magnetic parameters  
In interfacial systems spin waves propagating perpendicular to a saturating bias field, which are known as DamonEshbach modes, exhibit antisymmetric behaviour. CortesOrtuno2013 (); Moon2013 () To simulate this phenomenon, we refer to the method specified in Ref. Venkat2013, and use values from Table 4,

Define a thin stripe with the long axis in the direction and thickness along the direction.

Saturate and relax the sample using a sufficiently strong bias magnetic field. If the relaxation is done with the LLG equation it is possible to remove the precessional term and use a large damping to accelerate the relaxation.

Excite the system with a weak periodic field, based on a sinc function, in a small region at the centre of the stripe and applied in a specific direction ,
(13) We delay this signal by , and then excite the system during the time , saving the magnetization field every interval of duration .

Using the magnetization field files, we extract the dynamic components of the magnetization for a chain of magnetic moments along the direction across the middle of the stripe. The dynamic component is obtained by subtracting, from the excited spins, the components of the magnetic moments of the relaxed state obtained in step 2.

We save these components in a matrix, where every column is a magnetization component, or , of the spins across the spatial direction, and every row represents a saved time step saved in the previous step.

Perform a twodimensional spatialtemporal Fourier transform of the matrix, applying a Hanning windowing function Kumar2012 ().
We define a permalloy stripe, with magnetic parameters specified in Table 4, saturating the magnetization into the direction. For this sample we take into account dipolar interactions. To obtain a spectrum for positive and negative wave vectors, i.e., for waves propagating in opposite directions, we excite the system in a small region of width at the centre of the stripe, with a weak periodic signal based on the cardinal sine wave function. We excite DamonEshbach spin waves by applying the sinc field in the direction for a duration of ns. According to Table 4 we save steps to generate the spin wave spectrum.
The result of the spin waves simulation using Fidimag, after processing the data, is shown in Fig. 5. In the spectrum we compare the result using the theory of Moon et al. Moon2013 () for systems with interfacial DMI. The asymmetry in the spin wave depends on the DMI sign. To compare the theoretical curve with the data from the simulations we calculated the minimum in the dispersion relation for both curves. For the simulations we calculated the peaks with largest intensity from the spectrum and fit the data with a fourth order polynomial. The theory predicts that the minimum is located at with a frequency of . From the Fidimag simulation we estimate the minimum at and , which shows they are in good agreement. As an extra test we can notice from Ref. CortesOrtuno2013, that in systems with crystallographic classes or , the DamonEshbach spin waves will not be antisymmetric, however they would be for spin waves excited along the field direction.
In Fig. 5 we observe the presence of extra modes with a smaller intensity signal. These modes can be filtered by setting an exponential damping towards the boundaries of the stripe Venkat2018 () to avoid the reflection of spin waves. In the codes provided in the manuscript we implemented functions for the damping with a simple exponential profile that can be used to only obtain the main branch from the spectrum.
Using the OOMMF and MuMax3 codes for simulating spin waves in systems with DMI produced equivalent results to Fig. 5. With OOMMF simulations we obtained a minimum at and , while simulations performed with MuMax3 produce values of and , which is a slightly better approximation. Results for OOMMF and MuMax3 simulations, and details about the numerical interpolation to the curves are shown in Section S7 of the Supplementary Material.
Viii Conclusions
We have proposed four standard problems to validate the implementation of simulations of helimagnetic systems with DMI mechanisms found in crystals with , and symmetry class, where the former is also relevant in interfacial systems. The strength of the three DMI types we use in the problems can be quantified by a single DMI constant. For the onedimensional and twodimensional problems we test the boundary condition in confined geometries, which can be compared with analytical solutions. Moreover, profiles of different skyrmionic textures, which vary according to the DMI kind, are characterised by their radial profile, in particular at a distance from the skyrmion centre where , which we define as the skyrmion radius. Further, in order to test the effect of the DMI on the dynamics of the systems, we propose a problem based on the excitation of spin waves and the calculation of their spectrum. In this case, we analyse DamonEshbach spin waves in a stripe with interfacial DMI (or, equivalently, a crystal with symmetry), which is known for being antisymmetric, and compare the solution with analytical theory. Finally, we analyse an isolated skyrmion in a bulk material with symmetry class in a cylinder. In this sample the skyrmion profile propagates through the thickness and acquires an extra radial modulation. We notice that this modulation is nonexistent in a slice at the middle of the sample along the thickness direction and increases linearly towards the cylinder caps (normal to the direction). Additionally, it is greatest at the skyrmion radius (where ), decreases to zero at the skyrmion centre and towards the skyrmion boundary (in every slice), and is present at the cylinder boundary (normal to the radial direction) with an opposite orientation than the one within the skyrmion configuration.
Simulations in this study have been performed using codes based on the finite difference numerical technique. Since many of the problems are compared with semianalytical calculations the results can be also applied to finite element code simulations. Some finite element computations with our nonpublicly available software Finmag are shown in Section S8 of the Supplementary Material. In addition, we compared our data with the results from a nonpublic finiteelement code developed by R. Hertel, which is an entirely rewritten successor of the TetraMag software Hertel2007a (); Kakay2010 (). These results are also shown in Section S8, where we obtained an excellent quantitative agreement.
With this set of problems we intend to cover the functionality of the DMI interaction implemented in a micromagnetic code by testing boundary conditions, energy minimisation, which can be achieved using LLG dynamics or minimisation algorithms such as the conjugate gradient method, and spin dynamics. Overall, the micromagnetic codes used in our testings significantly agree with expected solutions and comparisons with the theory, thus our results substantiate studies based on micromagnetic simulations with the three codes we have tested. We hope this systematic analysis helps to promote the publication of codes in simulation based studies for their corresponding validation and reproducibility, and serve as a basis for more effective development of new simulation software.
For the realisation of some of the problems, we have implemented new DMI modules for MuMax3 Cortes2018 () and OOMMF OOMMFa (); OOMMFb (); OOMMFc (); Cortes2018 () that take advantage of the computer softwares framework, such as GPU implementation in MuMax3 or the robustness of OOMMF. We have used the Jupyter OOMMF (JOOMMF) interface to drive OOMMF and analyse data. JOOMMF2017 (). Scripts and notebooks to reproduce the problems and data analysis from this paper can be found in Ref. Cortes2018, .
Acknowledgements.
This work was financially supported by the EPSRC Programme grant on Skyrmionics (EP/N032128/1), EPSRCâs Centre for Doctoral Training in Next Generation Computational Modelling, http://ngcm.soton.ac.uk (EP/L015382/1) and EPSRCâs Doctoral Training Centre in Complex System Simulation (EP/G03690X/1), and the OpenDreamKit â Horizon 2020 European Research Infrastructure project (676541). We acknowledge useful discussions with the MuMax3 code team.Appendix A Finite difference discretisation for the DMI
Assuming a two dimensional film positioned in the plane, the energy density for the interfacial DMI used in this study is modeled as
(14)  
(15) 
The corresponding effective field of this interaction reads
(16)  
(17) 
When using finite differences, we can discretise the derivatives using a central difference at every mesh site. Thus, for example, the central difference for the first derivative of with respect to is
(18) 
where is the component of the closest magnetic moment at the mesh site in the direction and is the mesh discretisation in the direction. We can then, for every mesh site, collect all the terms related to the contribution to the field from its 4 mesh neighbours in the plane where the system is defined. For instance, the field contribution from the mesh neighbour is
(19)  
(20) 
The contribution for the other neighbours have the same structure except the denominator for the neighbours in the direction will have a factor of instead of (similar for ). In addition, the cross product is, in general, given by , with the unit vector from the mesh site directed towards the neighbour in the direction. Hence, the calculation for the field is similar than that of the Heisenberglike model, with an equivalent DMI vector of the form .
It is important to mention that an interfacial DMI described within the discrete spin model using the DMI vector , in the continuum limit leads to expression 15. However, when using finite differences for the continuum description of the system, the calculation of the field related to equation 15, which is similar to the atomistic model calculation, has an opposite sign for the equivalent DMI vector.
Similarly, for a class material, the finite differences discretisation leads to a calculation of the micromagnetic DMI field using a vector . In the case of symmetry, this vector is in the directions and in the directions.
Appendix B Cylindrical components
The cylindrical components of the magnetization are computed with a transformation matrix according to
(21) 
where is the azimuthal angle.
See pages 1 of SupplementaryMaterial.pdf See pages 2 of SupplementaryMaterial.pdf See pages 3 of SupplementaryMaterial.pdf See pages 4 of SupplementaryMaterial.pdf See pages 5 of SupplementaryMaterial.pdf See pages 6 of SupplementaryMaterial.pdf See pages 7 of SupplementaryMaterial.pdf See pages 8 of SupplementaryMaterial.pdf See pages 9 of SupplementaryMaterial.pdf See pages 10 of SupplementaryMaterial.pdf See pages 11 of SupplementaryMaterial.pdf See pages 12 of SupplementaryMaterial.pdf See pages 13 of SupplementaryMaterial.pdf See pages 14 of SupplementaryMaterial.pdf See pages 15 of SupplementaryMaterial.pdf See pages 16 of SupplementaryMaterial.pdf See pages 17 of SupplementaryMaterial.pdf See pages 18 of SupplementaryMaterial.pdf See pages 19 of SupplementaryMaterial.pdf
Footnotes
 This argument was discussed through an internal communication with the MuMax3 team. Free boundaries are going to be available in future MuMax3 releases.
References
 http://www.ctcms.nist.gov/ rdm/mumag.org.html.
 Najafi, M. et al. Proposal for a standard problem for micromagnetic simulations including spintransfer torque. J. Appl. Phys. 105, 113914 (2009).
 Venkat, G. et al. Proposal for a Standard Micromagnetic Problem: Spin Wave Dispersion in a Magnonic Waveguide. IEEE Trans. Magn. 49, 524–529 (2013).
 Baker, A. et al. Proposal of a micromagnetic standard problem for ferromagnetic resonance simulations. J. Magn. Magn. Mater. 421, 428 – 439 (2017).
 Dzyaloshinskii, I. A thermodynamic theory of “weak” ferromagnetism of antiferromagnetics. J. Phys. Chem. Solids 4, 241–255 (1958).
 Dzyaloshinskii, I. E. Theory of helicoidal structures in antiferromagnets. II. Metals. Soviet Physics JETP 19, 960–971 (1964).
 Moriya, T. Anisotropic superexchange interaction and weak ferromagnetism. Phys. Rev. 120, 91–98 (1960).
 Bogdanov, A. & Yablonskii, D. Thermodynamically stable ”vortices“ in magnetically ordered crystals. The mixed state of magnets. Zh. Eksp. Teor. Fiz 95, 178 (1989).
 Leonov, A. O. et al. The properties of isolated chiral skyrmions in thin magnetic films. New J. Phys. 18, 1–12 (2015).
 Crépieux, A. & Lacroix, C. DzyaloshinskyMoriya interactions induced by symmetry breaking at a surface. J. Magn. Magn. Mater. 182, 341–349 (1998).
 Fert, A., Cros, V. & Sampaio, J. Skyrmions on the track. Nat. Nanotechnol. 8, 152–156 (2013).
 Wiesendanger, R. Nanoscale magnetic skyrmions in metallic films and multilayers: a new twist for spintronics. Nature Reviews Materials 1, 16044 (2016).
 Mühlbauer, S. et al. Skyrmion lattice in a chiral magnet. Science 323, 915–919 (2009).
 Yu, X. Z. et al. Realspace observation of a twodimensional skyrmion crystal. Nature 465, 901–904 (2010).
 Milde, P. et al. Unwinding of a Skyrmion Lattice by Magnetic Monopoles. Science 340, 1076–1080 (2013).
 Romming, N. et al. Writing and Deleting Single Magnetic Skyrmions. Science 341, 636–639 (2013).
 Leonov, A. O. et al. Chiral Surface Twists and Skyrmion Stability in Nanolayers of Cubic Helimagnets. Phys. Rev. Lett. 117, 087202 (2016).
 MoreauLuchaire, C. et al. Additive interfacial chiral interaction in multilayers for stabilization of small individual skyrmions at room temperature. Nat. Nanotechnol. 11, 444–448 (2016).
 Pollard, S. D. et al. Observation of stable Néel skyrmions in cobalt/palladium multilayers with Lorentz transmission electron microscopy. Nat. Commun. 8 (2017).
 Shibata, K. et al. Temperature and Magnetic Field Dependence of the Internal and Lattice Structures of Skyrmions by OffAxis Electron Holography. Phys. Rev. Lett. 118, 087202 (2017).
 Nayak, A. K. et al. Magnetic antiskyrmions above room temperature in tetragonal Heusler materials. Nature 548, 561–566 (2017).
 Hoffmann, M. et al. Antiskyrmions stabilized at interfaces by anisotropic DzyaloshinskiiMoriya interactions. Nat. Commun. 8, 308 (2017).
 Nagaosa, N. & Tokura, Y. Topological properties and dynamics of magnetic skyrmions. Nat. Nanotechnol. 8, 899–911 (2013).
 Kiselev, N. S., Bogdanov, A. N., Schäfer, R. & Rößler, U. K. Chiral skyrmions in thin magnetic films: new objects for magnetic storage technologies? J. Phys. D: Appl. Phys. 44, 392001 (2011).
 BernandMantel, A. et al. The skyrmionbubble transition in a ferromagnetic thin film. ArXiv eprints (2017). eprint 1712.03154.
 Butenko, A. B., Leonov, A. A., Bogdanov, A. N. & Rößler, U. K. Theory of vortex states in magnetic nanodisks with induced DzyaloshinskiiMoriya interactions. Phys. Rev. B 80, 134410 (2009).
 Mruczkiewicz, M., Krawczyk, M. & Guslienko, K. Y. Spin excitation spectrum in a magnetic nanodot with continuous transitions between the vortex, Blochtype skyrmion, and Néeltype skyrmion states. Phys. Rev. B 95, 094414 (2017).
 fu Chen, Y. et al. Nonlinear gyrotropic motion of skyrmion in a magnetic nanodisk. J. Magn. Magn. Mater. – (2018).
 Rohart, S. & Thiaville, A. Skyrmion confinement in ultrathin film nanostructures in the presence of DzyaloshinskiiMoriya interaction. Phys. Rev. B 88, 184422 (2013).
 Bogdanov, A. & Hubert, A. The properties of isolated magnetic vortices. Phys. Status Solidi B 186, 527–543 (1994).
 Donahue, M. J. & Porter, D. G. OOMMF Userâs Guide, Version 1.0. Tech. Rep., National Institute of Standards and Technology, Gaithersburg, MD (1999).
 Vansteenkiste, A. et al. The design and verification of MuMax3. AIP Adv. 4, 107133 (2014).
 CortésOrtuño, D. et al. Fidimag. Zenodo doi:10.5281/zenodo.594647. Github: https://github.com/computationalmodelling/fidimag/ (2018).
 Rybakov, F., Borisov, A. & Bogdanov, A. Threedimensional skyrmion states in thin films of cubic helimagnets. Phys. Rev. B 87, 094424 (2013).
 Yosida, K. Theory of Magnetism, vol. 122 of Springer Ser. SolidState Sci. (SpringerVerlag Berlin Heidelberg, 1996), 1 edn.
 Landau, L. D. & Lifshitz, E. Statistical Physics, vol. 5 of Course of Theoretical Physics (ButterworthHeinemann, 1980), third edn.
 Bogdanov, A. N., Rößler, U. K., Wolf, M. & Müller, K.H. Magnetic structures and reorientation transitions in noncentrosymmetric uniaxial antiferromagnets. Phys. Rev. B 66, 214410 (2002).
 Lancaster, T. & Blundell, S. J. Quantum field theory for the gifted amateur (OUP Oxford, 2014).
 Sparavigna, A. Role of Lifshitz Invariants in Liq. Cryst. Materials 2, 674â698 (2009).
 Bogdanov, A. N. & Rößler, U. K. Chiral Symmetry Breaking in Magnetic Thin Films and Multilayers. Phys. Rev. Lett. 87, 037203 (2001).
 Yang, H., Thiaville, A., Rohart, S., Fert, A. & Chshiev, M. Anatomy of DzyaloshinskiiMoriya Interaction at Interfaces. Phys. Rev. Lett. 115, 267210 (2015).
 Leliaert, J. et al. Fast micromagnetic simulations on GPUrecent advances made with mumax 3. J. Phys. D: Appl. Phys. (2018).
 This argument was discussed through an internal communication with the MuMax3 team. Free boundaries are going to be available in future MuMax3 releases.
 Carey, R. et al. Hysteresis of nanocylinders with DzyaloshinskiiMoriya interaction. Appl. Phys. Lett. 109, 122401 (2016). eprint https://doi.org/10.1063/1.4962726.
 Beg, M. et al. Ground state search, hysteretic behaviour, and reversal mechanism of skyrmionic textures in confined helimagnetic nanostructures. Sci. Rep. 5, 17137 (2015).
 Camosi, L., Rougemaille, N., Fruchart, O., Vogel, J. & Rohart, S. Micromagnetics of antiskyrmions in ultrathin films. arXiv preprint arXiv:1712.04743 (2017).
 Buhrandt, S. & Fritz, L. Skyrmion lattice phase in threedimensional chiral magnets from Monte Carlo simulations. Phys. Rev. B 88, 195137 (2013).
 Charilaou, M. & Löffler, J. F. Skyrmion oscillations in magnetic nanorods with chiral interactions. Phys. Rev. B 95, 024409 (2017).
 Moon, J. J. H. et al. Spinwave propagation in the presence of interfacial DzyaloshinskiiMoriya interaction. Phys. Rev. B 88, 1–6 (2013).
 CortésOrtuño, D. & Landeros, P. Influence of the DzyaloshinskiiMoriya interaction on the spinwave spectra of thin films. J. Phys.: Condens. Matter 25, 156001 (2013).
 Belmeguenai, M. et al. Interfacial DzyaloshinskiiMoriya interaction in perpendicularly magnetized ultrathin films measured by Brillouin light spectroscopy. Phys. Rev. B 91, 180405 (2015).
 Di, K. et al. Asymmetric spinwave dispersion due to DzyaloshinskiiMoriya interaction in an ultrathin Pt/CoFeB film. Appl. Phys. Lett. 106, 052403 (2015). eprint https://doi.org/10.1063/1.4907173.
 Tacchi, S. et al. Interfacial DzyaloshinskiiMoriya Interaction in Films: Effect of the HeavyMetal Thickness. Phys. Rev. Lett. 118, 147201 (2017).
 Kumar, D., Dmytriiev, O., Ponraj, S. & Barman, A. Numerical calculation of spin wave dispersions in magnetic nanostructures. J. Phys. D: Appl. Phys. 45, 015001 (2012).
 Venkat, G., Fangohr, H. & Prabhakar, A. Absorbing boundary layers for spin wave micromagnetics. J. Magn. Magn. Mater. 450, 34 – 39 (2018). Perspectives on magnon spintronics.
 Hertel, R. Guided Spin Waves. In Handbook of Magnetism and Advanced Magnetic Materials, 1003–1020 (John Wiley & Sons, 2007).
 Kakay, A., Westphal, E. & Hertel, R. Speedup of FEM micromagnetic simulations with graphical processing units. IEEE Trans. Magn. 46, 2303–2306 (2010).
 CortésOrtuño, D. et al. Data set for: Proposal for a micromagnetic standard problem for materials with DzyaloshinskiiMoriya interaction. Zenodo doi:10.5281/zenodo.1174311. Github: https://github.com/fangohr/papersupplementstandardproblemdmi (2018).
 CortésOrtuño, D., Beg, M., Nehruji, V., Pepper, R. & Fangohr, H. OOMMF extension: DzyaloshinskiiMoriya interaction (DMI) for crystallographic classes T and O (Version v1.0). Zenodo doi:10.5281/zenodo.1196820. Github: https://github.com/joommf/oommfextensiondmicnv (2018).
 CortésOrtuño, D., Beg, M., Nehruji, V., Pepper, R. & Fangohr, H. OOMMF extension: DzyaloshinskiiMoriya interaction (DMI) for the crystallographic class Cnv (Version v1.0). Zenodo doi:10.5281/zenodo.1196417. Github: https://github.com/joommf/oommfextensiondmicnv (2018).
 CortésOrtuño, D., Beg, M., Nehruji, V., Pepper, R. & Fangohr, H. OOMMF extension: DzyaloshinskiiMoriya interaction (DMI) for the crystallographic class D2d (Version v1.0). Zenodo doi:10.5281/zenodo.1196451. Github: https://github.com/joommf/oommfextensiondmid2d (2018).
 Beg, M., Pepper, R. A. & Fangohr, H. User interfaces for computational science: A domain specific language for OOMMF embedded in Python. AIP Adv. 7, 056025 (2017). eprint https://doi.org/10.1063/1.4977225.