This Tutorial article focuses on magnetic phenomena and material systems that have gained significant importance since the original development of mumax3, but are challenging to simulate for users who rely solely on the originally provided examples. Alongside the physical background, we provide hands-on examples of advanced magnetic systems, including detailed explanations of complete mumax3 input files (13 in total, often showing different ways to achieve things), and highlighting potential pitfalls where applicable. Specifically, we explore two approaches to incorporate spin–orbit torques in mumax simulations, considering the trade-off between versatility and speed. We also examine complex multilayer material stacks, including synthetic antiferromagnets, demonstrating different implementation methods that again vary in speed, versatility, and realism. A key criterion for selecting the optimal simulation strategy is its suitability for modeling systems where the magnetization varies significantly in the third dimension. The material covered in this Tutorial paper includes content developed for the mumax3 workshop presented during the summer of 2020 within the context of the IEEE online spintronics seminar, along with additional new topics. Throughout the explanations, we ensure broad applicability beyond specific examples.

Magnetic nanostructures hold significant potential as low-power components for information and communication technologies (ICT). However, to harness this potential and establish themselves as competitive platforms in the post-CMOS (complementary metal–oxide–semiconductor) era, efficient methods for controlling their magnetization states are essential. One approach to exerting torque on the magnetization is through the interaction with an injected spin-polarized current, known as spin torque. In particular, spin–orbit torque (SOT) has garnered significant attention as a versatile mechanism for manipulating spins without the need for a second magnetic layer or external magnetic fields.1 This distinguishes SOT from spin-transfer torque (STT), where angular momentum is transferred via a spin-polarized current between two noncollinear magnetic domains or separate magnetic layers.2,3 As spin-transfer torques were already described in detail in earlier mumax documentation,4 this Tutorial focuses on spin–orbit torque-driven phenomena. Spin–orbit torques can emerge from several microscopic phenomena, which all have in common that spin–orbit interaction, i.e., the relativistic coupling between the spin and orbital angular momenta of individual electrons,5 plays a key role. Two classes of spin–orbit torques are usually distinguished.

The first type of spin–orbit torque emerges in systems that lack inversion symmetry. For centrosymmetric crystals, spin degeneracy of the electronic eigenstates is dictated by time reversal symmetry, even in the presence of spin–orbit coupling.6 However, for crystals without inversion symmetry, spin–orbit coupling causes the spin degeneracy of the electronic states to be lifted, leading to spin-momentum locking. This is called Dresselhaus or Rashba spin–orbit coupling, after those who first explored symmetry-constrained spin–orbit Hamiltonians, when the atomic positions of the bulk crystal itself or external potentials are responsible for the breaking of inversion symmetry, respectively.6–8 For a non-magnetic material, spin-up and spin-down states have an equal occupation, however, when an electric field is applied, the spin-momentum locking will lead to a net spin-polarization of the resulting charge current without the need of magnetic fields, known as the Rashba–Edelstein effect.9 

The second type of spin–orbit torque is widely known as the spin-Hall effect, predicted by Dyakonov and Perel in 1971,10 and experimentally verified in 2004,11,12 after receiving renewed attention.13,14 It consists of spin accumulation at the lateral boundaries of a conductor due to the coupling of charge and spin currents by the spin–orbit interaction.11,12

It has been shown that both types of spin–orbit torques give rise to a field-like and damping-like torque acting on local magnetic moments. Regardless of their qualitative similarity, there are important quantitative differences between both types of spin–orbit torque. Specifically, the spin-Hall effect gives rise to dynamics typically dominated by the damping-like torque.15,16 In contrast, Rashba–Edelstein effects are known to yield magnetization dynamics which is mostly driven by a field-like torque.17 

In general, both types of spin–orbit torques occur simultaneously. The magnitudes and signs of these torques strongly depend on the specific details of the heterostructures, such as the layer thicknesses and compositions,18 which can complicate the design of these structures. Therefore, having access to micromagnetic simulations becomes highly desirable as a complementary tool to experimental investigations of spin–orbit torques. Such simulations enable a quantitative assessment of spin–orbit torque-driven dynamics, offering insights into the microscopic working mechanisms and guiding the further development of magnetic and spintronic devices and technologies.

Moreover, when exploring experimentally relevant systems involving spin–orbit torques, particularly those related to the spin-Hall effect, dealing with complex multilayers becomes necessary. These multilayers provide an effective means of increasing spin–orbit torques by incorporating a larger number of symmetry-breaking interfaces. For example, stacking ferromagnetic (FM) and non-magnetic heavy metal layers (HM) in a [HM1/FM/HM2] × N pattern can enhance the torques. Consequently, studying stacked layers in micromagnetic simulations inevitably requires three-dimensional models. While effective two-dimensional models may suffice in some cases, there are instances where a full three-dimensional consideration is essential. This Tutorial highlights both approaches and discusses the necessary conditions for their validity.

Although the link between three-dimensional systems and spin–orbit torques has served as a thematic bridge, it is worth noting that three-dimensional magnetic textures are inherently fascinating in their own right.19,20 By presenting examples of full three-dimensional systems, this Tutorial provides insights into the emerging field of twisted three-dimensional magnetic structures.

This Tutorial paper celebrates the 10th anniversary of the first release (in 2013) of mumax3. Mumax is a powerful NVIDIA CUDA21 based, graphics-processing-unit (GPU) accelerated, micromagnetic software package which is freely available under GPLv3 licence.4,22,23 During this decade, which we consider to be an exceptionally long time for a scientific software package, (micro)magnetic research has advanced indomitably. Mumax’ versatility is one of the reasons for its continued adoption in the community, and yet several of the major recent magnetic systems are not straightforward to simulate by relying on the originally provided examples alone. To bridge this gap, this Tutorial aims to demonstrate how mumax3 can be effectively utilized to simulate spin–orbit torques along with other modern magnetic phenomena and material systems.

The Tutorial begins with a brief overview of micromagnetism that is implemented in mumax3 (Sec. II), providing readers with the necessary background to understand the subsequent discussions and examples. The focus then shifts to spin–orbit torques (Sec. III), showcasing how they can be incorporated into mumax3 simulations as custom effective fields. Additionally, an alternative implementation is presented, leveraging the qualitative similarity with the well-known Slonczewski spin-transfer torque.

To illustrate the practical application of these simulation techniques, several example problems are presented. These examples are carefully selected to be relevant to the field of spintronics, covering a range of complexities. They include a relatively simple skyrmion racetrack (Sec. IV), as well as more complex magnetic multilayer stacks, like synthetic antiferromagnets (Sec. V), which are often encountered in spin–orbit torque-based devices. These examples not only highlight the capability of mumax3 to simulate various spintronic systems but also provide valuable insights into their behavior.

Next, the Tutorial explores sophisticated spin structures with changing magnetization texture in the third dimension (Sec. VI) and the role of the stray fields on these magnetic arrangements (Sec. VII). This extension to three-dimensional magnetic textures is of great interest as it reflects the emerging field of (twisted) three-dimensional magnetic structures. Furthermore, the Tutorial provides a concise introduction to the high-frequency internal dynamics of magnetic quasi-particles (Sec. VIII), adding depth to the understanding of their response to external excitations, enabling researchers to make quantitative comparisons with experimental data of these systems.

Throughout the Tutorial, the mumax3 input files are thoroughly discussed, ensuring that readers have a comprehensive understanding of the simulation setup and parameters.

The complete input files are made available as supplementary material, facilitating practical implementation and further exploration.

The foundations of micromagnetic theory originate from the 1930s24 and has been cast in the framework used today in the 1960s.25,26 However, it is only recently that large-scale application to nontrivial problems became possible thanks to leaps in computing power, most notably the advent of GPU-accelerated calculations.4,23,27

Magnetization dynamics is described by the Landau–Lifshitz–Gilbert (LLG) equation,24,25
m ˙ = γ 0 m × H eff + α m × m ˙ ,
(1)
where γ 0 = 2.21 × 10 5 m A s 1 and α denotes the phenomenological Gilbert damping parameter. In micromagnetic simulations, the magnitude of the magnetization is normalized to the saturation magnetization, m ( r , t ) = M ( r , t ) / M sat, where m is a vector with unit length along the magnetization direction with m ˙ its time-derivative. The effective field H eff is the functional derivative of the energy functional,
H eff = 1 μ 0 δ E δ M ,
(2)
with E being the volume integral over the local energy densities ε,
E = ε d V .
(3)
H eff has multiple terms that either are intrinsic to the magnetic material or which are imposed externally by the sample geometry, the embedding in heterostructures, the applied charge or spin currents or external fields.

All energy terms (or derived effective fields) are user-specified by means of phenomenological parameters. Intrinsic materials properties include the exchange stiffness ( A ex28,29) and the magnetocrystalline anisotropy [ K u i ( i = 1 , 2 ), K c i ( i = 1 , 2 , 3 ), for uniaxial or cubic anisotropy, respectively, where higher-order terms are in many cases neglected30]. Furthermore, an effective Dzyaloshinskii–Moriya interaction (DMI) originates from spin–orbit coupling in non-centrosymmetric ferromagnetic films.31,32 An external magnetic field ( H ext) can also be applied. Finally, the demagnetization field, i.e., the stray field originating from the magnet with saturation magnetization M sat itself,33 is calculated from the magnetization and included in the total effective field. This term gives rise to the shape anisotropy.

A more detailed discussion on the various contributions to the micromagnetic energy density and effective field can be found in Refs. 23, 34, and 35.

Additionally, we would like to mention the capabilities of mumax3 in simulating systems at nonzero temperatures by incorporating a stochastic thermal field into the effective field.36,37 Furthermore, the simulation of elasto-magnetodynamics is possible using mumax3, as presented by Vanderveken et al.38,39

Spin-polarized currents will also exert additional torques on the magnetization, requiring additional terms to be added to Eq. (1).2,3 The terms originating from various spin-transfer and spin–orbit torques will be discussed in more detail in Sec. III.

Micromagnetic simulation packages, such as mumax3, enable the simulation of magnetization dynamics through the numerical integration of the (LLG) equation, Eq. (1). They also allow for static energy optimization, where the magnetization state that minimizes the total energy is determined. mumax3 offers three distinct commands to accomplish these tasks: run, relax, and minimize.

The run command performs the numerical integration of the LLG equation by employing an explicit Runge–Kutta scheme. It proceeds by calculating the magnetization state at each time step, using the current magnetization state as the initial condition, and solving Eq. (1). This process is repeated iteratively, advancing the simulation in time.

