Proposal for a micromagnetic standard problem for materials with Dzyaloshinskii-Moriya interaction

Understanding the role of the Dzyaloshinskii-Moriya 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 Bloch-like skyrmions, interfacial systems such as cobalt multilayers or Pd/Fe bilayers on top of Ir(111), known for stabilising N\'eel-like skyrmions, and, recently, alloys with a crystallographic symmetry where anti-skyrmions 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 so-called 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 well-defined 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.


Introduction
In computational science so-called 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 7 . 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 [1], spin waves [2] and ferromagnetic resonance [3]. Thus far, however, standard problems for materials with Dzyaloshinskii-Moriya interaction (DMI) have not been defined in the literature. Originally, the so-called DMI was phenomenologically described by Dzyaloshinskii [4,5] to explain the effect of weak ferromagnetism in antiferromagnets, and later it was theoretically explained by Moriya [6] 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 [7,8] or at the interface of a ferromagnet with a heavy metal [9][10][11]. 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 non-trivial vortex-like magnetic configurations. Although they were theoretically predicted almost thirty years ago [7], 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 [12][13][14][15][16][17][18][19][20]. The magnetic profile of a skyrmion changes according to the kind of DMI present in the material. Well-known skyrmionic textures are Néel skyrmions, Bloch skyrmions and antiskyrmions [10,11,20,21]. The former two are named according to the domain wall-like 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 [22][23][24] or vortices [25][26][27]. 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 semi-analytical results from published studies [28,29] where possible. In this context, we test these problems using three open-source micromagnetic codes, OOMMF [30], MuMax3 [31] and Fidimag [32]. In section 3 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 4, where the DMI has a distinctive influence on the boundary conditions. Then in section 5 we test the stabilisation of skyrmionic textures in a disk geometry for different kinds of DMI. In section 7 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 6 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 [33].

