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

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⁷ . 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–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–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 anti-skyrmions [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–24] or vortices [25–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 [4–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

$$\begin{eqnarray}&&{H}_{\mathrm{DM}}={\boldsymbol{D}}\cdot ({{\boldsymbol{S}}}_{1}\times {{\boldsymbol{S}}}_{2}),\end{eqnarray} \tag{ 1 }$$

where the vector ${\boldsymbol{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 ${\boldsymbol{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]

$$\begin{eqnarray}&&{{ \mathcal L }}_{{\rm{i}}{j}}^{(k)}={m}_{{\rm{i}}}\displaystyle \frac{\partial {m}_{j}}{\partial k}-{m}_{j}\displaystyle \frac{\partial {m}_{{\rm{i}}}}{\partial k}\,,\end{eqnarray} \tag{ 2 }$$

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

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]

$$\begin{eqnarray}&&{\bf{m}}=({m}_{x},{m}_{y},{m}_{z})=(\sin {\rm{\Theta }}\cos {\rm{\Psi }},\sin {\rm{\Theta }}\sin {\rm{\Psi }},\cos {\rm{\Theta }}),\end{eqnarray} \tag{ 3 }$$

where (Θ, Ψ) are spherical angles. In the general case, Θ = Θ(r, ϕ, z) and Ψ = Ψ(r, ϕ, z) with (r, ϕ, 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

$$\begin{eqnarray}&&E={\displaystyle \int }_{V}{\rm{d}}V\,\left(A\displaystyle \sum _{{i}\in \{x,y,z\}}{({\boldsymbol{\nabla }}{m}_{{i}})}^{2}-{K}_{{\rm{u}}}{({\bf{m}}\cdot \hat{z})}^{2},-{\mu }_{0}{\bf{M}}\cdot {\bf{H}}+{w}_{\mathrm{DMI}}\right),\end{eqnarray} \tag{ 4 }$$

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

$$\begin{eqnarray}&&{w}_{\mathrm{DMI}}=D({{ \mathcal L }}_{{zy}}^{(x)}+{{ \mathcal L }}_{{xz}}^{(y)}+{{ \mathcal L }}_{{yx}}^{(z)})\,=\,D{\bf{m}}\cdot ({\boldsymbol{\nabla }}\times {\bf{m}}).\end{eqnarray} \tag{ 5, 6 }$$

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

$$\begin{eqnarray}&&{w}_{\mathrm{DMI}}=D({{ \mathcal L }}_{{xz}}^{(x)}+{{ \mathcal L }}_{{yz}}^{(y)})\,=\,D({\bf{m}}\cdot {\boldsymbol{\nabla }}{m}_{z}-{m}_{z}{\boldsymbol{\nabla }}\cdot {\bf{m}}).\end{eqnarray} \tag{ 7, 8 }$$

For a C_nv symmetry with $n\leqslant 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}_{2d}$, the DMI energy density reads [7]

$$\begin{eqnarray}&&{w}_{\mathrm{DMI}}=D({{ \mathcal L }}_{{xz}}^{(y)}+{{ \mathcal L }}_{{yz}}^{(x)})\,=\,D{\bf{m}}\cdot \left(\displaystyle \frac{\partial {\bf{m}}}{\partial x}\times \hat{x}-\displaystyle \frac{\partial {\bf{m}}}{\partial y}\times \hat{y}\right).\end{eqnarray} \tag{ 9, 10 }$$

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, Ψ = ϕ + φ, with φ = ±π/2. In interfacial systems, the DMI has the structure of a C_nv symmetry class material, where Ψ = ϕ + φ with φ = 0, π. For the D_2d symmetry class an accepted solution is Ψ = −ϕ + π/2.

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 ${\bf{m}}=(\sin {\rm{\Theta }},0,\cos {\rm{\Theta }})$, or Ψ = π/2 (bulk), where ${\bf{m}}=(0,\sin {\rm{\Theta }},\cos {\rm{\Theta }})$. We therefore obtain the following differential equation after minimising the energy with a variational approach [28]

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}}^{2}{\rm{\Theta }}}{{\rm{d}}{x}^{2}}=\displaystyle \frac{\sin {\rm{\Theta }}\cos {\rm{\Theta }}}{{{\rm{\Delta }}}^{2}}\,{\rm{f}}{\rm{o}}{\rm{r}}-\frac{L}{2}\lt x\lt \frac{L}{2}\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{\rm{\Theta }}}{{\rm{d}}x}=\pm \displaystyle \frac{1}{\xi }\,{\rm{a}}{\rm{t}}\,x=-\frac{L}{2}\,{\rm{o}}{\rm{r}}\,x=\frac{L}{2},\qquad \end{eqnarray} \tag{ 12 }$$