The relax command in mumax3 is designed to find an energy minimum for the system (as described in detail in Vansteenkiste et al.4). It achieves this by disabling the precession term in the LLG equation [Eq. (1)], thereby considering only the damping term. The effective field resulting from the damping term directs the magnetization toward a direction with lower energy, although it may not necessarily be the direction of the absolute lowest energy.

The relaxation process consists of two steps. In the first step, the LLG dynamics are integrated until the total energy reaches the noise floor. This helps to establish an initial state closer to an energy minimum. In the second step, the maximum torque (which is less affected by noise) is monitored while gradually decreasing the maximum error tolerance (maxerr) in a series of steps until it reaches a value of 10 9. This unitless quantity is defined as the product of a torque (in units 1/time, as the magnetization is unitless) with the time step of the integrator. By progressively reducing the error tolerance, the relax command refines the magnetization state and attempts to converge toward an energy minimum.

Finally, the minimize command also aims to minimize the energy of the system. It employs a more efficient steepest descent method developed by Exl et al.40 This method is faster than the one used by the relax function. However, it may be less robust when starting far from an energy minimum. In such cases, the relax command is generally considered a safer choice to ensure convergence.

The effect of spin-transfer88 and spin–orbit torques41 can be described by the Landau–Lifshitz–Gilbert equation, provided that a contribution is added to the effective field,42,43
m ˙ = γ 0 m × ( H eff + H ST ) + α m × m ˙ .
(4)
H eff is the conventional effective field, originating from exchange, anisotropy, demagnetization and Dzyaloshinskii–Moriya interactions (see Sec. II). The presence of spin torque consists of a damping-like and field-like term,2,3
H ST = a J m × p + b J p .
(5)
Herein, a J and b J quantify the magnitude of the damping-like and field-like terms whereas the unit vector p characterizes the direction of the spin-polarization.
By straightforward mathematical manipulation, the right-hand side of the LLG equation can be rewritten as a torque-like expression. For this, both sides of Eq. (4) are cross-multiplied by m ×, leading to
m × m ˙ = γ 0 m × ( m × H eff ) + γ 0 a J m × p γ 0 b J m × ( m × p ) α m ˙ ,
(6)
where some standard vector algebra identities were used in addition to m m ˙ = 0. The latter is a consequence of the micromagnetic constraint that | m | = 1 remains constant. Finally, expression (6) can be substituted for the Gilbert damping term in Eq. (4), yielding
( 1 + α 2 ) m ˙ = γ 0 m × H eff γ 0 α m × ( m × H eff ) γ 0 a J m × ( m × p ) γ 0 α b J m × ( m × p ) γ 0 b J m × p + γ 0 α a J m × p .
(7)
The first line of this equation is the standard Landau–Lifshitz torque, including Gilbert damping. The second and third lines represent the additional contributions originating from spin torques, τ ST,
τ ST = γ 0 ( 1 + α 2 ) a J [ ( 1 + α ξ ) m × ( m × p ) + ( ξ α ) ( m × p ) ] ,
(8)
where the ratio of the magnitudes of the field-like and damping-like terms was introduced,
ξ = b J a J .
(9)

Although spin–orbit torques are caused by a current parallel to the plane of the magnetic material and Slonczewski spin-transfer torque by a current perpendicular to the plane of the magnetic layer, they both produce the same spin-polarization. The above equations for spin torque, expressed as effective field [Eq. (5)] or—equivalently—as torque [Eq. (8)], are, thus, suitable to describe spin-transfer torque as well as spin–orbit torque in a micromagnetic framework. The scalar parameters a J and b J and the vector p receive a different meaning for spin-transfer torque or spin–orbit torque.

In 1996, Slonczewski and Berger were the first to find that magnetic states can be modified by spin-polarized electric currents.2,3 This phenomenon was quickly valued and applied in spintronic components such as Magnetoresistive Random Access Memory (MRAM). Typical geometries of such devices consist of at least three layers: a fixed ferromagnetic layer that polarizes the current, a conducting non-magnetic spacer layer, and a free ferromagnetic layer on which the spin-transfer torque is exerted. It is only the latter layer whose magnetization is explicitly modeled. Elaborating on the original work of Slonczewski and Berger, it can be derived that the spin-transfer torque on the free layer can be modeled by44–47,
τ SL = β ϵ α ϵ 1 + α 2 ( m × ( p × m ) ) β ϵ α ϵ 1 + α 2 ( m × p ) ,
(10)
with
β = j z M sat e d ,
(11)
ϵ = P ( r , t ) Λ 2 ( Λ 2 + 1 ) + ( Λ 2 1 ) ( m p ) .
(12)
This equation is identical to the one implemented in OOMMF48 (Object Oriented MicroMagnetic Framework). Herein, the parameters Λ, d, and ϵ characterize the spacer layer and the interfaces whereas j z, P, and p characterize the spin-polarized current. d is the spacer layer thickness and ϵ was introduced to allow asymmetry between the properties of the fixed and free layers.46 Measuring the value of ϵ for experimental systems is far from trivial, but is described in detail in Ref. 49, where a value of ϵ / ϵ 0.3 is reported for the system under study, as determined by the ratio between the Slonczewski torque in the perpendicular ( m × p ) and in-plane [ m × ( p × m ) ] directions.50 

The value of Λ grasps deviations from the idealized scenario ( Λ = 1), as originally studied by Slonczewski, and is impacted by a spin-dependent resistance of the spacer layer, contact resistances, spin decoherence, etc. It is not straightforward to obtain a numerical value for Λ a priori due to the multiple origins of these nonidealities which are very system-dependent.51 Finally, j z, P, and p characterize the density of the electric current (whose direction is by definition opposite to the direction electrons move in), the spin-polarization, and polarization direction of the electrons, respectively. Consequently, p also represents the magnetization direction of the fixed layer.

If ϵ = ξ ϵ and Λ = 1 (i.e., ϵ = P / 2), then the Slonczewski torque, Eq. (10) adapts exactly Eq. (8) provided the following identification:
a J = β P 2 = j z 2 e d M sat P .
(13)

It has been shown that the spin-Hall and Rashba–Edelstein effects give rise to damping-like and field-like torques on the magnetization,15–17 suggesting that Eq. (5) is applicable for spin–orbit torques.

If it is assumed that the damping-like torque is predominantly caused by the spin-Hall effect, Eq. (13) can be modified accordingly. The magnitude of the spin-Hall angle α H gives, by definition, the maximal efficiency by which a charge current density ( j c) is converted into a transverse spin current density52 ( j s), i.e., j s = | α H j c | / e, hence substituting P j z for | α H j c | in Eq. (13) immediately yields the corresponding damping-like torque magnitude resulting from the spin-Hall effect,53 
a J = | j c | 2 e d M sat | α H | .
(14)
The direction of the injected magnetic moments, p, is determined by the sign of the spin-Hall angle, i.e.,
p = sgn ( α H ) j c × n ,
(15)
where n represents the surface normal of the interface between the ferromagnet and the heavy metal, pointing toward the ferromagnet.
Manchon and Zhang derived an expression for the field-like torque caused by the Rashba–Edelstein effect in a two-dimensional ferromagnetic film in an asymmetric stack.54 They applied the Boltzmann equation to the stationary eigenstates of the single-electron Rashba Hamiltonian featuring spin–orbit (with coupling strength α R) and exchange (with coupling strength J) interactions, leading to
b J = P | j c | μ B M sat α R ,
(16)
where P = J / ϵ F is the spin-polarization of the current in the ferromagnetic layer, with ϵ F being the Fermi level of the mean-field model and μ B the Bohr magneton.
The above expressions for a J and b J lead to
ξ = 2 e d P μ B α R | α H | .
(17)

In general, the parameters a J and b J show a dependence on the relative orientation of the magnetization and the polarization of the current in heterostructures, but in practice they can often be considered as constants for a given heterostructure when the magnetization angle barely changes.42 Values for a J and b J have, for instance, been derived experimentally by fitting harmonic Hall voltage measurements with analytic expressions following from Eq. (4).42 

Skyrmions are localized metastable magnetization textures with a non-trivial topology55 and can be thought of as magnetic “bubbles” in a uniformly magnetized background. They are currently considered in a plethora of applications,56–59 several of which are based on current-driven motion of skyrmions. Spin–orbit torques have been proposed as an effective means to transport magnetic skyrmions, because the interfaces which lead to the large DMI necessary to stabilize skyrmions also lead to a strong spin–orbit interaction. Compared to conventional spin-transfer-torque driven skyrmion propulsion, spin–orbit torque-driven propulsion is more robust, enabling significantly higher driving currents and hence higher skyrmion velocities to be reached, e.g., within a racetrack.60,61

The simulation script below is largely inspired by the work in Ref. 62 where regular spin-transfer torque is used to propel Néel skyrmions in a Co/Pt strip. The platinum layer induces chirality in the cobalt ferromagnet via interfacially induced DMI, as needed to stabilize the skyrmion (specified by Dind63 in J/m 2). Furthermore, a perpendicular magnetic anisotropy is included [specified by AnisU (direction) and Ku1 (magnitude), in J/m 3], following Ref. 62. Additional micromagnetic parameters are M sat (specified by Msat, in A/m), α (specified by alpha, dimensionless), and A ex (specified by Aex, in J/m). A Co strip of dimensions 256 × 64 × 1 nm 3 is simulated with a skyrmion in the initial magnetization,

Two equivalent implementations of spin–orbit torque are illustrated.

This section is accompanied by supplementary material mumax input file 1a.txt.

Spin–orbit torques are not explicitly implemented in mumax3, but because they can be written as an effective field [Eq. (5)], the custom fields functionality (see  Appendix) can be used.

Parameters corresponding to the Co/Pt interface with charge current directed along the negative x axis [see also Fig. 1(a)], as well as some fundamental constants are defined,

FIG. 1.

Skyrmion racetrack. Strip geometry (a), trajectory (b), and selected magnetizations (c) of a spin–orbit torque-driven skyrmion on a racetrack.

FIG. 1.

Skyrmion racetrack. Strip geometry (a), trajectory (b), and selected magnetizations (c) of a spin–orbit torque-driven skyrmion on a racetrack.

Close modal

Herein, SOT is added to the notations of xi_SOT and J_SOT to distinguish them from xi and J which refer to other (STT related) quantities in mumax3. In the considered case of the skyrmion racetrack, p is directed along the negative y axis [see Eq. (15)] for the given stack [ sgn ( α H ) > 0 for Pt] and current (see also Fig. 1). A phenomenological value of 2.0 is used for ξ [Eq. (9)]. Finally, a J and b J are obtained from Eqs. (14) and (9),