The Dzyaloshinskii-Moriya interaction
The Dzyaloshinskii-Moriya interaction [4][5][6] (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 S i of the form where the vector D depends on the induced orbital moments. According to the underlying crystal system, for two atoms it is usually possible to strongly constrain the direction of D through symmetry arguments [6,9,34]. 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) [35,36]. Systems featuring LIs range from Chern-Simons terms in gauge field theories [37], to chiral liquid crystals [38], but also include magnetic systems hosting DMIs. In this latter context, DMIs are described by inhomogeneous invariants in the form of spatial variations of the magnetisation m with the structure [35,36] where i, j, kä{x, y, z}. 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 D, 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 [9] and in interfacial systems [9,10,[39][40][41] 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 [41]. A specific example of this is interfacial DMI arising from the indirect exchange within a triangle composed of two spins and a non-magnetic atom with strong SO coupling [10,39,40]. From the atomistic description of interfacial DMI it is possible to derive expressions in the continuum based on LIs as shown in  [28,42].
In general, for bulk and interfacial systems the application of the LI-based continuum theory for inhomogeneous DMI, leads to the prediction of a rich variety of non-collinear magnetic structures such as vortex configurations [7,8,29,41]. An extended discussion on the theory of the DMI and further examples are discussed in sectionS1 of the supplementary material available online at stacks.iop.org/NJP/20/113015/ mmedia.

Micromagnetic model of the DMI and 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 [8,29] sin cos , sin sin , cos , 3 x y z where (Θ, Ψ) are spherical angles. In the general case, Θ=Θ(r, f, z) and Ψ=Ψ(r, f, z) with (r, f, z) 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 z-direction. For a chiral ferromagnet without inversion symmetry, we are going to describe stable magnetic configurations in confined geometries [28] when considering symmetric exchange and DMI interactions, an uniaxial anisotropy and in some cases, an applied field. Magnetostatic interactions will not be considered because the main idea of the standard problems is the verification of the DMI. The demagnetising field is usually non-trivial to model analytically and only adds an extra level of complexity. An exception will be made to the dynamics standard problem since the theoretical solution to the problem is well known. According to this, a general description of the energy of the magnetic system is given by where A is the exchange constant, K u is the uniaxial anisotropy constant, H the applied field and the last term is the DMI energy density, which can be written as a sum of LIs [7,8,29] (see section 2). For the sake of completeness, in the following we write the DMI expressions for different crystal symmetries in LIs and vector notation. In the case of thin systems the out of plane direction will correspond to the z-direction, which means the symmetry breaking originating at the surface normal to this direction. From the vector notation the DMI expressions can be generalised to any coordinate system, which can be relevant in a finite element formulation. For a material with symmetry classT orO, the DMI energy density is specified as For a thin film with interfacial DMI or a crystal with symmetry class C nv , with n>2, and located in the x-y plane, the energy density of the DMI is For a C nv symmetry with  n 2 the in-plane components of the DM vector of equation (1) are not uniquely defined [21,43] and it would be necessary more than one DMI constant in equation (7,8). For a crystal with symmetry class D d 2 , the DMI energy density reads [7] Axially symmetric magnetic configurations that are uniform along the z-direction can be found by substituting the magnetisation m (equation (3)) into equation (4), with Θ=Θ(r). Accepted solutions for Ψ are obtained according to the structure of the DMI [7,8,29]. For the T class material, Ψ=f+j, with j=±π/2. In interfacial systems, the DMI has the structure of a C nv symmetry class material, where Ψ=f+j with j=0, π. For the D 2d symmetry class an accepted solution is Ψ=−f+π/2.

One-dimensional case: edge tilting
In a one-dimensional magnetic system in the x-direction, where xä[−L/2, L/2], and at zero field, we can simplify the expression for the energy (equation (4)) using Θ=Θ(x) and Ψ=0 (interfacial), where = Q Q ( ) m sin , 0, cos , or Ψ=π/2 (bulk), where = Q Q ( ) m 0, sin , cos . We therefore obtain the following differential equation after minimising the energy with a variational approach [28] . According to our definition of DMIs in equations (5, 6) and (7,8), the positive sign in the boundary condition refers to the interfacial and D d 2 cases and the negative sign to the T class material. We numerically solve equations (11) and (12) using the shooting method. Solutions of these equations include a quasi-uniform configuration, where spins are tilted at the boundaries, and cycloids with different number of spin spirals. Although the quasi-uniform ordering is always an equilibrium configuration, the cycloid solutions depend on the DMI strength and the length of the system. For the chosen DMI magnitude of -3 mJ m 2 and length of 100nm the quasi-uniform configuration is a metastable state, because a single spiral cycloid has lower energy for lengths above a critical value around 50nm (see sectionS3 of the supplementary material). Because we are interested in comparing the tilting at the edges of the system, which is enhanced by a strong DMI, with the semi-analytical solution, we focus on the quasi-uniform solution and take advantage of a simplified evaluation of Θ at the boundary. For this, we refer to the alternative condition for Θ at x=−L/2 or x=L/2 derived in  [28], which reads The condition (13) follows directly from the quasi-uniform solution, where at the interior of the sample Θ is constant and equal to zero, as well as its first derivative, thus fixing the constant that arises after integrating equation (11). This assumption is valid when the system is large enough by a few Δ magnitudes, which sets the characteristic length scale for the edge tilting [28]. Alternatively, in  [28] the condition is justified when the system has sufficiently large anisotropy to avoid the formation of cycloids. Furthermore, we notice that the ratio of D and the critical magnitude is around 1.1, thus the constant arising from integrating equation (11) should tend to zero, as shown in  [28]. This condition can also be seen as a significantly large cycloid period. If the length L of the system is increased the formation of cycloids would be favoured [44]. With the simplified evaluation of Θ at the boundary it is straightforward to apply the shooting method. An alternative calculation to obtain general solutions of equations (11) and (12), as employed in  [44], requires a more careful analysis of the boundary conditions.
Depending on the chirality of the system, which can be observed from the simulations, we fix the condition arcsin and vary dΘ(−L/ 2)/dx until finding a solution that satisfies . The upper sign+ refers to the interfacial case and the bottom sign-to the bulk DMI case. In figures 1(a) and (b) we compare results from the theory and simulations of the one-dimensional problem, for systems with interfacial (C nv ) and bulk (T) DMI, respectively. For testing purposes, in every case we have specified micromagnetic parameters for an artificial material with a strong DMI and uniaxial anisotropy, as shown in table 1. This material has associated an exchange length of Additionally, the characteristic parameters from equation (13) are Δ=5.7 nm and ξ=5.7 nm, which means Θ≈±0.66 rad at the boundary or m z ≈±0.75. Simulations were performed with the finite difference OOMMF, Fidimag and MuMax3 software. In our examples we used a discretisation cell of 1×1×1 nm 3 volume, whose dimensions are well below the exchange length. The final magnetic configurations were obtained by relaxing an almost ferromagnetic state, i.e. m≈(0.00, 0.11, 0.99), using either the Landau-Lifshitz-Gilbert equation or a minimisation method. The profile of the z-component and either the x-component of the magnetisation, for the case of interfacial DMI, or the y-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 discretised according to the procedure shown in appendix A. For a one-dimensional system, and since we are using common magnetic parameters, the major difference between the m profiles of systems with different type of DMI, is the orientation of the spin rotation. Therefore, for a crystal with T or D d 2 symmetry, the profile of the m y component resemble the m x profile of the interfacial DMI case, which is according to the spin rotation favoured in the T and D d 2 symmetries. Accordingly, we only show the interfacial and bulk DMI solutions in figures 1(a) and (b), respectively. In the plots, data points from the simulations are compared with the solutions of equations (11) and (12). 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 figure 1). Specifically, in the interfacial case the average relative error (between the semi-analytical and simulation curves) for the m x component is about 3.8% and for m z is about 0.3%. Equivalent magnitudes are found for the bulk DMI system. Smaller errors are obtained when reducing the mesh spacing, for instance, the previous magnitudes reduce to 0.44% and 0.03% for m x and m z , respectively, when using a mesh discretisation of 0.1 nm in the x-direction. In the case of MuMax3, a similar agreement is found when imposing periodic boundary conditions along the y-direction of the one-dimensional system (with a single repetition across the width although the number should not affect the results since we are not considering dipolar interactions) because the DMI calculation is implemented with Neumann boundary conditions [31,45] rather than free boundaries 8 (see sectionS2 of the supplementary material for a comparison when not using periodic boundaries).