where ${\rm{\Delta }}=\sqrt{A/K}$ and $\xi =2A/D$. 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}_{2d}$ 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\,{{\rm{m}}{\rm{J}}{\rm{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

$$\begin{eqnarray}&&\sin {\rm{\Theta }}=\pm \displaystyle \frac{{\rm{\Delta }}}{\xi }.\end{eqnarray} \tag{ 13 }$$

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 ${D}_{c}=4\sqrt{{{AK}}_{{\rm{u}}}}/\pi$ 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 ${\rm{\Theta }}(-L/2)=\arcsin (\mp {\rm{\Delta }}/\xi )$ and vary dΘ(−L/ 2)/dx until finding a solution that satisfies ${\rm{\Theta }}(L/2)\,=\arcsin (\pm {\rm{\Delta }}/\xi )$. 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 ${L}_{\mathrm{ex}}=\sqrt{2A/({\mu }_{0}{M}_{{\rm{s}}}^{2})}\approx 5.3\ \mathrm{nm}$ and a helical length of ${L}_{D}=4\pi A/| D| \approx 54.5\ \mathrm{nm}$. 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³ 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}_{2d}$ 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}_{2d}$ 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⁸ (see section S2 of the supplementary material for a comparison when not using periodic boundaries).

Zoom In Zoom Out Reset image size

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}_{2d}$ where structures known as anti-skyrmions can be stabilised [20] (anti-skyrmions 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, ${m}_{z}=\cos {\rm{\Theta }}$, 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.

Zoom In Zoom Out Reset image size

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 ϕ is ${m}_{r}({r}_{\mathrm{sk}},\phi )=\sin ({\rm{\Theta }}({r}_{\mathrm{sk}}))\cos ({\rm{\Psi }}-\phi )=\cos ({\rm{\Psi }}-\phi )$. Therefore, ${m}_{{r}_{\mathrm{sk}}}^{{C}_{{nv}}}=1$, ${m}_{{r}_{\mathrm{sk}}}^{T}=0$ and ${m}_{{r}_{\mathrm{sk}}}^{{D}_{2d}}=\sin (2\phi )$. 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}_{2d}$ 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}_{z}| =1$, 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.

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 with T or O symmetry class, skyrmions can form long tubes by propagating their double-twist modulation along their symmetry axis [14, 33, 46, 49, 50], which we will assume is along the sample thickness. Moreover, it has been shown by Rybakov et al [33] 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 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 ${{ \mathcal L }}_{{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.

Zoom In Zoom Out Reset image size

We define an isolated skyrmion in a FeGe cylinder of $183\ \mathrm{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 stabilise a skyrmion tube modulated along the thickness of the sample. A characteristic parameter of FeGe is the helical length ${L}_{D}=4\pi A/| D| \approx 69.83\ \mathrm{nm}$.

Zoom In Zoom Out Reset image size

Results from Fidimag simulations are shown in 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_ϕ (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_ϕ 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 Ψ = Ψ(ϕ, 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₀ = 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.

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–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 direction ${\hat{x}}_{i}$,

  $$\begin{eqnarray}&&{{\bf{h}}}_{\mathrm{exc}}=h\,\mathrm{sinc}(2\pi {f}_{0}(t-{t}_{0})){\hat{x}}_{{\rm{i}}}.\end{eqnarray} \tag{ 14 }$$

  We delay this signal by t₀, 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 y-direction. 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⁻¹. The theory predicts that the minimum is located at k = −0.103 6 nm⁻¹ with a frequency of 12.409 8 GHz. From the Fidimag simulation we estimate the minimum at k = −0.1030 nm⁻¹ 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 $\xi =2D\gamma {(\pi {M}_{{\rm{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⁻², close to the magnitude of 3.00 mJ m⁻² predicted by the theory. As an extra test we can notice from  [52] that in systems with crystallographic classes T or ${D}_{2d}$, the Damon–Eshbach spin waves will not be antisymmetric, however they would be for spin waves excited along the field direction.

Zoom In Zoom Out Reset image size

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⁻¹, the estimated D magnitude significantly improves to D = 2.91 mJ m⁻² but the minimum now lies at (k, f) = (−0.106 nm⁻¹, 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⁻¹ and f = 12.1615 GHz. The DMI constant is estimated as 2.66 mJ m⁻². Simulations performed with MuMax3 produce values of k = −0.1031 nm⁻¹ and f = 12.6245 GHz for the spectrum minimum and a DMI constant of 2.62 mJ m⁻². 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.

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}_{2d}$ 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}_{z}=0$, 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–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].

This work was financially supported by the EPSRC Programme grant on Skyrmionics, www.skyrmions.ac.uk (EP/N032128/1), EPSRCs Centre for Doctoral Training in Next Generation Computational Modelling, http://ngcm.soton.ac.uk~(EP/L015382/1) and EPSRCs Doctoral Training Centre in Complex System Simulation~(EP/G03690X/1), and the OpenDreamKitHorizon 2020 European Research Infrastructure project~(676541).

We acknowledge useful discussions with the MuMax3 code team and L Camosi.

Assuming a two-dimensional film positioned in the x– y plane, the energy density w for the interfacial DMI used in this study is modeled as

$$\begin{eqnarray}\begin{array}{rcl}{w}_{\mathrm{DMI}} & = & D({{ \mathcal L }}_{{xz}}^{(x)}+{{ \mathcal L }}_{{yz}}^{(y)})\\ & = & D\left({m}_{x}\displaystyle \frac{\partial {m}_{z}}{\partial x}-{m}_{z}\displaystyle \frac{\partial {m}_{x}}{\partial x}+{m}_{y}\displaystyle \frac{\partial {m}_{z}}{\partial y}-{m}_{z}\displaystyle \frac{\partial {m}_{y}}{\partial y}\right).\end{array}\end{eqnarray} \tag{ A1 }$$

The corresponding effective field of this interaction reads

$$\begin{eqnarray}\begin{array}{rcl}{{\bf{H}}}_{\mathrm{DMI}} & = & -\displaystyle \frac{1}{{\mu }_{0}{M}_{s}}\displaystyle \frac{\delta {w}_{\mathrm{DMI}}}{\delta {\bf{m}}}\\ & = & -\displaystyle \frac{2D}{{\mu }_{0}{M}_{s}}\left(\displaystyle \frac{\partial {m}_{z}}{\partial x}\hat{x}+\displaystyle \frac{\partial {m}_{z}}{\partial y}\hat{y}-\left[\displaystyle \frac{\partial {m}_{x}}{\partial x}+\displaystyle \frac{\partial {m}_{y}}{\partial y}\right]\hat{z}\right).\end{array}\end{eqnarray} \tag{ A2 }$$

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 derivative of m_z (first term inside the round brackets of equation (A2)) with respect to x is

$$\begin{eqnarray}&&\displaystyle \frac{\partial {m}_{z}}{\partial x}\approx \displaystyle \frac{{m}_{z}(+x)-{m}_{z}(-x)}{2{\rm{\Delta }}x},\end{eqnarray} \tag{ A3 }$$

where m_z(±x) is the m_z component of the closest magnetic moment at the mesh site in the ±x-direction and Δx is the mesh discretisation in the x-direction. The other contributions from neighbours in the ±x-directions are given by the third term of equation (A2). We can then collect the field terms contributed by the 4 mesh neighbours (in the {+x, −x, +y, −y} directions) of every mesh site, by writing equation (A2) as

$$\begin{eqnarray}&&{{\bf{H}}}_{\mathrm{DMI}}={{\bf{h}}}_{\mathrm{DMI}}(+x)+{{\bf{h}}}_{\mathrm{DMI}}(-x)+{{\bf{h}}}_{\mathrm{DMI}}(+y)+{{\bf{h}}}_{\mathrm{DMI}}(-y).\end{eqnarray} \tag{ A4 }$$

For instance, the field contribution from the +x mesh neighbour is

$$\begin{eqnarray}\begin{array}{rcl}{{\bf{h}}}_{\mathrm{DMI}}(+x) & = & -\displaystyle \frac{2D}{{\mu }_{0}{M}_{s}}\displaystyle \frac{1}{2{\rm{\Delta }}x}({m}_{z}(+x)\hat{x}-{m}_{x}(+x)\hat{z})\\ & = & -\displaystyle \frac{2D}{{\mu }_{0}{M}_{s}}\displaystyle \frac{1}{2{\rm{\Delta }}x}([\hat{z}\times \hat{x}]\times {\bf{m}}(+x)).\end{array}\end{eqnarray} \tag{ A5 }$$

The contribution from the other neighbours have the same structure except the denominator for the neighbours in the y-direction will have a factor of 2Δ y instead of 2Δ x. We notice from equation (A5) that the cross product is, in general, given by $(\hat{z}\times {\hat{r}}_{{\rm{i}}{j}})\times {\bf{m}}$, with ${\hat{r}}_{{\rm{i}}{j}}$ the unit vector directed from the i mesh site towards the position of the neighbour in the j-direction. Hence, the calculation for the DMI field can be seen as that of the discrete spin model with an equivalent DM vector of the form ${{\bf{D}}}_{{\rm{i}}{j}}=(\hat{z}\times {\hat{r}}_{{\rm{i}}{j}})$, i.e. the field contribution from mesh site i can be computed as

$$\begin{eqnarray}&&{{\bf{H}}}_{\mathrm{DMI}}^{({\rm{i}})}=\displaystyle \sum _{j\in \{\pm x,\pm y\}}{\xi }_{j}\ {{\bf{D}}}_{{\rm{i}}{j}}\times {{\bf{m}}}_{j},\end{eqnarray} \tag{ A6 }$$

with ${\xi }_{j}=-2D{(2{\mu }_{0}{M}_{s}{\rm{\Delta }}{x}_{j})}^{-1}$ and Δ x_j as the mesh discretisation in the x_j-direction. The discretised form of the DMI field given in equation (A6) confirms the original mathematical structure of the DMI Hamiltonian from equation (1).

In the case of a T class material, the finite difference discretisation leads to a calculation of the micromagnetic DMI field with a vector ${{\bf{D}}}_{{\rm{i}}{j}}={\hat{r}}_{{\rm{i}}{j}}$. For a ${D}_{2d}$ symmetry this vector is $-{\hat{r}}_{{\rm{i}}{j}}$ for the neighbours in the x-directions and ${\hat{r}}_{{ij}}$ for the neighbours in the y-directions.

The cylindrical components of the magnetisation are computed with a transformation matrix according to

$$\begin{eqnarray}\left[\begin{array}{ccc}\cos \phi & \sin \phi & 0\\ -\sin \phi & \cos \phi & 0\\ 0 & 0 & 1\end{array}\right]\left[\begin{array}{c}{m}_{x}\\ {m}_{y}\\ {m}_{z}\end{array}\right]=\left[\begin{array}{c}{m}_{r}\\ {m}_{\phi }\\ {m}_{z}\end{array}\right],\end{eqnarray} \tag{ B1 }$$

where $\phi =\arctan (y/x)$ is the azimuthal angle.