Subsequently, both additional field contributions H add are defined and added to the effective field that is calculated by mumax3,

Herein, M_full represents the unnormalized magnetization, i.e., M sat m.

This section is accompanied by supplementary material mumax input file 1b.txt.

While mumax3 has a custom fields functionality, no custom torque functionality is available. However, as illustrated in Sec. III B, τ SOT has the same mathematical form as Slonczewski STT, τ S L, which is readily implemented in mumax3. Provided that the Slonczewski STT parameters are translated according to the derived recipe, the existing SST implementation can be directly used to simulate SOT. In practice, comparison of Eqs. (13) and (14) shows that | α H | must be identified with P. The vector p is specified by the keyword Fixedlayer. This leads to following declaration of parameters in the mumax3 input:

It should be remarked that the current J is specified here along the z direction which is due to the different geometries of Slonczewski spin-transfer torque (as implemented in mumax3) and spin–orbit torque. When the latter is simulated, this setting corresponds to a current along the negative x axis, as shown in Fig. 1.

The position of the skyrmion can be tracked by adding ext_bubblepos to the output table,

As expected, identical results are obtained for both simulations; however, the implementation using Slonczewski STT ran almost twice as fast as the second implementation using custom fields (2 vs 3.5 min using an NVIDIA Titan RTX GPU). This can be attributed to the non-simultaneous evaluation of the arithmetic operations in separate CUDA calls in the case of the custom fields. This contrasts the calculation of the Slonczewski STT, which is implemented in a single CUDA kernel. Figures 1(b) and 1(c) show the trajectory of the skyrmion on the racetrack. During the first 0.5 ns of the simulation, the skyrmion movement features a significant transverse component, moving about 10 nm downward while moving longitudinally along the racetrack, i.e., the skyrmion Hall effect.64 Subsequently, a reasonably longitudinal trajectory is followed where the repulsion from the edge of the simulation box and the skyrmion Hall effect equilibrate.65,66 Finally, shortly after 2 ns (not included in the simulation), the skyrmion collides with the edge nevertheless, leading to its annihilation. Several methods have been recently proposed to suppress the skyrmion Hall effect.67–72 In Sec. V, one particularly elegant solution will be simulated using mumax3.

In antiferromagnetically coupled magnetic layers, two oppositely magnetized skyrmions can be propelled as a single entity. As the sign of the skyrmion Hall angle depends on its magnetization direction, the transverse components of the force on the individual Néel skyrmions will cancel each other out in this synthetic antiferromagnet, so only the desired longitudinal motion remains.73–75 

Experimentally, synthetic antiferromagnets are fabricated by adding dedicated non-magnetic spacer layers of selected materials with well-tuned thickness between both ferromagnetic layer (see Fig. 2). For specific thicknesses this can result in an antiferromagnetic coupling due to the RKKY interaction. In analogy with the non-magnetic heavy metal layers, also the non-magnetic spacer layer is not explicitly modeled in mumax3. Below, we discuss two methods of simulating such systems. In the first implementation, we only consider the negative interlayer exchange and further neglect the spacer layer entirely (i.e., assume it has zero thickness), whereas in the second implementation, we explicitly add it as an empty space in the simulation to increase the accuracy of the magnetostatic field calculation between both layers. The latter results in a significantly slower simulation, and leads to almost identical results for the specific system simulated here, as shown below.

FIG. 2.

Synthetic antiferromagnet. Geometry (a), skyrmion position within the simulated geometry (b) and cut through (through the center of the wire, such that only half the skyrmion is visible) of the magnetization of a spin–orbit torque-driven skyrmion on a synthetic antiferromagnet for both implementations, i.e., touching layers (c) and separated layers (d). The zig-zag motion of the position within the simulation originates in the periodic boundary conditions used. In this case, both implementations give rise to the same skyrmion motion and identical magnetization structures within the magnetic layers. (c) and (d) were made with mumax-view,76 a freely available viewer for mumax output.

FIG. 2.

Synthetic antiferromagnet. Geometry (a), skyrmion position within the simulated geometry (b) and cut through (through the center of the wire, such that only half the skyrmion is visible) of the magnetization of a spin–orbit torque-driven skyrmion on a synthetic antiferromagnet for both implementations, i.e., touching layers (c) and separated layers (d). The zig-zag motion of the position within the simulation originates in the periodic boundary conditions used. In this case, both implementations give rise to the same skyrmion motion and identical magnetization structures within the magnetic layers. (c) and (d) were made with mumax-view,76 a freely available viewer for mumax output.

Close modal

This section is accompanied by supplementary material mumax input file 2a.txt.

Here, both ferromagnets are specified in mumax3 as adjacent regions that exhibit a mutual exchange interaction. Skyrmions of opposite magnetization are subsequently added,

Propulsion of the bilayer skyrmion is achieved by spin–orbit torque, implemented via the built-in Slonczewski spin-transfer torque as explained in Sec. IV B. The RKKY interaction between the two layers is modeled by changing the bulk exchange interaction. The ext _InterExchange function (value in J/m) locally overrules the pre-specified homogeneous value of A ex, leading -in this case due to the negative sign- to antiferromagnetic coupling between the two layers. This value is related to the RKKY surface energy density as follows:
A ex , RKKY = c z J RKKY 2 ,
(18)
with c z the cellsize in the z direction. The used value of A ex , RKKY = 5 × 10 13 J / m hence corresponds to a surface energy density of J RKKY = 0.001 J / m 2. The use of this functionality requires both layers to touch which is not achievable in the case where the antiferromagnetic coupling is achieved by an RKKY exchange interaction through a spacer layer of a well-tuned thickness. As a result, the magnetostatic field as calculated by mumax3 will not be physically correct. Often, the deviation will be negligible, but nonetheless, if necessary, mumax3 is capable of explicitly including an empty space between both ferromagnetic layers, as shown below. As shown in Ref. 77, this approximate way to model RKKY interactions is valid for thin layers (which is the case in our example), but fails if perpendicularly inhomogeneous states are formed in thick magnetic layers.

This section is accompanied by supplementary material mumax input file 2b.txt.

When both ferromagnetic layers are separated by a spacer layer, the interlayer exchange functionality cannot be used. An alternative strategy is hence needed to implement the antiferromagnetic RKKY exchange coupling.

First, both ferromagnetic layers are physically separated in the mumax3 input,

As can be seen from the setgridsize(256,64,3) command, the simulation now contains 3 layers: layers 0 and 2 are subsequently defined as a region which is added to the geometry, whereas layer 1 of the simulation box corresponds to the spacer layer. By default, mumax3 assumes that the entire simulation contributes to the thickness of the magnetic layers in calculating the spin torque. Here, this is not the case due to the presence of the spacer layer. The total thickness of the ferromagnetic layers must therefore be specified by the freelayerthickness keyword.

Subsequently, the exchange between both layers is explicitly implemented as an effective field:
H exch = 2 A ex μ 0 M sat i ( m i m ) r i 2 ,
(19)
where the sum ranges over the nearest neighbor cells in the other layer ( i , with magnetization m i) for a given cell ( m). r i is the associated distance, which boils down to the constant z component of the cellsize for the case at hand. It must be stressed that Eq. (19) only accounts for interlayer exchange, intralayer exchange is readily accounted for in Eq. (2). Equation (19) is specified in the mumax3 input and implemented via the custom fields functionality (see  Appendix). Note that mumax3 internally uses effective field terms in units of Tesla, which is why the custom field is also defined in units of Tesla,

The interlayer exchange is composed of two contributions, up and down. The specific implementation of up is visually explained in Fig. 3, exploiting the various operations that are available in mumax3. An analogous procedure is followed for down.

The position of the synthetic antiferromagnetic skyrmion is still monitored by adding the skyrmion position to the output table, which is scheduled to be written every 10 ps, and setting ext_BubbleMz=-1. This is necessary as ext_bubblepos only looks for the skyrmion position in the top layer, and the top layer has a background of m z=1, with an oppositely magnetized skyrmion, i.e., with m z = 1,

FIG. 3.

Graphical representation of how Eq. (19) can be implemented as custom field in mumax3. The vertical boxes correspond to the simulated system with a height of three cells. The dotted boxes correspond to additional zeros that are inserted as a result of a shift operation.

FIG. 3.

Graphical representation of how Eq. (19) can be implemented as custom field in mumax3. The vertical boxes correspond to the simulated system with a height of three cells. The dotted boxes correspond to additional zeros that are inserted as a result of a shift operation.

Close modal

Figure 2 shows the trajectory of the skyrmion along the racetrack for both described implementation methods.

As anticipated, the skyrmions now move longitudinally without converging towards the lateral edges of the nanostrip. Due to the used periodic boundary conditions, the skyrmion disappears at x = 128 nm to re-appear at x = 128 nm as illustrated in the figure.

Again, the implementation with the custom fields exchange coupling requires more time (7 min vs 86 sec) due to the different CUDA kernels that are used consecutively. Additionally, in this simple example, the spacer layer had the same thickness as the 2 FM layers. In reality, this is not always the case, and as explained in detail below, this will significantly slow down the simulations further due to the need to add more cells along the thickness direction.

This section is accompanied by supplementary material mumax input files 3a.txt, 3b.txt, 3c.txt, and 3d.txt.

Multilayer structures have always been challenging for micromagnetic investigations. The main bottleneck is the sheer number of cells that a three-dimensional sample description requires and the inherent limitation of finite-difference simulations to only support constant cell sizes throughout the simulation volume. This usually leads to a vast increase in the number of cells as the least common divisor defines the cell dimensions. For example, if a material is composed of [Pt(5 nm)/Co(1 nm)/Ta(3 nm)] × 15, the thickness of a cell must be 1 nm or less to accurately capture its geometry, yielding a total of at least 135 layers. Trying to simulate such systems becomes unfeasible.

Luckily, in many cases, the magnetic texture is not, or only minorly varying in the third dimension or there are layers that do not carry any magnetic moment. A typical example would be the multilayer skyrmion tube (or string) that consists of several individual skyrmion spin textures stacked repeatedly in each magnetic layer. For a typical thin-film multilayer system, this also includes many non-magnetic layers as spacers (see Fig. 4). Such a system can be simplified to achieve a massive speed up of the simulation and reducing the required calculation time from days or weeks down to manageable durations. Ideally, we want to simplify our system as shown in Fig. 4. Here, we generalized the stack to two non-magnetic heavy metal layers (HM1 and HM2) and a ferromagnetic source layer (FM) sandwiched between them. The thickness of the magnetic layer is t m, the total repetition thickness, including both non-magnetic layers, is t r. This base layer will be repeated N times. We now want to replace this full stack by a configuration, where we only have magnetic layers of thickness t m , thus being able to reduce all individual layers into a single one with thickness N t m .