Two-dimensional case: magnetisation profile of a skyrmion
It has been shown in  [28,46,47] 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 C nv crystal symmetry. (ii)The so-called bulk DMI, which favours Bloch spin rotations and is found in systems with symmetry class T or O, such as FeGe. And, recently, (iii)a DMI found in systems with symmetry class D d 2 where structures known as anti-skyrmions can be stabilised [20] (antiskyrmions have also been found in interfacial systems but they are best described within a discrete spin formalism [21,43,48], where DM vectors are determined by symmetry rules or ab initio calculations). These three DMI mechanisms can be described by a combination of LIs with a single DMI constant.
We propose a two-dimensional cylindrical system of 50 nm radius and 2 nm thickness to test the stabilisation of skyrmions using the three aforementioned DMIs, using magnetic parameters from the artificial material of section 4 (see table 2).
We summarise in figure 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 figure 2(a), which shows a good agreement between them, thus we only plot results from Fidimag in figure 2(b) (see sectionS4 in the supplementary material for a more detailed comparison). Simulations were specified with cell sizes of 2 nm×2 nm×2 nm volume. The three skyrmions (see figure 2(c)) are energetically equivalent and have a similar configuration, specifically, the magnetisation 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 [28], by minimising the right-hand side of equation (4). We compare the out of plane component of the spins, = Q m cos z , with that of the simulations by extracting the data from the spins along the disk diameter, which we show in figure 2(a). As in the one-dimensional 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 r sk by finding the value of r where m z (r)=0, and plot the radial component of the spins m r (see appendix B) located at a distance r sk from the disk centre. Since spins are in-plane at r=r sk , then Θ=π/2 and the radial component (see appendix B and equation (3)) as a function of f is . According to this, we see in figure 2(b) the simulated skyrmion radial profiles at r=r sk for the C nv , T and D d 2 symmetry class materials, which agree with the theory (shown in dashed lines and curves).
In figure 2(c) we illustrate the three different configurations. The radial component of the magnetisation is shown with a colourmap and the out-of-plane component is shown in grayscale, where white means = | | m 1 z , 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 two-dimensional 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 r sk ≈22. 03 nm. For the calculation of the skyrmion radius in the simulation results, we use a third order spline interpolation of the m z profile from figure 2(a). According to this, OOMMF and Fidimag give a radius of 21.87 nm and MuMax3 produces a radius of 22.1nm, which is a slightly better approximation to the theoretical result.    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 L D . The 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 m z =0 in every slice normal to z (see figure 3(c)). The modulation has its origin in the  ( ) yx z invariant of equation (5, 6) [16]. This term imposes spatial modulations across the z-direction of the sample, which corresponds, in our case, to the thickness direction.
We define an isolated skyrmion in a FeGe cylinder of 183 nm diameter and 21 nm thickness with its long axis (the thickness) in the z-direction (see the system illustrated in figure 4). The parameters and dimensions of the system are detailed in table 3. We simulate the cylinder using finite differences with cells of dimension 1 nm×1 nm×1 nm, thus obtaining the (x, y) coordinates of the cell centres in the [−91, 91] nm range and the z-coordinates in the [−10, 10] nm range, remembering that the magnetisation is assumed constant within a finite difference cell which spans a range of±0.5 nm in every spatial direction around its centre. Relaxing the system with an initial state resembling a Bloch skyrmion [47], and an applied magnetic field of B z =0.4 T, we  figure 3. We obtain a skyrmion tube along the sample thickness with a radius that slightly varies along the z-direction (we define as the radius where m z =0 in a x-y plane-cut). In a slice of the cylinder at z=0, we compute a skyrmion radius of r sk ≈15. 3 nm, and this value decreases towards the top and bottom surfaces down to 15.0 nm, which is negligible in the scale of the chosen mesh discretisation. We emphasize that we are not considering the demagnetising field in this problem, which can enhance this effect.
We analyse the magnetisation field profiles for different slices (different z) by plotting the components of the spins located across a layer diameter (or from the skyrmion centre), (x, y)=(0, 0) nm, up to the sample boundary, (x, y)=(91.5, 0) nm, as shown by the red dashed line in every snapshot of figure 3(c). Lines in figures 3(a) and(b) are a cubic spline interpolation of the data points, which we extrapolate up to x=91.5 nm since the outermost cell centre lies at x=91 nm. To quantify the radial modulation of spins, and because of the axial symmetry of skyrmions, we calculate the cylindrical components m r and m f (see appendix B). Consistent with the results of  [33], figure 3(a) reveals that at the middle of the sample, i.e. at the x-y plane-cut located at z=0, 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 z=±10 nm, spins obtain an extra radial component that increases linearly with the z-distance. We illustrate this effect in figure 3(b) for the bottom layer at z=−10 nm. This same effect occurs at the top layer, but with the radial component pointing inwards towards the skyrmion centre, thus m r looks like a mirror image of that of figure 3(b). Furthermore, we notice that towards the cylinder caps, after the point where m z =1, the radial component towards the boundary changes sign as m f does. Snapshots with a zoomed view of the sample for the bottom, middle and top layers are shown in figure 3(c). We show spins at the skyrmion boundary, where m z =0, coloured according to their radial component, and with the background coloured according to the m z component.
The linear dependence of the radial component m r as a function of z towards the surfaces is shown in figure 4, where we plot m r as a function of z at three different (x, y=0) positions in every layer: the centre, x=0 nm, close to the skyrmion radius r sk (according to the discretisation of the mesh), x=15 nm, and at the sample boundary, x=91 nm. These spatial positions are shown as dots in the cylinder system at the top of figure 4, with lines denoting where the data is being extracted. From the curves of figure 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 z=0 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 two-dimensional model, similar to the one used in sections 3 and 5, to describe its profile. In performing this comparison (see sectionS5 in the supplementary material) we noticed that the skyrmion in the cylinder system has a larger skyrmion radius than the model predicts. In  [33] an approximate solution is provided as an ansatz for the Ψ angle, which is based on a linear dependence on z. 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 Θ=Θ(r, z) and Ψ=Ψ(f,z), but it generates a non-trivial set of nonlinear equations to be minimised. Because of the consistent skyrmion size across z it is likely that the dependence on z in the Θ angle only appears as a weak term or a constant, which differentiates the solution from that of the two-dimensional model.
Testing this problem using the OOMMF code we obtain equivalent results with the same skyrmion radius size. In the case of MuMax3, we found slight differences when comparing the radial component m r of the magnetisation at different layers, with respect to the magnitudes computed with Fidimag, being the largest differences close to the skyrmion radius at x 0 =15 nm with an error of approximately 0.007. Furthermore, results from MuMax3 simulations produce a skyrmion with larger radii compared to Fidimag and OOMMF, with magnitudes of approximately 15.8nm at to z=0 and 15.6nm at the cylinder caps. The most likely reason for this effect is the different boundary conditions imposed to the DMI calculation, as explained in section 4. 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S6 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S7 of the supplementary material.

