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.

## I. INTRODUCTION

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] $ \xd7 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.

### A. Motivation and outline

This Tutorial paper celebrates the 10th anniversary of the first release (in 2013) of mumax3. Mumax is a powerful NVIDIA CUDA^{21} 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.

## II. MICROMAGNETISM IN MUMAX3

The foundations of micromagnetic theory originate from the 1930s^{24} 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}

^{24,25}

All energy terms (or derived effective fields) are user-specified by means of phenomenological parameters. Intrinsic materials properties include the exchange stiffness ( $ A ex$^{28,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 neglected^{30}]. 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}

### A. Dynamics and statics: Run, relax, and minimize

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 \u2212 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.

## III. SPIN TORQUES IN MICROMAGNETISM

### A. Effective field

^{88}and spin–orbit torques

^{41}can be described by the Landau–Lifshitz–Gilbert equation, provided that a contribution is added to the effective field,

^{42,43}

^{2,3}

### B. Torque

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.

### C. Parameters for spin-transfer torques

^{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 by

^{44–47}

^{,}

^{48}(Object Oriented MicroMagnetic Framework). Herein, the parameters $\Lambda $, $d$, and $ \u03f5 \u2032$ 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 $ \u03f5 \u2032$ was introduced to allow asymmetry between the properties of the fixed and free layers.

^{46}Measuring the value of $ \u03f5 \u2032$ for experimental systems is far from trivial, but is described in detail in Ref. 49, where a value of $ \u03f5 \u2032/\u03f5\u22480.3$ is reported for the system under study, as determined by the ratio between the Slonczewski torque in the perpendicular $( m\xd7 p)$ and in-plane $[ m\xd7( p\xd7 m)]$ directions.

^{50}

The value of $\Lambda $ grasps deviations from the idealized scenario ( $\Lambda =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 $\Lambda $ *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.

### D. Parameters for spin–orbit torques

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.

^{52}( $ j s$), i.e., $ j s= | \alpha H j c |/e$, hence substituting $P j z$ for $ | \alpha H j c |$ in Eq. (13) immediately yields the corresponding damping-like torque magnitude resulting from the spin-Hall effect,

^{53}

^{54}They applied the Boltzmann equation to the stationary eigenstates of the single-electron Rashba Hamiltonian featuring spin–orbit (with coupling strength $ \alpha R$) and exchange (with coupling strength $J$) interactions, leading to

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}

## IV. SPIN–ORBIT TORQUE-DRIVEN SKYRMION PROPULSION IN A FERROMAGNETIC STRIP

Skyrmions are localized metastable magnetization textures with a non-trivial topology^{55} 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 Dind^{63} 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), $\alpha $ (specified by alpha, dimensionless), and $ A ex$ (specified by Aex, in J/m). A Co strip of dimensions $256\xd764\xd71 nm 3$ is simulated with a skyrmion in the initial magnetization,

Two equivalent implementations of spin–orbit torque are illustrated.

### A. Implementation 1: Effective field

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,

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( \alpha H)>0$ for Pt] and current (see also Fig. 1). A phenomenological value of $\u2212$2.0 is used for $\xi $ [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$.

### B. Implementation 2: SOT as additional torque

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, $ \tau SOT$ has the same mathematical form as Slonczewski STT, $ \tau 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 $| \alpha 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.

### C. Result

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.

## V. SPIN–ORBIT TORQUE-DRIVEN SKYRMION PROPULSION IN A SYNTHETIC ANTIFERROMAGNET

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.

### A. Implementation 1: Touching layers

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,

### B. Implementation 2: Separated 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.

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.

### C. Result

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=\u22121$,

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=\u2212128$ 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.

## VI. COMPLEX MULTILAYER STRUCTURES

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

### A. Effective medium approach

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)] $ \xd7 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 \u2032$, thus being able to reduce all individual layers into a single one with thickness $N t m \u2032$.

^{78}and results in the rather intuitive scaling,

*not*scaled and remains unchanged in this set of scaling laws. Since spin–orbit torques (with $ b J=\xi =0$ in the case of pure damping-like spin-Hall torque) follow the form

^{78}From Eq. (14),

### B. Comparison between effective 2D and full 3D simulation

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 $\u2212z$ 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 $\u2212z$ directions,

### C. Dealing with highly symmetric initial states

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,

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.

## VII. STRAY FIELD EFFECTS

### A. Thick magnetic strips and bulk materials

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\xd7250\xd7200 nm 3$ material dimension discretized by $64\xd764\xd750$ 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.

### B. Low DMI multilayers with many layers

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.

### C. Skyrmion Halbach array

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\xd7150\xd721 nm 3$ gridsize discretized by $1\xd71\xd71 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.

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.

## VIII. CHIRAL BOBBERS AND INTERNAL DYNAMICS

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\pi f maxt)$ is a rectangular shaped spectrum ranging from $\u2212 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 $\alpha =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 $ \tau s t e p<1/(2\u22c550 GHz)=1\xd7 10 \u2212 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 $\delta m z(t)= m z(0)\u2212 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 () . |
---|---|---|---|

0 | 1.7089844e-09 | $\u2212$1.1935763e-09 | $\u2212$0.9543505 |

1.0015…e-11 | 1.5597873e-09 | $\u2212$1.1935763e-09 | $\u2212$0.9543331 |

2.0009…e-11 | 1.6411675e-09 | $\u2212$8.6805557e-10 | $\u2212$0.95431864 |

3.0009…e-11 | 1.7903646e-09 | $\u2212$1.3020833e-09 | $\u2212$0.9543042 |

… | … | … | … |

t (s) . | mx () . | my () . | mz () . |
---|---|---|---|

0 | 1.7089844e-09 | $\u2212$1.1935763e-09 | $\u2212$0.9543505 |

1.0015…e-11 | 1.5597873e-09 | $\u2212$1.1935763e-09 | $\u2212$0.9543331 |

2.0009…e-11 | 1.6411675e-09 | $\u2212$8.6805557e-10 | $\u2212$0.95431864 |

3.0009…e-11 | 1.7903646e-09 | $\u2212$1.3020833e-09 | $\u2212$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.

## IX. CONCLUSIONS

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.

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.

## SUPPLEMENTARY MATERIAL

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

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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).

## DATA AVAILABILITY

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

### APPENDIX A: CUSTOM FIELDS FUNCTIONALITY

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.

## REFERENCES

*International Series of Monographs on Physics*(The Clarendon Press, 1930).

^{88}which describe, e.g., skyrmion motion through a nanotrack driven by a current flowing

*along*the track.

*Spin Dynamics in Confined Magnetic Structures III*, edited by B. Hillebrands and A. Thiaville (Springer Berlin Heidelberg, Berlin, 2006), pp. 225–308.