To achieve this feat, we must re-write the LLG equation into an effective form that captures the energetics and dynamics of a multilayer system in only one dimension. To do this, we change the field and magnetization parameters of the LLG into primed variables. Starting from the explicit version of Eq. (4), this reads
m ˙ = γ 0 ( 1 + α 2 ) m × H eff γ ( 1 + α 2 ) α m × ( m × H eff ) .
(20)
The transformed effective medium equation then reads
m ˙ = γ 0 ( 1 + α 2 ) m × H eff γ ( 1 + α 2 ) α m × ( m × H eff ) .
(21)
Here, we assume that m ( x , y , i ) depends only on the position ( x , y ) in the plane of the individual layers. Additionally, when investigating the variation of the magnetization in different layers, the repetition i can also be added (though care must be taken to not add exchange coupling between adjacent layers if the full stack consisted of spacers without exchange coupling). These new primed parameters will then have to be able to describe the system in a stack of a homogeneous material in the form shown in Fig. 4. We, thus, remove the necessity to describe the system by small cells in the third dimension and can use only one cell of thickness N t m. The derivation of all micromagnetic terms can be readily found in the literature78 and results in the rather intuitive scaling,
M sat M sat = A ex A ex = K eff K eff = D D = t m t m = t m t r .
(22)
As a result, multiplying the respective magnetic terms with t m / t r will scale the interactions to allow for a simplified, all magnetic structure. Herein, K eff is the effective anisotropy. The uniaxial anisotropy usually used in micromagnetics can be obtained by
K u = K eff + K demag = K eff + μ 0 2 M sat 2 .
(23)
Note that the anisotropy field H K defined from K eff = μ 0 2 H K M sat is not scaled and remains unchanged in this set of scaling laws. Since spin–orbit torques (with b J = ξ = 0 in the case of pure damping-like spin-Hall torque) follow the form
τ SH = γ 0 1 + α 2 a J ( m × ( m × p ) + α ( m × p ) ) ,
(24)
where a J also scales with the distance of the individual moment from the spin accumulating interface d, and, thus, with the thickness of the magnetic layer t m, the obtained scaling factor can simply be applied in those cases by scaling the effective Hall angle strength to the total number N of driving layers and t r.78 From Eq. (14),
a J = j 2 e d M sat α H N t r .
(25)
As a result, if the assumption that the magnetization does not vary from layer to layer holds, the simulations should effectively reproduce the results obtained when simulating a full stack.
FIG. 4.

Principle of the effective medium model. A stack containing several magnetic and non-magnetic layers can be compressed into an effective magnetic material that reproduces energetics and dynamic behavior of the full material system but only requires homogeneous layers that can be compressed into a single cell in thickness. Reproduced with permission from Woo et al., Nat. Mater. 15, 501–506 (2016).78 Copyright 2016 Springer Nature Limited.

FIG. 4.

Principle of the effective medium model. A stack containing several magnetic and non-magnetic layers can be compressed into an effective magnetic material that reproduces energetics and dynamic behavior of the full material system but only requires homogeneous layers that can be compressed into a single cell in thickness. Reproduced with permission from Woo et al., Nat. Mater. 15, 501–506 (2016).78 Copyright 2016 Springer Nature Limited.

Close modal

To demonstrate the equivalence of the effective medium approach and the full 3D simulation, we will investigate the example of a disk-shaped structure with four uniformly magnetized stripe domains as initial state, alternating between z and + z directions. Starting from this state, we will relax the magnetization toward a local energy minimum state and compare the magnetic textures found using both approaches.

The input files for the effective 2D model is built as follows. First, the scaling factors, the rescaled parameters, grid, and cell-sizes are defined,

We can see that the grid only contains a single cell in the z direction. The scaling is implemented according to Eqs. (22) and (23). It corresponds to a stack consisting of 15 repetitions of a HM/FM/HM trilayer with one magnetic layer of 0.9 nm and a thickness of one trilayer of 8.1 nm, each. Next, the geometry is set to a 800 nm disk and four regions of equivalent dimensions are defined. Uniform magnetization is set as the initial state in these areas alternating in the + z and z directions,

After minimization, the magnetization can fall into a series of energetically identical ground states as shown in Fig. 5(a) and 5(b), which indeed accurately reproduce the results obtained from the full 135-layer simulation, as shown in Fig. 5(c).

If the initial state is chosen to be a perfectly symmetrical arrangement of stripe domains, for which, due to their nature, both energy and torque minimization algorithms will be unable to find the direction toward a distinctive ground state. To break this symmetry, the system was run for 500 ns before using the minimize command (see Sec. II A). As the system resides in a saddle point of the energy landscape, small changes in the magnetization state introduced by the finite numerical precision can build up during this time, eventually breaking the symmetry of the initial state. This process can be observed in Fig. 6, where a highly symmetric state (corresponding to a saddle point in the energy landscape) prevails for more than 200 ns [panel (a)]. Only after the symmetry is broken will the minimizer finally find a suitable relaxed configuration, as shown in panel (b). A more efficient method to overcome such issues can be achieved by slightly breaking the symmetry. This might be achieved by introducing a small external magnetic field,

FIG. 5.

Results of the effective medium 2D simulation with one cell in thickness [panels (a) and (b)] and the full 135 layer 3D simulation, panel (c). As clearly visible, states (b) and (c) are practically identical with small variations on the order of the domain wall width. The simulation illustrates the possibility to obtain different states with the same code (a) and (b) when the underlying geometry is sufficiently symmetric. Re-plotted using simulation parameters used in Ref. 78.

FIG. 5.

Results of the effective medium 2D simulation with one cell in thickness [panels (a) and (b)] and the full 135 layer 3D simulation, panel (c). As clearly visible, states (b) and (c) are practically identical with small variations on the order of the domain wall width. The simulation illustrates the possibility to obtain different states with the same code (a) and (b) when the underlying geometry is sufficiently symmetric. Re-plotted using simulation parameters used in Ref. 78.

Close modal
FIG. 6.

Metastability of highly symmetric states. (a) shows a time evolution of the magnetization highlighting the eventual formation of asymmetries in the system as deviations from the initial state are introduced due to the finite precision of the calculation; note the curvature in the red central domain wall which appears after 250 ns. Before these asymmetries form, the minimize solver results in an incorrect state. As an example, this is shown in (b) (left) for minimizing after letting the system evolve for 5 ns. Only after sufficient time when the numerical error has been accumulated, or after the application of a symmetry-breaking magnetic field, will the solver finally find a ground state as shown in (b) (right).

FIG. 6.

Metastability of highly symmetric states. (a) shows a time evolution of the magnetization highlighting the eventual formation of asymmetries in the system as deviations from the initial state are introduced due to the finite precision of the calculation; note the curvature in the red central domain wall which appears after 250 ns. Before these asymmetries form, the minimize solver results in an incorrect state. As an example, this is shown in (b) (left) for minimizing after letting the system evolve for 5 ns. Only after sufficient time when the numerical error has been accumulated, or after the application of a symmetry-breaking magnetic field, will the solver finally find a ground state as shown in (b) (right).

Close modal

Another possibility is to include a small amount of randomness to the initial magnetization. Note that this may not work well in multilayer simulations where the randomness can average out over a series of layers,

An other important issue one needs to be cautious about arises from compressing all layers in the thickness direction into one effective cell. For this approach to be physical, the assumption that the spin texture is identical in each layer must be met. This is not necessarily true in thick stacks or films in which long-range effects of stray fields become crucial with increasing thicknesses. The arising phenomena will be discussed in Sec. VII.

This section is accompanied by supplementary material mumax input file 4.txt.

In a thick magnetic film with interfacial DMI, which favors Néel domain walls at the surface, the internal stray fields at a domain wall can become strong enough to cause a reorientation of the wall toward a flux closure Bloch configuration and, thus, reduce the domain wall energy. The dipolar fields will generally favor Bloch wall configurations in the center of the film and Néel orientation at its surface, leading to mixed states within the material. This effect can be easily demonstrated by simulating a material with relatively low DMI and sizable saturation magnetization. An example is given below, with a 250 × 250 × 200 nm 3 material dimension discretized by 64 × 64 × 50 cells. The initial magnetization is set by defining a cylindrical region with a 100 nm diameter which is oppositely magnetized to the ferromagnetic background.

It should be noted that in bulk materials, a bulk DMI can occur instead of an interfacial DMI. This is particularly true for materials lacking inversion symmetry, such as B20 materials. However, in the context discussed here, we do not consider these bulk DMI effects. Nonetheless, it is worth mentioning that a mixed state can still be achieved in such systems, where Néel endcaps are formed at the top and bottom interfaces. The relaxed magnetization of the model system is shown in Fig. 7. It presents a cross section of the skyrmion tube and demonstrates how the skyrmion tubes’ domain wall changes its sense of rotation in the thickness direction due to the interplay between DMI and stray fields. As stated earlier and shown by the zoomed-in part in the figure, the stray fields cause a flux closure configuration to minimize the stray field energy. The choice whether to compress this magnetization into a two-dimensional simulation or to accurately model the three-dimensional variation in the magnetization is dictated by the phenomena under study (i.e., whether internal structural changes are relevant for the investigated effects or not) but unfortunately can greatly impact the simulation time for thick structures.

This section is accompanied by supplementary material mumax input file 5.txt.

Another example showcasing the impact of the demagnetizing field on the magnetic structure in stacked multilayers is depicted in Fig. 8. In this case, the presence of stray fields induces the formation of a magnetic vortex, characterized by a chirality vector in an in-plane direction. This differs from the typical vortices observed in thin-film materials, where the chirality vector points out of the film’s plane. The position of the vortex core can be shifted toward the edges of the system and eventually expelled if there is a sufficient amount of interfacial DMI present. The presence of interfacial DMI favors one Néel orientation over the other. Once the critical DMI threshold, which depends on the specific material, is surpassed, the system behaves consistently with the assumptions of an effective medium layer, as discussed in Sec. VI. However, below this threshold, it becomes necessary to simulate the system using multiple layers in order to accurately capture its dynamics. Neglecting this effect can lead to incorrect or unphysical results, particularly in scenarios involving domain wall and skyrmion motion driven by Hall torques.