Dynamics: asymmetric spin wave propagation in the presence of DMI
Analysing 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 [51,52]. 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 [53][54][55]. Accordingly, a standard problem based on spin waves offers the possibility to test the DMI influence on the spin dynamics of the system.
In interfacial systems spin waves propagating perpendicular to a saturating bias field, which are known as Damon-Eshbach modes, exhibit antisymmetric behaviour [51,52]. To simulate this phenomenon, we refer to the method specified in  [2] and use values from table 4, 1. Define a thin stripe with the long axis in the x-direction and thickness along the z-direction.
2. 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.
3. 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 directionx i , We delay this signal by t 0 , and then excite the system during the time τ, saving the magnetisation field every interval of duration Δτ.
4. Using the magnetisation field files, we extract the dynamic components of the magnetisation for a chain of magnetic moments along the x-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.

5.
We save these components in a matrix, where every column is a magnetisation component, m x , m y or m z , of the spins across the spatial x-direction, and every row represents a saved time step saved in the previous step.
6. Perform a two-dimensional spatial-temporal Fourier transform of the matrix, applying a Hanning windowing function [56].
We define a stripe with magnetic parameters specified in table 4, saturating the magnetisation into the ydirection. For this sample we take into account dipolar interactions because full theoretical models that consider this interaction, exchange, anisotropy and an applied field, are well known in the literature [51,52]. Including the effect of the stray field also helps to test its effect on a system with DMI. 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 2 nm width at the centre of the stripe, with a weak periodic signal based on the cardinal sine wave function. We excite Damon-Eshbach spin waves by applying the sinc field in the x-direction for a duration of τ=4 ns. According to table 4 we save τ/Δτ=4000 steps to generate the spin wave spectrum.
The result of the spin waves simulation using Fidimag with cell sizes of 2 nm×2 nm×1 nm volume, after processing the data, is shown in figure 5. In the spectrum we compare the result using the theory of Moon et al [51] 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 second order polynomial using the data points in the range of k values in the interval [−0.299, 0.016]nm −1 . The theory predicts that the minimum is located at k=−0.103 6 nm −1 with a frequency of   12.409 8 GHz. From the Fidimag simulation we estimate the minimum at k=−0.1030 nm −1 and f=12.1495 GHz, which shows they are in good agreement. Furthermore, using the fit to the simulation data we can estimate the DMI constant from the frequency asymmetry and compare it with the value specified in the simulation. The asymmetry [51,52] is defined as Δ f=f (k)−f (−k)=ξ k, with x g p = -( ) D M 2 s 1 (see table 4 for the parameters). Therefore, computing the slope of the frequency difference from opposite wave vector numbers, and diving by the corresponding factor, we obtain the DMI constant. From the Fidimag data, we obtain a DMI of 2.67 mJ m −2 , close to the magnitude of 3.00 mJ m −2 predicted by the theory. As an extra test we can notice from  [52] that in systems with crystallographic classes T or D d 2 , the Damon-Eshbach spin waves will not be antisymmetric, however they would be for spin waves excited along the field direction.
In figure 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 [57] to avoid the reflection of spin waves, or by setting periodic boundary conditions. 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. Furthermore, a better approximation of the estimated D value from the spin wave asymmetry could be obtained by using longer simulation times τ and longer stripe lengths in order to improve the resolution of the spin wave spectrum. In addition, performing the polynomial fit using a slightly larger range of k points around the spectrum minimum can also help since the curve fit is closer to the theoretical curve. However, since the spectrum signal away from the minimum is weaker, some error is introduced to the perfect second order curve tendency and the curve fit in this case misrepresents the data points around the minimum. For example, using a range of k points in the interval [−0.299, 0.076] nm −1 , the estimated D magnitude significantly improves to D=2.91 mJ m −2 but the minimum now lies at (k,f )=(−0.106 nm −1 , 11. 93 GHz).
Using the OOMMF and MuMax3 codes for simulating spin waves in systems with DMI produced equivalent results to figure 5. With OOMMF simulations we obtained a minimum at k=−0.1029 nm −1 and f=12.1615 GHz. The DMI constant is estimated as 2.66 mJ m −2 . Simulations performed with MuMax3 produce values of k=−0.1031 nm −1 and f=12.6245 GHz for the spectrum minimum and a DMI constant of 2.62 mJ m −2 . These approximations can be improved using a finer mesh discretisation, smaller time steps and longer relaxation times for the spin waves after the excitation. Results for OOMMF and MuMax3 simulations, and details about the numerical interpolation to the curves are shown in sectionS8 of the supplementary material.

Conclusions
We have proposed four standard problems to validate the implementation of simulations of helimagnetic systems with DMI mechanisms found in crystals with C nv , T and D d 2 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 one-dimensional and two-dimensional 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 r from the skyrmion centre where = m 0 z , 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 Damon-Eshbach spin waves in a stripe with interfacial DMI (or, equivalently, a crystal with C nv 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 T 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 non-existent in a slice at the middle of the sample along the thickness direction and increases linearly towards the cylinder caps (normal to the z-direction). Additionally, it is greatest at the skyrmion radius (where m z =0), 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 semi-analytical calculations the results can be also applied to finite element code simulations. Some finite element computations with our non-publicly available software Finmag are shown in sectionS9 of the supplementary material. In addition, we compared our data with the results from a non-public finite-element code developed by R Hertel, which is an entirely rewritten successor of the TetraMag software [58,59]. These results are also shown in sectionS9, 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 [60] and OOMMF [60][61][62][63] 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 [64]. Scripts and notebooks to reproduce the problems and data analysis from this paper can be found in  [60].