FIG. 7.

A cross section of a skyrmion tube to demonstrate stray field effects on domain wall configurations in thick materials with low DMI. The domain walls’ magnetization at the cross-sectional area is shown. The magnetization direction of the domain wall rotates inside the material. Hence, the chirality of the skyrmion tube changes in every layer. The top left inset shows a zoomed area with a flux closure, which minimizes the stray field energy.

FIG. 7.

A cross section of a skyrmion tube to demonstrate stray field effects on domain wall configurations in thick materials with low DMI. The domain walls’ magnetization at the cross-sectional area is shown. The magnetization direction of the domain wall rotates inside the material. Hence, the chirality of the skyrmion tube changes in every layer. The top left inset shows a zoomed area with a flux closure, which minimizes the stray field energy.

Close modal
FIG. 8.

Stray field effects on domain wall configurations in multilayer systems. At zero DMI, the stray fields will tend to form a flux closure within the multilayer stack, leading to a rotating magnetization vortex state throughout the layers with opposite Néel configuration at the edges. A nonzero DMI will favor one of these orientations and move the vortex core out of the material if sufficiently strong.79,80 (a) shows the frontal view of the material stack. (b) shows a close-up of the middle part. Both figures were made with mumax-view,76 a freely available viewer for mumax output.

FIG. 8.

Stray field effects on domain wall configurations in multilayer systems. At zero DMI, the stray fields will tend to form a flux closure within the multilayer stack, leading to a rotating magnetization vortex state throughout the layers with opposite Néel configuration at the edges. A nonzero DMI will favor one of these orientations and move the vortex core out of the material if sufficiently strong.79,80 (a) shows the frontal view of the material stack. (b) shows a close-up of the middle part. Both figures were made with mumax-view,76 a freely available viewer for mumax output.

Close modal

When dealing with rotating spin structures, an interesting stray field phenomenon may be found in configurations forming the so-called Halbach array. This effect, coined after the physicist Klaus Halbach, is known for specific rotated arrangements of permanent magnets since the 1980s.81 It has been used for a variety of technical applications that require high and strongly localized fields, one very common example being the fridge magnet. This trivially appearing array of magnetic elements usually demonstrates the interesting peculiarity of only sticking to the fridge on one side. This is because the arrangement of magnets in a chiral fashion enhances the stray field lines on one side while suppressing them on the other, depending on whether the in-plane direction of the magnetization forming the chiral spiral helps or counteracts the flux closure of the out-of-plane components.

A depiction of the configuration is shown in Fig. 9(a). On the nanoscale, this effect can be achieved, e.g., by Néel domain walls or skyrmion textures as depicted in Figs. 9(b) and 9(c). These figures specifically show the difference of stray fields on the same surface but with reversed chirality. Generally, skyrmions in multilayers or SAFs as described in Sec. V do not necessarily need to have the same chirality. Experimentally, the sign of DMI can be tuned in each layer individually by surface-engineering. At the same time, this can also alter the globally observed stray field by the means of the Halbach effect as just described. One consequence of this consideration can be seen very clearly when simulating RKKY-coupled skyrmion tubes in a SAF structure with different chiralities. For this, let us consider a 150 × 150 × 21 nm 3 gridsize discretized by 1 × 1 × 1 nm 3 cells, with two regions comprising of ten layers each, and a 1 nm gap in-between. Again, the material parameters of Co/Pt strips will be used for each region, whereas this time two sets of simulations will be initialized. One with the DMI of both regions having the same sign and one with opposite signs. For the sake of simplicity, both configurations are shown in the same code snippet. When simulating, each configuration must be initialized in its own simulation file.

FIG. 9.

Representation of the Halbach effect at different scales. (a) illustrates a specific rotating arrangement of permanent magnets and the emerging stray field. The stray field is stronger on the top surface for this specific Halbach array. Similar effects on nanoscale can appear for magnetization textures with similar rotational senses. An example of this is shown in (b) and (c) for Néel skyrmion textures with counter clockwise and clockwise rotational senses, respectively. A cross section of the magnetization texture and the emerging stray fields are shown for both chiralities at the top interface. Figure modified and reused with permission from Bassirian et al., APL Mater. 10, 101107 (2022).82 Copyright 2022 Author(s), licensed under a Creative Commons Attribution (CC BY) license.

FIG. 9.

Representation of the Halbach effect at different scales. (a) illustrates a specific rotating arrangement of permanent magnets and the emerging stray field. The stray field is stronger on the top surface for this specific Halbach array. Similar effects on nanoscale can appear for magnetization textures with similar rotational senses. An example of this is shown in (b) and (c) for Néel skyrmion textures with counter clockwise and clockwise rotational senses, respectively. A cross section of the magnetization texture and the emerging stray fields are shown for both chiralities at the top interface. Figure modified and reused with permission from Bassirian et al., APL Mater. 10, 101107 (2022).82 Copyright 2022 Author(s), licensed under a Creative Commons Attribution (CC BY) license.

Close modal

Similar to Sec. V B, the RKKY-coupled layers are chosen not to touch each other. This is not a necessity to show effects of the Halbach array. Moreover, it also serves as an example how the approach described in Sec. V B must be adapted when dealing with SAFs comprised of layers with more than one cell in the thickness direction,

Before relaxing the system, the initial spin structures of a Néel skyrmion in each layer is implemented. The specific polarity and vorticity must be chosen according to the ferromagnetic background and the sign of DMI,

As can be seen, this approach is almost identical as described by Fig. 3. We will skip a detailed graphical representation at this point, as the changes can be easily understood when following the same logic as presented in the aforementioned figure.

The results are shown in the bottom parts of Fig. 10. When simulating two RKKY coupled layers with the same sign of DMI, the diameters of the skyrmionic structures differ in the two layers. This can be understood by asymmetric stray fields due to the Halbach effect. The skyrmion in the bottom layer has high stray fields crossing the top interface, acting as an external field source, while the skyrmion in the top layer has rather low stray fields at the bottom interface. However, when changing the sign of DMI of the bottom layer, hence changing the rotational sense of the skyrmions’ domain walls, the emerging stray fields of both layers will become symmetric. This will lead to coupled skyrmions of identical size. A comprehensive illustration showing the stray fields of each layer is given in Fig. 10. Although the change in diameter seems small, it can be crucial when studying internal modes. If the RKKY coupling is lower than assumed in this example, higher discrepancies between the skyrmion diameters can appear.

This is yet another example showing that stray fields can be important sometimes. It demonstrates how such delicate interactions as the interplay between RKKY and stray field energies can be captured by the demagnetization module of mumax3.

FIG. 10.

Skyrmion Halbach fingerprint in antiferromagnetically coupled skyrmions. Depending on the chirality of each individual skyrmion in the coupled layers, the mutual stray field impact can be either antisymmetric or symmetric. These two cases are shown in panels (a) and (b), respectively. The top figures in both panels display the stray field line density for each layer of the SAF system individually on a cross-sectionial view of the layer. The bottom figures show a cross section of the coupled structure, according to the code in the main text. A slight difference in skyrmion diameters can be observed in top and bottom layers for the asymmetric case.

FIG. 10.

Skyrmion Halbach fingerprint in antiferromagnetically coupled skyrmions. Depending on the chirality of each individual skyrmion in the coupled layers, the mutual stray field impact can be either antisymmetric or symmetric. These two cases are shown in panels (a) and (b), respectively. The top figures in both panels display the stray field line density for each layer of the SAF system individually on a cross-sectionial view of the layer. The bottom figures show a cross section of the coupled structure, according to the code in the main text. A slight difference in skyrmion diameters can be observed in top and bottom layers for the asymmetric case.

Close modal

This section is accompanied by supplementary material mumax input file 7.txt.

Following the spirit of the last section, it is conceivable to go from spin structures that are isotropic in the thickness direction toward structures with even more variations in their profile in the third dimension. This opens a door to a world of new magnetic quasi-particles.19 One example, the chiral bobber,83 will be simulated in this section and an introduction to the simulation of their internal dynamics will be given. The applied methods can be generalized and easily utilized to simulate internal dynamics of other spin structures.

Chiral bobbers can be understood as skyrmion textures, which slowly decrease in diameter when going toward the bulk of the material. The decreasing size eventually terminates inside the material forming a so-called Bloch point. A point where all surrounding spins point toward or away from, leading to its description as a quasi-monopole. Bobbers can be found in simulations of various (usually thick layers of) materials, where they can appear in coexistence with regular skyrmion tubes. They can also be deliberately evoked with a useful trick by putting a uniformly magnetized layer with fixed spins or sufficiently high perpendicular magnetic anisotropy in proximity to a chiral layer. When setting up the chiral layer with a skyrmion as its initial state, energy minimization can immediately generate a bobber,

The result of above code is shown in Fig. 11(a). Next, we can study some internal dynamics. There are periodic motions that appear when driving the spin structure out of its equilibrium state. A frequency analysis often reveals specific resonant modes, that can vary with different material parameters and also depend on the resonating spin texture. The internal resonant modes of skyrmions, for example, can be categorized in breathing and gyrating motions.84 The latter is a periodic gyrating motion of the skyrmions’ center, mostly in a circular fashion, whereas the breathing motion is a periodic change in the skyrmion size.85 The breathing mode resonance finds practical use in a method proposed to stabilize the velocity of a skyrmion in a notched racetrack.86 The skyrmion’s velocity automatically adjusts itself to synchronize the time necessary to move between two notches with the inverse of its breathing mode frequency. By doing so, it ensures that the skyrmion size is minimized precisely when it passes through the constrictions, leading to enhanced stability of its motion. Although a combination of different internal modes often appears in concert, for the sake of this Tutorial, we will continue to focus on breathing modes. One possibility to induce these is by the application of a periodically changing external magnetic field in an out-of-plane configuration. Instead of applying a sinusoidal pulse with one frequency, we choose to apply a sinus cardinalis (sinc). The Fourier transformation of sinc ( 2 π f max t ) is a rectangular shaped spectrum ranging from f max to + f max. Hence, utilization of a sinc function is a feasible way of applying a range of frequencies simultaneously. Following code can be added to study the internal breathing dynamics:

To study high-frequency responses which are often in the GHz range, the Gilbert damping is set to a sufficiently low value. Here, we choose α = 0.01, however, smaller values might be needed for the investigation of even higher frequencies. The amplitude of the externally applied sinc pulse is chosen to be 0.5 mT, which is enough to induce internal breathing modes. The frequency is set as f max = 50 GHz, therefore, to resolve responses up to 50 GHz; the time step needs to be chosen according to the Nyquist–Shannon theorem τ s t e p < 1 / ( 2 50 GHz ) = 1 × 10 11 s. This value is set as the time step in which the results will be saved to the output directory. The file “table.txt” will then contain a table of the following form:

Since the breathing dynamics corresponds to a change in diameter, its effect can be directly quantified by changes in the average out-of-plane magnetization, as given by the fourth column in the table denoted by “mz.” A plot of the change in average out-of-plane magnetization δ m z ( t ) = m z ( 0 ) m z ( t ) is shown in Fig. 11(b). As can be seen, the system is driven out of equilibrium and subsequently relaxes in a periodic motion. The Fourier transform is presented in Fig. 11(c). It shows three specific peaks at 7.84, 44 and 50 GHz, whereas the last one might be shifted toward higher frequencies and would need further investigation by the application of pulses with higher frequency.

t (s)mx ()my ()mz ()
1.7089844e-09  1.1935763e-09  0.9543505 
1.0015…e-11 1.5597873e-09  1.1935763e-09  0.9543331 
2.0009…e-11 1.6411675e-09  8.6805557e-10  0.95431864 
3.0009…e-11 1.7903646e-09  1.3020833e-09  0.9543042 
… … … … 
t (s)mx ()my ()mz ()
1.7089844e-09  1.1935763e-09  0.9543505 
1.0015…e-11 1.5597873e-09  1.1935763e-09  0.9543331 
2.0009…e-11 1.6411675e-09  8.6805557e-10  0.95431864 
3.0009…e-11 1.7903646e-09  1.3020833e-09  0.9543042 
… … … … 

To get a better understanding how these breathing modes look in real space, one could apply the same Fourier analysis for each simulation cell rather than for all cells on average. For this, we have saved the necessary files with autosave(m , 1e-11). After Fourier analysis of each cell, we can pick the amplitude of each cell for a given frequency and make a color-coded plot of the result. This will visualize which cells have contributed to the periodic motion of a specific frequency. A cross section of this is shown for the main breathing mode at 7.84 GHz in the inset of Fig. 11(c), illustrating the localization at the domain wall region.

In conclusion, a decade after its initial release, this Tutorial paper serves as a comprehensive resource to simulating modern magnetic systems in mumax3. It offers an in-depth exploration and a collection of exemplary input files to guide researchers in simulating various magnetic material systems.

FIG. 11.

Chiral bobber breathing modes. (a) shows a chiral bobber, which basically comprises of a skyrmion texture at the surface with continuously decreasing diameter toward the bulk. (b) displays the breathing mode as a difference in the out-of-plane magnetization δ m z ( t ) after application of a sinus cardinalis (sinc)-shaped external magnetic field in the out-of-plane direction. The latter is shown in the inset. (c) shows the Fast Fourier Transformation (FFT) spectrum of the breathing dynamics. The resonance mode at 7.84 GHz is further analyzed in the inset. This shows the FFT results of individual cells on a cross section of the bobber. The color-code highlights which cells contribute the most to this specific resonance mode.

FIG. 11.

Chiral bobber breathing modes. (a) shows a chiral bobber, which basically comprises of a skyrmion texture at the surface with continuously decreasing diameter toward the bulk. (b) displays the breathing mode as a difference in the out-of-plane magnetization δ m z ( t ) after application of a sinus cardinalis (sinc)-shaped external magnetic field in the out-of-plane direction. The latter is shown in the inset. (c) shows the Fast Fourier Transformation (FFT) spectrum of the breathing dynamics. The resonance mode at 7.84 GHz is further analyzed in the inset. This shows the FFT results of individual cells on a cross section of the bobber. The color-code highlights which cells contribute the most to this specific resonance mode.

Close modal

Our discussions encompassed the integration of spin–orbit torques into mumax simulations either as a custom field or as a reparameterization of the Slonczewski spin-transfer torque. Additionally, we delved into the simulation of multilayer magnetic structures, spanning from synthetic antiferromagnets with or without explicitly added spacer layers to complex multilayer stacks. In certain cases, these stacks can be simplified to pseudo-2D simulations, while in others, they exhibit sophisticated spin structures with varying magnetization texture along the third dimension.

Our selection of these topics is not only driven by their current significance but also by the fact that the discussed solutions provide some universal recipes. We, therefore, anticipate that the examples presented here, particularly through the utilization of the custom fields functionality, will empower mumax users to contribute to the continuous advancement of magnetism research as new magnetic systems and phenomena emerge in the years to come.

The complete input files corresponding to all examples presented in the manuscript are made available in the supplementary material.

Part of the contents of this tutorial paper were developed for the mumax3 workshop https://mumax.ugent.be/mumax3-workshop/ presented in summer 2020 within the context of the IEEE online spintronics seminars, and for which we gratefully acknowledge support from the NVIDIA corporation with the donation of 2 Titan RTX GPU’s. J. L. is supported by the Research Foundation—Flanders (FWO) through a senior postdoctoral research fellowship (12W7622N). B.V.W. and J.J. have received funding from the European Union’s Horizon 2020 FET-Open program under Grant Agreement No. 861618 (SpinENGINE). P.B. acknowledges financial support from EPSRC (UK) through a DTG. Funding by the International Max Planck Research School for Science and Technology of Nano-Systems and the Max Planck Center for Quantum Materials is gratefully acknowledged. Finally, we wish to thank the (micro)magnetic community for the trust placed, and contributions made to mumax3 and hope that the content presented here will find meaningful applications and further propel the progress of the micromagnetic research field.

The authors have no conflicts to disclose.

Jonas J. Joos: Investigation (equal); Methodology (equal); Visualization (equal); Writing – original draft (lead); Writing – review & editing (equal). Pedram Bassirian: Investigation (equal); Methodology (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Pieter Gypens: Investigation (equal); Methodology (equal); Writing – review & editing (equal). Jeroen Mulkers: Conceptualization (equal); Investigation (equal); Methodology (equal); Software (equal); Writing – review & editing (equal). Kai Litzius: Investigation (equal); Methodology (equal); Supervision (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). Bartel Van Waeyenberge: Funding acquisition (equal); Supervision (equal); Writing – review & editing (equal). Jonathan Leliaert: Conceptualization (equal); Funding acquisition (equal); Investigation (equal); Methodology (equal); Supervision (lead); Writing – review & editing (equal).

All data presented in this paper can be generated using the input files in the supplementary material.

As first introduced in the supplementary material of Ref. 87 in the context of antiferromagnetic interfaces, mumax3 allows one to define custom scalar or vector quantities. These can subsequently be adopted in the effective field that is evaluated during simulated dynamics or energy minimization. Custom quantities are defined as shown in the code box below. Some standard arithmetic operations are available as well,

The latter operation, shifted, allows to shift a quantity across different cells in a given direction, specified by the three integers. This is a useful functionality to, e.g., calculate spatial derivatives by finite differences.

User-defined custom fields (in units of Tesla) are automatically accounted for during energy optimization and simulated magnetization dynamics,

Additionally, if the user has defined the corresponding energy densities, they can be added to the simulation to be evaluated efficiently on the GPU and the results can be added to the output table as follows. Note that mumax internally only uses the custom fields so that users that are not interested in outputting the energy densities are not obliged to define them,

When multiple arithmetic operations on custom fields are combined, there will in general be no CUDA kernel implemented for the combination. As a consequence, the different operations are ran consequently by CUDA. This is clearly not an optimal situation, but it is justified by the benefits of the possibility to have user-defined effective fields and energy densities.

1.
Q.
Shao
,
P.
Li
,
L.
Liu
,
H.
Yang
,
S.
Fukami
,
A.
Razavi
,
H.
Wu
,
K.
Wang
,
F.
Freimuth
,
Y.
Mokrousov
,
M. D.
Stiles
,
S.
Emori
,
A.
Hoffmann
,
J.
Åkerman
,
K.
Roy
,
J. P.
Wang
,
S. H.
Yang
,
K.
Garello
, and
W.
Zhang
, “
Roadmap of spin–orbit torques
,”
IEEE Trans. Magn.
57
,
1
39
(
2021
).
2.
J. C.
Slonczewski
, “
Current-driven excitation of magnetic multilayers
,”
J. Magn. Magn. Mater.
159
,
L1
L7
(
1996
).
3.
L.
Berger
, “
Emission of spin waves by a magnetic multilayer traversed by a current
,”
Phys. Rev. B
54
,
9353
9358
(
1996
).
4.
A.
Vansteenkiste
,
J.
Leliaert
,
M.
Dvornik
,
M.
Helsen
,
F.
Garcia-Sanchez
, and
B.
Van Waeyenberge
, “
The design and verification of MuMax3
,”
AIP Adv.
4
,
107133
(
2014
).
5.
P. A. M.
Dirac
, “The principles of quantum mechanics,” in International Series of Monographs on Physics (The Clarendon Press, 1930).
6.
G.
Dresselhaus
, “
Spin-orbit coupling effects in zinc blende structures
,”
Phys. Rev.
100
,
580
586
(
1955
).
7.
E. I.
Rashba
, “
Properties of semiconductors with an extremum loop. I. Cyclotron and combinational resonance in a magnetic field perpendicular to the plane of the loop
,”
Sov. Phys. Solid State
2
,
1224
(
1960
).
8.
Y. A.
Bychkov
and
E. I.
Rashba
, “
Oscillatory effects and the magnetic susceptibility of carriers in inversion layers
,”
J. Phys. C Solid State Phys.
17
,
6039
6045
(
1984
).
9.
V. M.
Edelstein
, “
Spin polarization of conduction electrons induced by electric current in two-dimensional asymmetric electron systems
,”
Solid State Commun.
73
,
233
235
(
1990
).
10.
M. I.
Dyakonov
and
V. I.
Perel
, “
Current-induced spin orientation of electrons in semiconductors
,”
Phys. Lett. A
35
,
459
460
(
1971
).
11.
Y. K.
Kato
,
R. C.
Myers
,
A. C.
Gossard
, and
D. D.
Awschalom
, “
Observation of the spin Hall effect in semiconductors
,”
Science
306
,
1910
1913
(
2004
).
12.
J.
Wunderlich
,
B.
Kaestner
,
J.
Sinova
, and
T.
Jungwirth
, “
Experimental observation of the spin-Hall effect in a two-dimensional spin-orbit coupled semiconductor system
,”
Phys. Rev. Lett.
94
,
047204
(
2005
).
13.
J. E.
Hirsch
, “
Spin Hall effect
,”
Phys. Rev. Lett.
83
,
1834
1837
(
1999
).
14.
S.
Zhang
, “
Spin Hall effect in the presence of spin diffusion
,”
Phys. Rev. Lett.
85
,
393
396
(
2000
).
15.
L. Q.
Liu
,
O. J.
Lee
,
T. J.
Gudmundsen
,
D. C.
Ralph
, and
R. A.
Buhrman
, “
Current-induced switching of perpendicularly magnetized magnetic layers using spin torque from the spin Hall effect
,”
Phys. Rev. Lett.
109
,
096602
(
2012
).
16.
L. Q.
Liu
,
C. F.
Pai
,
Y.
Li
,
H. W.
Tseng
,
D. C.
Ralph
, and
R. A.
Buhrman
, “
Spin-torque switching with the giant spin Hall effect of tantalum
,”
Science
336
,
555
558
(
2012
).
17.
P. M.
Haney
,
H. W.
Lee
,
K. J.
Lee
,
A.
Manchon
, and
M. D.
Stiles
, “
Current induced torques and interfacial spin-orbit coupling: Semiclassical modeling
,”
Phys. Rev. B
87
,
174411
(
2013
).
18.
J.
Kim
,
J.
Sinha
,
M.
Hayashi
,
M.
Yamanouchi
,
S.
Fukami
,
T.
Suzuki
,
S.
Mitani
, and
H.
Ohno
, “
Layer thickness dependence of the current-induced effective field vector in Ta|CoFeB|MgO
,”
Nat. Mater.
12
,
240
245
(
2013
).
19.
B.
Göbel
,
I.
Mertig
, and
O. A.
Tretiakov
, “
Beyond skyrmions: Review and perspectives of alternative magnetic quasiparticles
,”
Phys. Rep.
895
,
1
28
(
2021
).
20.
C.
Donnelly
,
M.
Guizar-Sicairos
,
V.
Scagnoli
,
S.
Gliga
,
M.
Holler
,
J.
Raabe
, and
L. J.
Heyderman
, “
Three-dimensional magnetization structures revealed with x-ray vector nanotomography
,”
Nature
547
,
328
331
(
2017
).
21.
P.
Vingelmann
and
F. H.
Fitzek
and,
NVIDIA
, “CUDA, Release: 10.2.89” (2020).
22.
See https://mumax.github.io/index.html for “mumax 3: GPU-accelerated Micromagnetism.”
23.
J.
Leliaert
,
M.
Dvornik
,
J.
Mulkers
,
J.
De Clercq
,
M. V.
Milošević
, and
B.
Van Waeyenberge
, “
Fast micromagnetic simulations on GPU–Recent advances made with mumax 3
,”
J. Phys. D: Appl. Phys.
51
,
123002
(
2018
).
24.
L.
Landau
and
E.
Lifshitz
, “
On the theory of the dispersion of magnetic permeability in ferromagnetic bodies
,”
Phys. Z. Sowjet.
8
,
153
(
1935
).
25.
T.
Gilbert
, “
A phenomenological theory of damping in ferromagnetic materials
,”
IEEE Trans. Magn.
40
,
3443
3449
(
2004
).
26.
W. F.
Brown
,
Micromagnetics
(
Interscience Publishers
,
1963
).
27.
J.
Leliaert
and
J.
Mulkers
, “
Tomorrow’s micromagnetic simulations
,”
J. Appl. Phys.
125
,
180901
(
2019
).
28.
W.
Heisenberg
, “
Zur theorie des ferromagnetismus
,”
Z. Phys.
49
,
619
636
(
1928
).
29.
W.
Döring
,
Ferromagnetism/Ferromagnetismus
(
Springer
,
1966
), pp. 341–437.
30.
J. B.
Collings
,
R.
Rama-Eiroa
,
R. M.
Otxoa
,
R. F. L.
Evans
, and
R. W.
Chantrell
, “
Generalized form of the magnetic anisotropy field in micromagnetic and atomistic spin models
,”
Phys. Rev. B
107
,
064413
(
2023
).
31.
I. E. J.
Dzyaloshinskii
, “
A thermodynamic theory of ‘weak’ ferromagnetism of antiferromagnetics
,”
J. Phys. Chem. Solids
4
,
241
255
(
1958
).
32.
T.
Moriya
, “
Anisotropic superexchange interaction and weak ferromagnetism
,”
Phys. Rev.
120
,
91
98
(
1960
).
33.
J. D.
Jackson
,
Classical Electrodynamics
(
John Wiley & Sons
,
1975
).
34.
D.
Kumar
and
A. O.
Adeyeye
, “
Techniques in micromagnetic simulation and analysis
,”
J. Phys. D: Appl. Phys.
50
,
343001
(
2017
).
35.
C.
Abert
, “
Micromagnetics and spintronics: Models and numerical methods
,”
Eur. Phys. J. B
92
,
120
(
2019
).
36.
A.
Lyberatos
,
D. V.
Berkov
, and
R. W.
Chantrell
, “
A method for the numerical simulation of the thermal magnetization fluctuations in micromagnetics
,”
J. Phys.: Condens. Matter
5
,
8911
8920
(
1993
).
37.
J.
Leliaert
,
J.
Mulkers
,
J.
De Clercq
,
A.
Coene
,
M.
Dvornik
, and
B.
Van Waeyenberge
, “
Adaptively time stepping the stochastic Landau-Lifshitz-Gilbert equation at nonzero temperature: Implementation and validation in MuMax 3
,”
AIP Adv.
7
,
125010
(
2017
).
38.
F.
Vanderveken
,
J.
Mulkers
,
J.
Leliaert
,
B.
Van Waeyenberge
,
B.
Sorée
,
O.
Zografos
,
F.
Ciubotaru
, and
C.
Adelmann
, “
Confined magnetoelastic waves in thin waveguides
,”
Phys. Rev. B
103
,
054439
(
2021
).
39.
F.
Vanderveken
,
J.
Mulkers
,
J.
Leliaert
,
B.
Van Waeyenberge
,
B.
Sorée
,
O.
Zografos
,
F.
Ciubotaru
, and
C.
Adelmann
, “
Finite difference magnetoelastic simulator
,”
Open Res. Eur.
1
,
35
(
2021
).
40.
L.
Exl
,
S.
Bance
,
F.
Reichel
,
T.
Schrefl
,
H.
Peter Stimming
, and
N. J.
Mauser
, “
LaBonte’s method revisited: An effective steepest descent method for micromagnetic energy minimization
,”
J. Appl. Phys.
115
,
17D118
(
2014
).
41.
In this article, we focus on cases in which the spin polarization (direction) of the current remains constant as the current traverses the magnetic layers. Practically, this means that we do not consider Zhang-Li-type STT’s,88 which describe, e.g., skyrmion motion through a nanotrack driven by a current flowing along the track.
42.
M.
Hayashi
,
J.
Kim
,
M.
Yamanouchi
, and
H.
Ohno
, “
Quantitative characterization of the spin-orbit torque using harmonic Hall voltage measurements
,”
Phys. Rev. B
89
,
144425
(
2014
).
43.
F.
Büttner
,
I.
Lemesh
,
M.
Schneider
,
B.
Pfau
,
C. M.
Günther
,
P.
Hessing
,
J.
Geilhufe
,
L.
Caretta
,
D.
Engel
,
B.
Krüger
,
J.
Viefhaus
,
S.
Eisebitt
, and
G. S. D.
Beach
, “
Field-free deterministic ultrafast creation of magnetic skyrmions by spin–orbit torques
,”
Nat. Nanotechnol.
12
,
1040
1044
(
2017
).
44.
J. C.
Slonczewski
, “
Currents and torques in metallic magnetic multilayers
,”
J. Magn. Magn. Mater.
247
,
324
338
(
2002
).
45.
J. C.
Slonczewski
, “
Currents, torques, and polarization factors in magnetic tunnel junctions
,”
Phys. Rev. B
71
,
024411
(
2005
).
46.
J.
Xiao
,
A.
Zangwill
, and
M. D.
Stiles
, “
Boltzmann test of Slonczewski’s theory of spin-transfer torque
,”
Phys. Rev. B
70
,
172405
(
2004
).
47.
A.
Meo
,
C. E.
Cronshaw
,
S.
Jenkins
,
A.
Lees
, and
R. F. L.
Evans
, “
Spin-transfer and spin-orbit torques in the Landau–Lifshitz–Gilbert equation
,”
J. Phys.: Condens. Matter
35
,
025801
(
2022
).
48.
M.
Donahue
and
D.
Porter
, “OOMMF User’s Guide 1.2a5” (2012).
49.
J. C.
Sankey
,
Y.-T.
Cui
,
J. Z.
Sun
,
J. C.
Slonczewski
,
R. A.
Buhrman
, and
D. C.
Ralph
, “
Measurement of the spin-transfer-torque vector in magnetic tunnel junctions
,”
Nat. Phys.
4
,
67
71
(
2008
).
50.
These terms mix up due to the damping α, but in this case the effect is neglegible as α was determined to be < 0.01 for the system described in Ref. 49.
51.
M. D.
Stiles
and
J.
Miltat
, “Spin-transfer torque and dynamics,” in Spin Dynamics in Confined Magnetic Structures III, edited by B. Hillebrands and A. Thiaville (Springer Berlin Heidelberg, Berlin, 2006), pp. 225–308.
52.
T.
Kimura
,
Y.
Otani
,
T.
Sato
,
S.
Takahashi
, and
S.
Maekawa
, “
Room-temperature reversible spin Hall effect
,”
Phys. Rev. Lett.
98
,
156601
(
2007
).
53.
A. V.
Khvalkovskiy
,
V.
Cros
,
D.
Apalkov
,
V.
Nikitin
,
M.
Krounbi
,
K. A.
Zvezdin
,
A.
Anane
,
J.
Grollier
, and
A.
Fert
, “
Matching domain-wall configuration and spin-orbit torques for efficient domain-wall motion
,”
Phys. Rev. B
87
,
020402
(
2013
).
54.
A.
Manchon
and
S.
Zhang
, “
Theory of nonequilibrium intrinsic spin torque in a single nanomagnet
,”
Phys. Rev. B
78
,
212405
(
2008
).
55.
N.
Nagaosa
and
Y.
Tokura
, “
Topological properties and dynamics of magnetic skyrmions
,”
Nat. Nanotechnol.
8
,
899
911
(
2013
).
56.
A.
Fert
,
N.
Reyren
, and
V.
Cros
, “
Magnetic skyrmions: Advances in physics and potential applications
,”
Nat. Rev. Mater.
2
,
1
15
(
2017
).
57.
C.
Back
,
V.
Cros
,
H.
Ebert
,
K.
Everschor-Sitte
,
A.
Fert
,
M.
Garst
,
T.
Ma
,
S.
Mankovsky
,
T. L.
Monchesky
,
M.
Mostovoy
,
N.
Nagaosa
,
S. S. P.
Parkin
,
C.
Pfleiderer
,
N.
Reyren
,
A.
Rosch
,
Y.
Taguchi
,
Y.
Tokura
,
K.
von Bergmann
, and
J.
Zang
, “
The 2020 skyrmionics roadmap
,”
J. Phys. D: Appl. Phys.
53
,
363001
(
2020
).
58.
K.
Everschor-Sitte
,
J.
Masell
,
R. M.
Reeve
, and
M.
Kläui
, “
Perspective: Magnetic skyrmions—Overview of recent progress in an active research field
,”
J. Appl. Phys.
124
,
240901
(
2018
).
59.
R.
Msiska
,
J.
Love
,
J.
Mulkers
,
J.
Leliaert
, and
K.
Everschor-Sitte
, “
Audio classification with skyrmion reservoirs
,”
Adv. Intell. Syst.
5
,
2200388
(
2023
).
60.
R.
Tomasello
,
E.
Martinez
,
R.
Zivieri
,
L.
Torres
,
M.
Carpentieri
, and
G.
Finocchio
, “
A strategy for the design of skyrmion racetrack memories
,”
Sci. Rep.
4
,
6784
(
2014
).
61.
A.
Fert
,
V.
Cros
, and
J.
Sampaio
, “
Skyrmions on the track
,”
Nat. Nanotechnol.
8
,
152
156
(
2013
).
62.
J.
Sampaio
,
V.
Cros
,
S.
Rohart
,
A.
Thiaville
, and
A.
Fert
, “
Nucleation, stability and current-induced motion of isolated magnetic skyrmions in nanostructures
,”
Nat. Nanotechnol.
8
,
839
844
(
2013
).
63.
J.
Mulkers
,
B.
Van Waeyenberge
, and
M. V.
Milošević
, “
Effects of spatially-engineered Dzyaloshinskii-Moriya interaction in ferromagnetic films
,”
Phys. Rev. B
95
,
144401
(
2017
).
64.
W.
Jiang
,
X.
Zhang
,
G.
Yu
,
W.
Zhang
,
X.
Wang
,
M.
Benjamin J.
,
J.
Pearson
,
X.
Cheng
,
O.
Heinonen
,
K. L.
Wang
,
Y.
Zhou
,
A.
Hoffmann
, and
S. E.
te Velthuis
, “
Direct observation of the skyrmion Hall effect
,”
Nat. Photonics
13
,
162
169
(
2017
).
65.
S.
Rohart
and
A.
Thiaville
, “
Skyrmion confinement in ultrathin film nanostructures in the presence of Dzyaloshinskii-Moriya interaction
,”
Phys. Rev. B
88
,
184422
(
2013
).
66.
X.
Zhang
,
G. P.
Zhao
,
H.
Fangohr
,
J. P.
Liu
,
W. X.
Xia
,
J.
Xia
, and
F. J.
Morvan
, “
Skyrmion-skyrmion and skyrmion-edge repulsions in skyrmion-based racetrack memory
,”
Sci. Rep.
5
,
7643
(
2015
).
67.
S.
Huang
,
C.
Zhou
,
G.
Chen
,
H.
Shen
,
A. K.
Schmid
,
K.
Liu
, and
Y.
Wu
, “
Stabilization and current-induced motion of antiskyrmion in the presence of anisotropic Dzyaloshinskii-Moriya interaction
,”
Phys. Rev. B
96
,
144412
(
2017
).
68.
K.-W.
Kim
,
K.-W.
Moon
,
N.
Kerber
,
J.
Nothhelfer
, and
K.
Everschor-Sitte
, “
Asymmetric skyrmion Hall effect in systems with a hybrid Dzyaloshinskii-Moriya interaction
,”
Phys. Rev. B
97
,
224427
(
2018
).
69.
B.
Göbel
,
A.
Mook
,
J.
Henk
,
I.
Mertig
, and
O. A.
Tretiakov
, “
Magnetic bimerons as skyrmion analogues in in-plane magnets
,”
Phys. Rev. B
99
,
060407
(
2019
).
70.
B.
Göbel
,
A. F.
Schäffer
,
J.
Berakdar
,
I.
Mertig
, and
S. S. P.
Parkin
, “
Electrical writing, deleting, reading, and moving of magnetic skyrmioniums in a racetrack device
,”
Sci. Rep.
9
,
12119
(
2019
).
71.
R.
Zarzuela
,
V. K.
Bharadwaj
,
K.-W.
Kim
,
J.
Sinova
, and
K.
Everschor-Sitte
, “
Stability and dynamics of in-plane skyrmions in collinear ferromagnets
,”
Phys. Rev. B
101
,
054405
(
2020
).
72.
X.
Zhang
,
J.
Xia
,
L.
Shen
,
M.
Ezawa
,
O. A.
Tretiakov
,
G.
Zhao
,
X.
Liu
, and
Y.
Zhou
, “
Static and dynamic properties of bimerons in a frustrated ferromagnetic monolayer
,”
Phys. Rev. B
101
,
144435
(
2020
).
73.
X.
Zhang
,
Y.
Zhou
, and
M.
Ezawa
, “
Magnetic bilayer-skyrmions without skyrmion Hall effect
,”
Nat. Commun.
7
,
10293
(
2016
).
74.
R.
Tomasello
,
V.
Puliafito
,
E.
Martinez
,
A.
Manchon
,
M.
Ricci
,
M.
Carpentieri
, and
G.
Finocchio
, “
Performance of synthetic antiferromagnetic racetrack memory: Domain wall versus skyrmion
,”
J. Phys. D: Appl. Phys.
50
,
325302
(
2017
).
75.
R.
Msiska
,
D. R.
Rodrigues
,
J.
Leliaert
, and
K.
Everschor-Sitte
, “
Nonzero skyrmion Hall effect in topologically trivial structures
,”
Phys. Rev. Appl.
17
,
064015
(
2022
).
76.
“MuMax-view: A web app to visualize OVF files. The GUI allows to scroll through multiple ovf files and visualises the magnetization either as cuboids or arrows. Among other visual options, it allows the user to choose between the mumax3 or red-white-blue gradient color scheme. It works on both 2D and 3D magnetizations and allows to cut through the structures to generate images of cross sections. It is freely available on https://mumax.ugent.be/mumax-view/.
77.
D.
Suess
,
S.
Koraltan
,
F.
Slanovc
,
F.
Bruckner
, and
C.
Abert
, “
Accurate finite-difference micromagnetics of magnets including RKKY interaction: Analytical solution and comparison to standard micromagnetic codes
,”
Phys. Rev. B
107
,
104424
(
2023
).
78.
S.
Woo
,
K.
Litzius
,
B.
Krüger
,
M.-Y.
Im
,
L.
Caretta
,
K.
Richter
,
M.
Mann
,
A.
Krone
,
R. M.
Reeve
,
M.
Weigand
,
P.
Agrawal
,
I.
Lemesh
,
M.-A.
Mawass
,
P.
Fischer
,
M.
Kläui
, and
G. S. D.
Beach
, “
Observation of room-temperature magnetic skyrmions and their current-driven dynamics in ultrathin metallic ferromagnets
,”
Nat. Mater.
15
,
501
506
(
2016
).
79.
W.
Legrand
,
J.-Y.
Chauleau
,
D.
Maccariello
,
N.
Reyren
,
S.
Collin
,
K.
Bouzehouane
,
N.
Jaouen
,
V.
Cros
, and
A.
Fert
, “
Hybrid chiral domain walls and skyrmions in magnetic multilayers
,”
Sci. Adv.
4
,
eaat0415
(
2018
).
80.
N.
Kerber
,
D.
Ksenzov
,
F.
Freimuth
,
F.
Capotondi
,
E.
Pedersoli
,
I.
Lopez-Quintas
,
B.
Seng
,
J.
Cramer
,
K.
Litzius
,
D.
Lacour
,
H.
Zabel
,
Y.
Mokrousov
,
M.
Kläui
, and
C.
Gutt
, “
Faster chiral versus collinear magnetic order recovery after optical excitation revealed by femtosecond XUV scattering
,”
Nat. Commun.
11
,
6304
(
2020
).
81.
K.
Halbach
, “
Design of permanent multipole magnets with oriented rare earth cobalt material
,”
Nucl. Instrum. Methods
169
,
1
10
(
1980
).
82.
P.
Bassirian
,
T.
Hesjedal
,
S. S. P.
Parkin
, and
K.
Litzius
, “
Breathing mode dynamics of coupled three-dimensional chiral bobbers
,”
APL Mater.
10
,
101107
(
2022
).
83.
F. N.
Rybakov
,
A. B.
Borisov
,
S.
Blügel
, and
N. S.
Kiselev
, “
New type of stable particle like states in chiral magnets
,”
Phys. Rev. Lett.
115
,
117201
(
2015
).
84.
M.
Mochizuki
, “
Spin-wave modes and their intense excitation effects in skyrmion crystals
,”
Phys. Rev. Lett.
108
,
017601
(
2012
).
85.
J.-V.
Kim
,
F.
Garcia-Sanchez
,
J. A.
Sampaio
,
C.
Moreau-Luchaire
,
V.
Cros
, and
A.
Fert
, “
Breathing modes of confined skyrmions in ultrathin magnetic dots
,”
Phys. Rev. B
90
,
064410
(
2014
).
86.
J.
Leliaert
,
P.
Gypens
,
B.
Van Waeyenberge
,
J.
Mulkers
, and
M.
Milošević
, “
Coupling of the skyrmion velocity to its breathing mode in periodically notched nanotracks
,”
J. Phys. D: Appl. Phys.
52
,
024003
(
2019
).
87.
J.
De Clercq
,
J.
Leliaert
, and
B.
Van Waeyenberge
, “
Modelling compensated antiferromagnetic interfaces with MuMax 3
,”
J. Phys. D: Appl. Phys.
50
,
425002
(
2017
).
88.
S.
Zhang
and
Z.
Li
, “
Roles of nonequilibrium conduction electrons on the magnetization dynamics of ferromagnets
,”
Phys. Rev. Lett.
93
,
127204
(
2004
).

Supplementary Material