Micromagnetics and spintronics: Models and numerical methods

Computational micromagnetics has become an indispensable tool for the theoretical investigation of magnetic structures. Classical micromagnetics has been successfully applied to a wide range of applications including magnetic storage media, magnetic sensors, permanent magnets and more. The recent advent of spintronics devices has lead to various extensions to the micromagnetic model in order to account for spin-transport effects. This article aims to give an overview over the analytical micromagnetic model as well as its numerical implementation. The main focus is put on the integration of spin-transport effects with classical micromagnetics.


Introduction
The micromagnetic model has proven to be a reliable tool for the theoretical description of magnetization processes on the micron scale. In contrast to purely quantum mechanical theories, such as density functional theory, micromagnetics does not account for distinct magnetic spins nor nondeterministic effects due to collapse of the wave function. However, micromagnetics integrates quantum mechanical effects that are essential to ferromagnetism, like the exchange interaction, with a classical continuous field description of the magnetization in the sense of expectation values. The main assumption of this model is that the organizing forces in the magnetic material are strong enough to keep the magnetization in parallel on a characteristic length scale λ well above the lattice constant a where S i{j and x i{j are distinct spins and their positions, respectively. For a homogeneous density of spins, the discrete distribution of magnetic moments S i is well approximated by a continuous vector density M pxq such that renders approximately true for arbitrary volumes Ω of the size λˆλˆλ and larger with 1 Ω being the indicator function of Ω. The continuous magnetization M pxq a e-mail: claas.abert@univie.ac.at has a constant norm due to the homogeneous density of spins and can thus be written in terms of a unit vector field mpxq M pxq " M s mpxq with |mpxq| " 1, where M s is the spontaneous magnetization. In the case of zero temperature, which is often considered for micromagnetic modeling, M s is the saturation magnetization which is a material constant. While m and M have to be strictly distinguished, both are referred to as magnetization throughout this work for the sake of simplicity. Due to the combination of classical field theory with quantum mechanical effects, micromagnetics is often referred to as semiclassical theory. Opposed to the macroscopic Maxwell equations, the micromagnetic model resolves the structure of magnetic domains and domain walls. This enables accurate hysteresis computations of macroscopic magnets, since hysteresis itself is the direct result of field-induced domain generation and annihilation. While static hysteresis computations are very important for the development of novel permanent magnets [1], another application area for micromagnetics is the description of magnetization dynamics on the micron scale. The time and space-resolved description of magnetic switching processes and domain-wall movements are essential for the development of novel storage and sensing technologies such as magnetoresistive random-access memory (MRAM) [2], sensors for read heads [3], and angle sensors [4]. Besides the manipulation of the magnetization with external fields, the interaction of spin polarized electric currents with the magnetization plays an increasing role for novel devices. Several extensions to the classical micromagnetic theory have been proposed in order to account for these spintronics effects.
A lot of articles and books have been published on analytical [5][6][7] as well as numerical [1,8,9] micromagnetics. This review article is supposed to serve two purposes. First, it is meant to give a comprehensive yet compact overview over the micromagnetic theory, both on an analytical as well as a numerical level. The article describes the most commonly used discretization strategies, namely the finite-difference method (FDM) and the finite-element method (FEM), and discusses advantages, disadvantages and pitfalls in their implementation. The second and main purpose of this article, however, is to give an overview over existing models for spin transport in the context of micromagnetics. This article reviews the applicability and the limits of these models and discusses their discretization.

Energetics of a ferromagnet
The total energy of a ferromagnetic system with respect to its magnetization is composed by a number of contributions depending on the properties of the respective material. While some of these contributions, like the demagnetization energy and the Zeeman energy, can be described by classical magnetostatics, other contributions like the exchange energy and the magnetocrystalline anisotropy energy have a quantum mechanical origin. This section aims to give an overview over typical energy contributions and their representation in the micromagnetic model.

Zeeman energy
The energy of a ferromagnetic body highly depends on the external field H zee . The corresponding energy contribution is often referred to as Zeeman energy. According to classical electromagnetics the Zeeman energy of a magnetic body Ω m is given by E zee "´µ 0 ż Ωm M s m¨H zee dx (4) with µ 0 being the vacuum permeability.

Exchange energy
The characteristic property of ferromagnetic materials is its remanent magnetization, i.e. even for a vanishing external field, a ferromagnetic system can have a nonvanishing macroscopic magnetization. In a system where the spins are coupled by their dipole-dipole interaction only, the net magnetization always vanishes for a vanishing external field as known from classical electrodynamics [10]. However, in ferromagnetic materials, the spins are subject to the so-called exchange interaction. For two localized spins, this quantum mechanical effect energetically favors a parallel over an antiparallel spin alignment. The origin of this energy contribution can be attributed to the Coulomb energy of the respective two-particle system, typically consisting of two electrons. Depending on the spin-alignment, the two-particle wave function is either symmetric or anti-symmetric leading to higher expectation value of the distance, and thus a lower expectation value of the Coulomb energy, in case of an parallel alignment. The classical description of the exchange interaction is given by the Heisenberg model. Details on its derivation can be found in any textbook on quantummechanics, e.g. [11]. The Heisenberg formulation of the exchange energy of two unit spins s i and s j is defined as with J being the so-called exchange integral. With respect to the continuous magnetization field mpxq, the exchange energy associated with all couplings of a single spin site x is given by where the index i runs over all exchange coupled spin sites at positions x`∆x i , and J i denotes the exchange integral with the respective spin. Expression (7) is obtained by application of the unit-vector identity pn 1´n2 q 2 " 2´2n 1¨n2 and performing Taylor expansion of lowest order. The exchange integral J i highly depends on the distance of the spin sites. Hence, significant contributions to the exchange energy are only provided by nearby spins, usually next neighbors. The transition from the discrete Heisenberg model to a continuous expression for the total exchange energy is done by integrating (7) while considering a regular spin lattice, i.e. a regular spacing of the spin sites x as well as identical J i and ∆x i for each site.
In the most general form this procedure yields (8) where the coefficients of the matrix A jk depend on the crystal structure and the resulting exchange couplings of the spins in the magnetic body. The term C results from the integration of the constant part of (7) and is usually neglected since it does not depend on m and thus only gives a constant offset to the energy without changing the physics of the system. The matrix A can always be diagonalized by a proper choice of coordinate system [12] which yields For cubic and isotropic lattice structures, the exchange coupling constants A j simplify further to the scalar exchange constant A which results in the typical micromagnetic expression for the exchange energy (10) where p∇mq 2 " ř i,j pBm i {Bx j q 2 is to be understood as a Frobenius inner product. Although derived for localized spins and isotropic lattice structures, this energy expression turns out to accurately describe a large number of materials including band magnets and anisotropic materials [13]. This is explained by the fact, that (10) exactly represents the lowest order phenomenological energy expression that penalizes inhomogeneous magnetization configurations.

Demagnetization energy
The demagnetization energy accounts for the dipoledipole interaction of a magnetic system. This energy contribution, that is also referred to as magnetostatic energy or stray-field energy, owes its name to the fact that magnetic systems energetically favor macroscopically demagnetized states if they are subject to dipole-dipole interaction only. For a continuous magnetization M " M s m the demagnetization energy can be derived from classical electromagnetics. Assuming a vanishing electric current j e " 0, Maxwell's macroscopic equations reduce to ∇¨B " 0 (11) ∇ˆH dem " 0 (12) where the magnetic flux B can be written in terms of the magnetic field H dem and the magnetization M B " µ 0 pH dem`M q.
According to (12), the magnetic field H dem is conservative and thus has a scalar potential H dem "´∇u. With these definitions (11)- (12) can be reduced to the single equation which is solved in the whole space R 3 . Assuming a localized magnetization configuration, the boundary conditions are given in an asymptotical fashion by upxq " Op1{|x|q for |x| Ñ 8 (15) which is referred to as open boundary condition since the potential is required to drop to zero at infinity. The defining equation (14) is often transformed in Poisson's equation However, in contrast to the original equation (14), the divergence on the right-hand side of (16) may become singular at the boundary of the magnetic material in the case of a localized magnetization. In this case, (16) is well defined only in a distributional sense. The potential u can be expressed in terms of an integral equation by considering the well-known fundamental solution to the Laplacian that naturally fulfills the required open boundary conditions [10] upxq "´1 4π For a localized magnetization with sharp boundaries, this solution, like (16), suffers from a singular divergence at the boundary of the magnetic body. However, this singularity is integrable as demonstrated in the following. Consider a finite magnet with |M pxq| " M s for x P Ω m surrounded by a thin transition shell Ω t where the magnetization continuously decays to zero, see Figure 1. In this case, the integration over R 3 in (17) can be reduced to an integration over Ω m Y Ω t since the magnetization, and thus the integrand of (17), vanishes outside the magnet. Further, the integral is split into integration over Ω m and Ω t , and the integral over Ω t is transformed with Green's theorem (18) where ds 1 denotes the areal measure to x 1 and n is an outward-pointing normal vector. In order to obtain the potential for an ideal magnet with a sharp transition of the magnetic region Ω m to the air region, we consider the limit of a vanishing transition region Ω t Ñ 0. In this case, the right-hand side of (18) reduces to the boundary integral. Furthermore, the boundary integral vanishes for the outer boundary of Ω t , because of a vanishing magnetization M . The inner boundary, however, coincides with the boundary of the magnetic region BΩ m except for its orientation. A complete integral form for the magnetic scalar potential of an ideal localized magnet in region Ω m accordingly reads In analogy to the integral equation for the electric field, the terms ρ "´∇¨M and σ " M¨n are often referred to as magnetic volume charges and magnetic surfaces charges, respectively. An alternative integral expression for the potential is obtained by applying Green's theorem to (19) upxq " 1 4π Starting from this formulation, the demagnetization field H dem can be expressed as a convolution H dem pxq "´∇upxq " ż ΩmÑ px´x 1 qM px 1 q dx 1 (21) with the so-called demagnetization tensorÑ given bỹ According to classical electrodynamics, the energy connected to the demagnetization field is given by where the factor 1{2 accounts for the quadratic dependence of the energy on the magnetization M . The competition of the demagnetization energy and the exchange energy leads to the formation of magnetic domains. A graphic example for this effect is the magnetic vortex structure depicted in Figure 2. In order to minimize surface charges σ, that contribute to the demagnetization energy, the magnetization field aligns parallel with edges and surfaces, which explains the curl-like configuration. The exchange energy favors a parallel alignment of the magnetization which leads to the creation of four distinct, almost homogeneously magnetized, triangular domains, each aligned with one of the square's edges. A perfect in-plane curl configuration, which completely avoids surface charges, is very unfavorable with respect to the exchange energy because it leads to a singularity in the center of the curl. In order to reduce the exchange energy, the magnetization rotates out-of-plane in a distinct area around the center of the vortex, called the vortex core.

Crystalline anisotropy energy
Another important contribution to the total free energy of a magnet is the anisotropy energy that favors the parallel alignment of the magnetization to certain axes referred to as easy axes. The origin of this energy lies in the spin-orbit coupling either due to an anisotropic crystal structure or due to lattice deformation at material interfaces [13]. Depending on the symmetry of these anisotropies, the respective material will exhibit one or multiple easy axes. These axes are undirected and thus the energy does not depend on the sign of the magnetization For the simplest case of a single easy axis, the anisotropy energy is given by where e u is a unit vector parallel to the easy axis and K u1 and K u2 are the scalar anisotropy constants. This phenomenological expression is obtained by symmetry considerations. For a uniaxial anisotropy, the energy may depend only on the angle between the magnetization and easy axis m¨e u . Furthermore, only even powers in m¨e u are considered in order to fulfill condition (24). Uniaxial anisotropy typically occurs in materials with a hexagonal or tetragonal crystal structure, e.g. cobalt.
Materials with cubic lattice symmetry such as iron, which has a body-centered cubic structure, exhibit three easy axes e i which are pairwise orthogonal Like for the uniaxial anisotropy, the expression for the cubic anisotropy energy is developed as series in magnetization components along the easy axes up to sixth order where m i " e i¨m is the projection of the magnetization m on the anisotropy axis e i . Only contributions compatible with the symmetry condition (24) are considered. Moreover, the resulting expression is required to be constant under permutation of magnetization components m i in order to have a cubic symmetry.
While the magnetization prefers to align parallel to the respective axes in the case of positive anisotropy constants K u1 , K u2 , K c1 , K c2 , the magnetization avoids a parallel configuration for negative anisotropy constants. In the case of uniaxial anisotropy, this leads to an easyplane anisotropy. In the case of cubic anisotropy, this leads to four easy axes as is the case for nickel which has a face-centered cubic structure.
Both, equations (25) and (27) hold for magnetic anisotropies in bulk material. If magnetic anisotropy is caused at material interfaces, either due to lattice deformation or due to the electric band structure, the energy depends on the magnetization configuration m at this interface only. The energy for such a surface anisotropy is obtained by similar expressions as (25) and (27). However, instead of integrating over the magnetic volume Ω m the integration in this case has to be carried out over the respective interface BΩ m only.
Although being derived only phenomenologically, the expressions (25) and (27) have proven to describe anisotropy effects with high accuracy. For many application it is even sufficient to consider the lowest order contributions only and setting the higher order constants K u2 and K c2 to zero.

Antisymmetric exchange energy
As discovered by Dzyaloshinskii [14] and Moriya [15], neighboring spins can be subject to an antisymmetric exchange interaction in addition to the regular exchange interaction discussed in Section 2.2. This effect, that is often referred to as Dzyaloshinskii-Moriya interaction (DMI), is caused by the spin-orbit coupling in certain material systems. The general antisymmetric exchange energy of two spins s i and s j is given as where the vector d ij depends on the symmetry of the system. A typical system that gives rise to DMI is a magnetic layer with an interface to a heavy-metal layer. In this case, the antisymmetric exchange between two neighboring magnetic spins near the interface is mediated by a single atom in the heavy metal layer and the vector d ij is given as where d is a scalar coupling constant, ∆x " ∆x{|∆x| is a unit vector pointing from spin site i to spin site j, and e d is the interface normal. The transition to continuum theory is done similar to the exchange interaction. Namely, in a first step, the energy of the couplings for a single spin site is expressed in terms of the continuous magnetization field m and the magnetization at the neighboring site mpx`∆xq is expanded in powers of ∆x to the lowest order where the vector identity paˆbq¨pcˆdq " pa¨cqpb¨dqṕ a¨dqpb¨cq was used and the summation is carried out over the coupled neighboring spins. Performing integration and assuming isotropic coupling d i as well as an isotropic lattice spacing ∆x i , similar to the exchange interaction in Section 2.2, yields the continuous expression for the total antisymmetric exchange energy for interface DMI. The scalar coupling constant D i depends on the coupling constants d i as well as the relative positions ∆x i . Another class of materials exhibiting DMI are magnetic bulk materials lacking inversion symmetry [16,17]. For these materials, the coupling vector d ij is given as which results in the following energy E x for the couplings of a single spin site x Again, assuming isotropic coupling d i and lattice spacing ∆x i results in the continuous formulation for the energy with the coupling constant D b depending on the atomistic coupling constants d i and the lattice spacing ∆x i . Besides the prominent interface and bulk DMI, further antisymmetric exchange couplings are defined by Lifshitz invariants [18,19]. The antisymmetric exchange counteracts the regular exchange energy that favors homogeneous magnetization configurations and penalizes domain walls. It gives rise to a magnetization configuration called skyrmion,  see Figure 3. A skyrmion is characterized by a continuous rotation of the magnetization field on any lateral axis crossing the center. It has topological charge meaning that the skyrmion configuration is not continuously transformable into the homogeneous ferromagnetic state.

Interlayer-exchange energy
The magnetic layers of a multilayer structure may be exchange coupled even when separated by a nonmagnetic spacer layer. This coupling, which was first proposed by Ruderman and Kittel [20], is mediated by the conducting electrons of the nonmagnetic layer. The coupling constant A shows oscillatory behavior with respect to the thickness of the spacer layer, i.e. depending on its thickness, the coupling of the magnetic layers may be either ferromagnetic or antiferromagnetic. This effect was described in a more generalized theory by Kasuya [21] and Yosida [22] and is referred to as Ruderman-Kittel-Kasuya-Yosida (RKKY) interaction.
In the continuous approximation, the interaction is assumed to couple the interface between one magnetic layer and the spacer layer Γ 1 and the interface between the other magnetic layer and the spacer layer Γ 2 . These interfaces are assumed to have equal size. Integration of the Heisenberg interaction (5) yields the continuous expression for the interlayer-exchange energy A mpxq¨mrP pxqs ds (38) with A being the exchange constant whose sign and strength depend on the thickness of the spacer layer and P : Γ 1 Ñ Γ 2 being an isomorphism that maps any point on Γ 1 to its nearest point on Γ 2 . The RKKY interaction is often exploited in order to build so-called synthetic antiferromagnets, see Figure 4. For this purpose, the thickness of the spacer layer is chosen such that A is negative which results in an antiferromagnetic coupling of the magnetic layers. Synthetic antiferromagnets are important for applications due to their stability and lack of strayfield.

Other energy contributions and effects
While this work will focus on the energy contributions introduced above, there are numerous additional energy contributions and other effects that may play important roles in certain systems [5]. For instance, the effect of magnetostriction coupled the mechanical properties of magnetic materials to its magnetization configuration [23,24]. If magnetic systems are subject to charge currents, eddy currents [25,26] and the Oersted field [27] need to be considered.
Another important area of research is finite-temperature micromagnetics. Various approaches have been proposed in order account for temperature effects in micromagnetics, among them Langevin dynamics [28][29][30] and the Landau-Lifshitz-Bloch equation [31][32][33]. However, this comprehensive topic is out of the scope of this article.

Static micromagnetics
Static micromagnetics is the theory of stable magnetization configurations and hence a valuable model for the investigation of material properties such as hysteresis. The prerequisite for a stable magnetization configuration is a minimum of the total free energy E of the system with respect to its magnetization m. In order to be a valid micromagnetic solution, the solution m is further required to fulfill the micromagnetic unit-sphere constraint min Epmq subject to |mpxq| " 1.
Since the solution variable m is a continuous vector field, variational calculus is applied in order to solve for an energetic minimum. A necessary condition for a minimum is a functional differential δE that vanishes for arbitrary test functions v P V m with V m being the function space of the magnetization m An alternative formulation for this condition can be stated in terms of the functional derivative δE{δm that is defined as where the function space V 0 m Ă V m includes only functions of V m that vanish on the boundary vpBΩ m q " 0. That said, depending on the considered energy E, the differential δE as defined in (40) in general differs from the left-hand side of (41) by a boundary integral, i.e.
This means that the knowledge of the functional derivative (41) is not sufficient in order to solve the minimization problem (39). In general, additional boundary conditions defined by Bpmq have to be considered. All variational considerations so far do not account for the unit-sphere constraint |m| " 1. This constraint can be incorporated by a Lagrange multiplier technique where the modified functional E λ pm, λ, µq, given by is minimized with respect to both the magnetization m and the Lagrange multiplier fields λ and µ implementing the constraint on the volume and surface respectively. The solution to this minimization problem is again obtained by variational calculus where the variations of the solution variables v m , v λ and v µ can be treated separately where V λ and V µ are appropriate function spaces for the variation of λ and µ, respectively. Expanding (44) considering the definition (43) yields Since (48) has to vanish for arbitrary v m P V m , it also vanishes for test functions with vanishing boundary values v m P V 0 m . Hence the functional derivative of the energy has to fulfill δE δm "´2λm (49) in Ω m which is required to hold for arbitrary λ P V λ . This condition is satisfied if and only if δE{δm is parallel to m and hence which is exactly Brown's condition [5]. Testing (48) with functions that are defined on the boundary only vpΩ m zBΩ m q " 0 and considering the surface Lagrange multiplier µ in the same manner as above yields the additional boundary condition Moreover, inserting (43) into (45) yields " 0 (54) which is required to hold for arbitrary v λ P V λ and thus represents the micromagnetic constraint |m| 2 " 1. Due to the interface Lagrange multiplier µ, this constraint is further specifically enforced on the boundary by (46). In the following, the functional derivatives and boundary conditions for the energy contributions introduced in Section 2 will be discussed in detail.

Zeeman energy
The energy differential δE{δm for the Zeeman energy is obtained by variation of (4) which yields The variation does not give rise to any additional boundary integral. Hence, the derivative and boundary term B for the Zeeman energy are given by

Exchange energy
The differential δE ex for the exchange energy is derived from (10) resulting in Here, integration by parts is performed in order to eliminate spatial derivatives of the test functions v m . The resulting volume integral is of the same form as the integral in (41) which enables the identification of the functional derivative δE ex {δm. However, this necessary step also gives rise to a surface integral and thus to a boundary term B. The resulting derivative and boundary term for the exchange energy read where (62) can be simplified to δE ex {δm "´2A∆m if A is assumed constant throughout the magnetic region Ω m .

Demagnetization energy
The differential for the demagnetization energy is obtained similarly to the differential for the Zeeman energy. A decisive difference to the Zeeman energy, however, is the linear relation of the demagnetization field H dem pM q to the magnetization M . The variation of the magnetization therefore leads to an additional factor of 2 which results in the differential Consequently the derivative and boundary term for the demagnetization energy are given by

Anisotropy energy
For the uniaxial anisotropy (25) the derivative and boundary terms are given by δE aniu δm "´2K u1 e u pe u¨m q´4K u2 e u pe u¨m q 3 (68) B " 0 (69) and for the cubic anisotropy (27) the respective terms read For interface anisotropy contributions, the variational derivative δE{δm obviously vanishes and the influence of the energy contribution reduces to the boundary term with pxq being the respective areal energy density.

Antisymmetric exchange energy
Similar to the exchange energy, the variation of the antisymmetric exchange energy (33) needs to be transformed by partial integration in order to eliminate spatial derivatives of the test functions v m The resulting variational derivative and the boundary term for the interface DMI energy read B "´D i pe dˆn qˆm.
For the antisymmetric bulk exchange (37) the differential is given by which leads to the following variational derivative and boundary term

Energy minimization with multiple contributions
In order to minimize the total energy of a system subject to multiple energy contributions both Brown's condition (50) and the boundary condition (51) have to be fulfilled for the composite energy functional. Namely, if a system is subject to the exchange energy (10) and the demagnetization energy (23), the respective conditions for an energy minimum read with (84) being the "classical" micromagnetic boundary condition. Spatial derivatives of the magnetization m are always orthogonal to the magnetization due to the micromagnetic unit-sphere constraint. Hence, the boundary condition (84) is usually simplified to Bm{Bn " 0.
If the system is additionally subject to the antisymmetric exchange (33), both Brown's condition (83) and the boundary condition (84) are supplemented with the respective contributions. The resulting boundary condition reads where the cross product mˆB was again neglected by orthogonality arguments. Depending on the considered energy contributions, this boundary condition is supplemented by additional terms. Hence, adding energy contributions does not add additional boundary conditions, but changes the single boundary condition instead.

Dynamic micromagnetics
In Micromagnetics, magnetization dynamics is described by the Landau-Lifshitz (LL) equation that was originally proposed in [34]. This equation describes the spatially resolved motion of the magnetization in an effective field. Due to problems with the dissipative term, an alternative formulation was derived by Gilbert [35,36]. Both formulations are completely equivalent under proper parameter transformation. However, for the purpose of distinction, the latter is usually referred to as Landau-Lifshitz-Gilbert (LLG) equation or, alternatively, as Gilbert or implicit form of the LL equation.
The LLG can be derived by means of classical Lagrangian mechanics by the choice of an appropriate action S. Due to the micromagnetic unit-sphere constraint |m| " 1, the magnetization field m can be described by means of spherical coordinates mpxq "˜s inrθpxqs cosrφpxqs sinrθpxqs sinrφpxqs cosrθpxqs¸(

86)
with the polar angle θ and the azimuthal angle φ. For the sake of readability we omit the spatial dependence of fields in the following. According to Hamilton's principle, the temporal evolution of any field is given as the path with stationary action δS " 0 with the action S defined as where L is the so-called Lagrangian which, in turn, is given by with T being the kinetic energy density and V being the potential energy density. The potential energy density V is naturally given by the free energy E " ş V dx whose contributions are introduced in Section 2. However, the choice of the kinetic energy T is not immediately clear. Due to the unit-sphere constraint, the motion of the magnetization is restricted to rotations. Hence, it seems reasonable to assume a kinetic energy similar to that of a rotating rigid body with I being an inertial tensor and Ω being the angular velocity vector. In the rigid-body picture, the magnetization in a certain point mpxq is represented by a cylindrical stick with one end fixed at the coordinate origin, see Figure 5. We introduce the fixed-body frame with coordinate axes e 1 1 , e 1 2 , e 1 3 as shown in Figure 5 and mark vectors with a prime whose coordinates are expressed in terms of this new basis. In the fixed-body frame, the magnetization is trivially given as Due to the cylindrical symmetry of the rigid-body representation of the magnetization, it is clear that the inertia tensor is diagonal in the fixed-body frame and thus reads The angular velocity vector in the fixed-body frame is given as where the spherical coordinates θ and φ are complemented by the angle ψ that describes the rotation of magnetization's stick representation around its symmetry axis. In order to derive the LLG from the general kinetic energy density (89), two assumptions are required. The first assumption is that of vanishing moments of inertia I 1 " I 2 " 0. This assumption is reasonable since the magnetization stick has no mass in a classical sense. Hence, the rotation of a magnetic moment in an external field is expected to stop instantaneously if the external field is switched off rapidly. With this assumption, the kinetic energy (89) reduces to The second assumption is that the angular momentum of the rotation around the symmetric axis L 3 is connected to the saturation magnetization by the relation where γ 3 is the electron's gyromagnetic ratio. This relation is reasonable since the saturation magnetization takes the place of the magnetic moment in the continuous theory of micromagnetics and the magnetic moment is generated by the spin, i.e. the angular momentum connected to the symmetry axis. Inserting (92) into (93) and further using the relation (94) yields the following expression for the variation of the integrated kinetic energy with respect to the azimuthal angle φ Applying the same procedure to the variation with respect to the polar angle θ yields The angular velocity (92) can be simplified by setting ψ " 0. This angle describes the rotation of the magnetization's stick representation around its symmetry axis. Hence, this assumption does not introduce any losses of generality [37].
This simplification leads to the following relations for the time derivatives of magnetization coordinates and their respective variations Inserting into the kinetic-energy variations (98) and which, in the fixed-body frame, can be summarized in the vector valued variation (105) in the fixed-body frame where the magnetization is given as m " p0, 0, 1q, see (90). In order to compute the variation of the action δS, the variation of the kinetic energy (105) has to be complemented by the variation of the potential energy which is given as where the boundary term B is the same as introduced in Section 3 and hence depends on the particular choice of energy contributions. Putting together the variation of the kinetic energy (105) and the potential energy (108) results in the variation of the action where B t m¨m vanishes due to the micromagnetic unitsphere constraint that requires any derivative of the magnetization to be perpendicular on m. The equation of motion (111) describes the magnetization dynamics without energy losses. However, realistic systems are expected to lose magnetic energy by conversion to e.g. phonons, eddy currents [5]. In the framework of Lagrangian mechanics, dissipative processes are described by a Rayleigh function D. The time evolution of the magnetization m subject to the dissipative function D is then given as The Rayleigh function D is usually chosen to be proportional to the square of the time derivative of the variable of motion. Choosing D " αM s {p2γ e qpB t mq 2 in the case of magnetization dynamics yields where α ě 0 is a dimensionless damping parameter. Inserting the variation of the action (109) and once more considering variations v m P V 0 m that vanish on the boundary BΩ m results in the equation of motion By introducing the effective field defined as this equation can be turned into the well-known Gilbert form of the LLG with γ " µ 0 γ e « 2.2128ˆ10 5 m{As being the reduced gyromagnetic ratio. This equation of motion is completed by the boundary condition mˆB " 0 which is obtained by varying v m on the boundary for (109) and by considering the same cross product with m that is applied to obtain (111). Note that this boundary condition resembles the static micromagnetic boundary condition (51). From (109) it is clear that this boundary condition has to hold at all times. This semi-implicit formulation introduced by Gilbert can be transformed into an explicit form by inserting the complete right-hand side of (117) into B t m on the righthand side of (117). Applying basic vector algebra and considering m¨B t m " 0 and m¨m " 1 yields B t m "´γ 1`α 2 mˆH eff´α γ 1`α 2 mˆpmˆH eff q (118) which, apart from the definition of the parameters γ and α, equals the original equation introduced by Landau and Lifshitz. While the presented derivation of the LLG is not rigorous, it should be noted that the required assumptions, namely vanishing moments of inertia I 1 " I 2 " 0 and the connection of the remaining moment of inertia with the saturation magnetization M s " γ e I 3 Ω 3 , are physically reasonable. The strict application of variational calculus not only yields the LLG but also its boundary conditions depending on the contributions to the effective field H eff .
A more detailed investigation on the LLG as derived from a Lagrangian is presented in [38] where it is also shown that the kinetic contribution to the Lagrangian (93) is equivalent to the assumption introduced in the original work by Gilbert [35] and earlier by Döring [6]. A full quantum mechanical description of a spin subject to exchange interaction, anisotropy and Zeeman field is given in [39] where the LLG equation is also obtained in a limit case. Figure 6 illustrates the damped precessional motion of the magnetization m in an effective field H eff as described by the LLG. As noted in Section 1, the main assumption of micromagnetics is the constant modulus of the magnetization field m. This property is conserved by the LLG as can be seen by considering the time derivative of the squared magnetization

Properties of the Landau-Lifshitz-Gilbert equation
Inserting the right-hand side of (118) yields B t |m| 2 " 0 and hence also B t |m| " 0. In classical micromagnetics, the energy connected to the magnetization, as defined by the sum of the energy contributions introduced in Section 2, may only change due to external fields varying in time or due to the energy dissipation modeled by the damping term. In the case of an effective field that does not explicitly depend on the time t the time derivative of the energy is given by Inserting (118) for B t m yields The value of the integral is zero in the case of an energy minimum, see Brown's condition (50), and positive otherwise. Hence, for a positive damping constant α ą 0, the energy of a magnetic system is a non-increasing function in time. In this case the LLG is said to have Lyapunov structure [40,41]. In the special case of no damping α " 0 the right-hand side of (124) vanishes and the LLG has Hamiltonian structure, i.e. it preserves the energy.

Spintronics in micromagnetics
The term spintronics summarizes all effects caused by the interaction of electrons with solid state devices due to their spin rather than their charge. For magnetic systems, this particularly covers the origin of spinpolarized currents and their impact on the magnetization configuration. The term spintronics was coined in the 1980s when the giant magnetoresistance (GMR) was discovered by Fert and coworkers [42] and Grünberg and coworkers [43]. Exploiting the spin of electrons, in addition to their charge, adds extra degrees of freedom and allows for the development of novel devices especially in the areas of storage and sensing technology.
In the semiclassical picture of micromagnetics, the spin polarization of an electric current is described by a threedimensional vector field p. If a polarized electric current passes a magnetic region, it exerts a torque on the magnetization. This so-called spin torque, similar to the torque generated by a magnetic field, can be split into a fieldlike contribution T field and a dampinglike contribution T damp , see Figure 7. The fieldlike torque has the same form as Fig. 7. Spin-torque contributions for neglectable damping α ! 1 with respect to the reference polarization p. The fieldlike torque T field leads to a precessional motion of the magnetization m around p. The dampinglike torque T damp leads to a direct relaxation of m toward p. the torque generated by a regular effective-field contribution, i.e. it leads to a damped precessional motion as described in Section 4. Since the Gilbert damping α is usually small α ! 1, the magnetization dynamics caused by the fieldlike torque are dominated by the precessional part. In contrast, the dynamics caused by the dampinglike torque are dominated by the direct rotation of the magnetization toward the polarization and accompanied by a small precessional contribution. Depending on the origin of the polarized current, the torque is either referred to as spin-transfer torque or spin-orbit torque.

Spin-transfer torque in multilayers
A typical device that exploits spin-transfer torque, consists of two magnetic layers separated by a nonmagnetic spacer layer and sandwiched with two nonmagnetic leads, see Figure 8. If passed by an electric current, the conducting electrons are subject to scattering processes depending on the spin configuration of the conducting electrons. Even if the applied current has a net spin polarization of zero, these spin-dependent scattering processes lead to a non-vanishing spin-polarization distribution across the multilayer. In particular, the interfaces between magnetic and non-magnetic regions act as scattering sites due to the rapid transition of magnetization. A simplified illustration of the scattering processes for different magnetization configurations of the multilayer, i.e. antiparallel and parallel, is given in Figure 8.
In the case of an antiparallel configuration, a scattering process takes place at the first interface that is passed by the conducting electrons, see Figure 8a. Due to this scattering, the FM1 layer acts as a spin filter and the majority of the conducting electrons that reach the FM2 layer carry the polarization of FM1. The spin-transfer torque will cause the switching of the FM2 magnetization if exceeding a critical strength. In addition, a scattering process at the first FM2 interface will reflect electrons with the polarization of FM1 leading to a stabilization of the FM1 magnetization configuration.
For a parallel configuration, the same scattering as for the antiparallel case appears at the first interface of FM1, see Figure 8b. This leads to a spin polarization parallel to the magnetization configuration of FM1 in the spacer layer between FM1 and FM2. However, if the spacer layer has sufficient thickness, the spin polarization reduces due to spin-flip events. At the first interface of FM2 the recovered electrons with antiparallel polarization to the magnetization of FM2 are scattered back to FM1. The scattered electrons exert a torque on the magnetization of FM1 and can switch it if exceeding a critical strength. The magnetization of FM2 on the other hand is stabilized by the electrons polarized by the FM1 layer. Possible applications of this torque mechanism, that was first investigated in by Slonczewski [44], Berger [45] and Waintal et al. [46], are the spin-transfer torque magnetoresistive random access memory (STT MRAM) [2,47] and spin torque oscillators (STO) [48,49]. A comprehensive theoretical overview over spin-transfer torque is given in a work by Ralph and Stiles [50]. A very popular model for the description of the magnetization dynamics in spin-transfer-torque devices is the model proposed by Slonczewski [51]. This model uses the macrospin approach, where the magnetic region subject to spin torque, also referred to as free layer, is described by a single spin m. The current is assumed to be polarized by another magnetic layer, referred to as polarizing layer, whose magnetization is described by the vector p. The motion of the free-layer magnetization m is described by the extended LLG where the torque T consists of a dampinglike and fieldlike contribution T " T damp`T field . According to the model of Slonczewski these contributions are given by where the dimensionless functions η damp and η field describe the angular dependence of the torque strength with ϑ being the angle between m and p. By comparing the torque contributions (127) and (128) with the effective-field term in the LLG (126), the torque can be expressed by means of an effective field contribution H T given by Inserting into the LLG (126) and transforming the LLG into the explicit form (118) yields where the vector identity mˆrmˆpmˆpqs "´mˆp was used. From this formulation it is clear that both the dampinglike torque T damp and the fieldlike torque T field contribute to the precessional motion as well as the dampinglike motion. However, this intermixing highly depends on the Gilbert damping α. In the limit case of vanishing α, the LLG simplifies to where the fieldlike torque exclusively contributes to the precessional motion and the dampinglike torque exclusively contributes to the dampinglike motion. This limit demonstrates the unique feature of the dampinglike torque to facilitate a dampinglike motion independently from the Gilbert damping α. In the original work of Slonczewski, the expression is derived as angular dependence of the torque for symmetric systems, i.e. systems with two identical magnetic layers. This expression is valid for both the dampinglike and the fieldlike torque, but in general requires a different set of model parameters P and Γ. The dimensionless parameters P and Γ depend on geometry and materials of the complete system and describe the polarization strength and the angular asymmetry of the STT, respectively. A more general expression for the angular dependence is introduced in [52] as This expression accounts for asymmetric devices with two different magnetic layers. Again, the free parameters q`, q´, A and B are dimensionless and depend on the geometry and material composition of the complete stack. The model of Slonczewski is often used to describe STT devices with one hard-magnetic layer acting as spin polarizer and one soft magnetic layer that is subject to the spin torque induced by the hard magnetic layer. For these devices the magnetization dynamics in the hard magnetic layer, referred to as reference layer or pinned layer, are neglectable and LLG is solved for the soft magnetic layer, referred to as free layer, only [53]. However, the model can also be used to describe the bidirectional coupling of the magnetization configuration in both magnetic layers [54].
The macrospin approach as proposed in the original work of Slonczewski is accurate for small systems below the single-domain limit. Systems of this size are dominated by the exchange interaction and hence act much like a single spin. With growing size, other energy contributions such as the demagnetization energy gain influence leading to the generation of magnetic domains, which renders the macrospin approximation useless. For thin film structures with large lateral dimensions but small thicknesses below the exchange length, the generalization of the macrospin model is straightforward. In this case, the magnetization m and the polarization p in (127) and (128) are functions of the lateral position x in the multilayer stack and the tilting angle ϑ is computed accordingly, see Figure 9. The torque contributions (127) and (128) are usually applied as volume terms with the polarization p assumed to be independent of the perpendicular position in the free layer.
However, the spin-transfer torque is considered to be a surface effect rather than a volume effect. For free-layer thicknesses below the exchange length, the treatment as volume effect does not affect the torque in a qualitative fashion, since the magnetization can be assumed constant across the free layer in this case. Yet the strength of the torque needs to be scaled with the reciprocal free-layer thickness 1{d in order to account for the surface nature of the effect.
While the generalization from a macro-spin model to a spatially resolved model allows for the description of various multilayer devices, the model of Slonczewski has a number of shortcomings. The presented generalization neglects lateral diffusion of the spins which might introduce inaccuracies for strongly inhomogeneous magnetization configurations. Moreover, the treatment as volume term is only justified for free-layer thicknesses below the exchange length. Another disadvantage of this model are the free parameters introduced in (132) and (133) which depend on the geometry and material parameters of the complete system in a nontrivial fashion. A comprehensive overview over Slonczewski-like models is given in a work by Berkov and Miltat [55].

Spin-transfer torque in continuous media
Spin-transfer torque cannot only be exploited in magnetic multilayer structures, but also in continuous single phase magnets. In these systems, magnetic domains take over the role of the distinct layers in multilayer stacks, i.e. they act as spin polarizer. While the spin torque in multilayers acts on the surfaces of the magnetic layers to the spacer layer, the spin torque in single phase magnets acts in regions of high magnetization gradients, i.e. domain walls. A simple picture for the origin of spin torque in single phase magnets is given in Figure 10. While the conducting electrons pass the magnetic region, they pick up the polarization from the local magnetization and carry this polarization in the direction of the electron flow where they exert a torque. This mechanism can be used to move domain walls and complete domain structures with electric currents. In contrast to field induced domain-wall motion, the spin-torque moves magnetic domain walls always in the direction of electron motion regardless of the nature of the wall, e.g. head-to-head, tail-to-tail. This property is exploited, for example, by the magnetic racetrack memory proposed by Parkin et al. [56].
An established model for the description of spin torque in continuous magnets is the model proposed by Zhang and Li [57]. In this model, the torque contribution T to the LLG is given as where ξ describes the degree of nonadiabacity according to [57] and b is given as with β being the dimensionless polarization rate of the conducting electrons, µ B being the Bohr magneton and The Zhang-Li model delivers reasonable results for the description of current driven domain-wall motion. However, the description of spin-torque is purely local since the torque only depends on first derivatives of the magnetization. This means that the diffusion of spin polarization is completely neglected. Consequently, the model is not suited for the description of spin torque in multilayers, since this requires the transport of spin across a nonmagnetic spacer layer. Moreover, the lack of diffusion also introduces inaccuracies in systems with highly inhomogeneous magnetization configurations.

Spin-diffusion
Both, the model of Slonczewski introduced in Section 5.1 and the model of Zhang and Li introduced in Section 5.2 are applicable only for specific material systems and magnetization configurations. Also, both models neglect the diffusion of the spin polarization to some extent. A more general approach to spin torque considers the torque generated by a vector field s referred to as spin accumulation where J denotes the coupling strength of the spin accumulation s and the magnetization m. The torque definition leads to the extended LLG The spin accumulation spxq describes the deviation of the conducting electron's polarization compared to the equilibrium configuration at vanishing charge current j e " 0. That said, by definition s is zero if no current is applied to the system. Several variations of the spin-diffusion model have been proposed for the computation of the spin accumulation s. According to [58] the dynamics of s are given by where τ sf denotes the spin-flip relaxation time which is a material parameter. While the LLG (137) is defined in the magnetic region Ω m only, the spin accumulation s is generated also in nonmagnetic regions such as the leads or the spacer layer of an STT device and hence has to be solved in the complete sample region Ω. Consequently, (138) is potentially defined in composite media and thus all material parameters, such as J and τ sf , may vary spatially. The matrix-valuedj s in (138) denotes the spin current defined bỹ with u being the electric potential, µ B being the Bohr magneton and e being the elementary charge. The variables C 0 , D 0 and β are material parameters. C 0 is connected to the electric conductivity σ by σ " 2C 0 and D 0 denotes the material's diffusion constant. β is a dimensionless constant that denotes the rate of polarized conducting electrons in magnetic materials. If the distribution of the electric potential u is known, the coupled system (137) and (138) can be solved in order to simultaneously resolve the dynamics of the magnetization mptq and the spin diffusion sptq. However, the distribution of the electric potential u might not be known upfront. Specifically, the charge current j e in the diffusion model is defined as which suggests a strong coupling of the electric potential u with m and s. If the charge current distribution j e is known, the spin-diffusion dynamics (138) can be solved by inserting (140) into (139) via the gradient of the electric potential ∇u, which results in the spin-current definitioñ (141) Inserting into (138) directly yields the spin-accumulation dynamics for a given magnetization configuration m. The resulting magnetization dynamics due to the spin torque can be resolved by simultaneous solution of (138) and the LLG (137).
However, in most realistic systems, the spin accumulation relaxes two orders of magnitude faster than the magnetization configuration [57]. If the quantity of interest is the magnetization dynamics rather than the spin-accumulation dynamics, this difference in time scales can be exploited in order to simplify the model. Namely, the spin-accumulation can be assumed to instantaneously relax when the magnetization changes. In this case, the spin-accumulation does no longer explicitly depend on the time t, but only on the magnetization m. The defining equation for the spin accumulation spmq is derived by setting B t s " 0 in (138) which results in Inserting the definition of the spin current (139) yields a linear partial differential equation of second order in s. Typical solutions for the spin accumulation s in STT devices as well as domain walls are depicted in Figure 11. By application of this simplified model, the treatment of the spin-accumulation s in the context of dynamical micromagnetics becomes similar to the treatment of effective-field contributions. Instead of performing a coupled time integration on both the magnetization m and the spin accumulation s as required by (138), the spin accumulation is defined by the magnetization m only. The previous methods require the knowledge of the charge-current distribution j e for the computation of the spin accumulation s. In a regular shaped sample with a homogeneous conductivity σ, the charge-current density may be assumed constant in a first-order approximation. However, in a magnetic system subject to spin-polarized currents, Ohm's law gives only one contribution to the total conductivity which may also depend on the magnetization configuration. In order to accurately account for magnetization dependent resistance effects, neither the electric potential u nor the charge current j e should be prescribed in the sample. Instead, these entities should be solution variables like the spin accumulation s. In order to set up such a self-consistent model, the source equation (142) has to be complemented by a source equation for the charge current, which is naturally given by the continuity equation Inserting the current definitions (139) and (140) in the source equations (142) and (143) yields a system of linear partial differential equations that can be solved for the solution pair ps, uq with the magnetization m being an input to the system. In this self-consistent model, instead of prescribing the charge current j e or potential u in the complete sample, boundary conditions are used to define the potential or current inflow on specific interfaces.

Boundary conditions for the spin-diffusion model
Since the partial differential equations for the solution of the spin-diffusion model are of second order in both solution variables u and s, boundary conditions are required in order to find a unique solution. The boundary conditions of the electric potential u directly correspond to the voltage or charge current applied to the system through electric contacts. A typical multilayer system with contact regions Γ 1 and Γ 2 is depicted in Figure 12. By choosing either Dirichlet or Neumann conditions for u on a part of the boundary, a constant electric potential or chargecurrent inflow can be prescribed on this part, respectively. For example, in order to simulate a constant potential u 0 at contact Γ 1 , a Dirichlet condition is applied directly to the potential u u " u 0 on Γ 1 .
Applying a constant current inflow j 0 on contact Γ 2 is achieved by applying the Neumann condition j e¨n ",´2 where n denotes the outward pointing normal to Γ 2 . In order to complete the set of boundary conditions for u, all parts of the sample's boundary which are not used as contacts are treated with homogeneous Neumann conditions The spin accumulation s is treated with homogeneous Neumann conditions on the complete boundary which is equivalent to a no-flux condition on the spin currentj s¨n " 0 for systems as depicted in Figure 12, where the contacts belong to the boundary of the nonmagnetic region Ω n . This equivalence is obtained by multiplying (141) with the boundary normal n and inserting the Neumann condition which yields the boundary flux This spin-current flux is nonzero only at boundaries with both nonvanishing charge-current flux j e¨n and nonvanishing magnetization m.
A vanishing charge-carrier flux usually implies a vanishing spin-current flux since the spin is transported by the charge carriers. This makes the no-flux condition on the spin current a reasonable choice for parts of the boundary that do not serve as electric contacts. For electric contacts, the no-flux condition is reasonable only if the thickness of the respective nonmagnetic regions exceeds the spin-flip relaxation length. In this case, the polarization of the current and thus the spin flux can be assumed zero at the contact, see e.g. Figure 11a. If this is not the case, the contact should be treated with a Robin condition instead where the exponential decay of the spin accumulation in the nonmagnetic lead region is taken into account. While the boundary conditions on the electric potential u are relevant only for the self-consistent treatment of u and s, the boundary conditions for the spin accumulation s also hold for the simplified spin-diffusion models with prescribed electric potential or charge current, respectively.

Spin-orbit torque
Several extensions to the spin-diffusion model introduced in Section 5.3 have been proposed in order to account for further spintronics effects. An important class of effects describes the conversion of charge currents to spin currents and vice versa due to spin-orbit coupling. The origin for this conversion is the polarization dependent deflection of the conducting electrons either due to material impurities [59,60] or due to intrinsic asymmetries of the material [61,62]. Depending on the direction of the current conversion, this spin-orbit effect is either referred to as spin-Hall effect or inverse spin-Hall effect. The spin-Hall effect describes the conversion of charge currents into spin currents, see Figure 13a, while the conversion of spin currents into charge currents is referred to as inverse spin-Hall effect, see 13b. Incorporating these effects into the spindiffusion model is done by extending the original current definitions (140) and (139) according to Dyakonov [63]. The extended current definitions j 1 e andj 1 s are defined in terms of the original current definitions and read where index notation was used. Here, ijk is the Levi-Cevita tensor and θ SH is the dimensionless spin-Hall angle. Inserting j 1 e andj 1 s into the source equations (143) and (142) yields a self-consistent spin-diffusion model including spin-Hall effects. Typically, the spin-Hall effects are exploited in multilayer structures with heavy-metal layers which are subject to spin-orbit coupling and neighboring magnetic layers where the spin-polarized currents interact with the magnetization configuration. Similar to the considerations in the original spin-diffusion model, equations (151) and (152) together with the respective source equations are solved in the complete structure including magnetic and nonmagnetic regions, using spatially varying material parameters in order to account for the different material properties.

Material parameters in the spin-diffusion model
The spin-diffusion model introduces a number of material parameters to the set of parameters required by classical micromagnetics. Depending on the exact formulation of the spin-diffusion model, different sets of parameters are used. The spin-flip relaxation time τ sf and the exchange coupling of the spin-accumulation and the magnetization J are often specified in terms of the characteristic length scales λ sf and λ J defined by Alternatively, the exchange strength J may be quantified by the characteristic time τ J " {J. While classical micromagnetics is often used to describe single-phase materials, the spin-diffusion model is usually used to solve the spin accumulation in composite systems. In order to account for such systems, that expose different material properties in different regions, all material parameters in the governing equations are scalar fields rather than constants. It should be noted that none of the material parameters, both in classical micromagnetics as well as in the spin-diffusion model, are subject to spatial differentiation. Hence, the material-parameter fields may comprise jumps across material interfaces without compromising the mathematical formulation of the model.

Spin dephasing
In addition to the spin-flip relaxation and the exchange coupling in the original spin-diffusion model (138), an additional term for the description of spin-dephasing was proposed in several works [64][65][66] B t s "´∇¨j s´s τ sf´sˆm τ J´mˆp where τ φ is the spin-dephasing time. This extended spindiffusion model can be solved similar to the original model either in nonequilibrium or equilibrium and either for prescribed charge current or self-consistently.

Valet-Fert model
An alternative to the spin-diffusion model introduced in Section 5.3 was introduced by Valet and Fert in [67]. The originally one-dimensional model for collinear magnetization configurations was generalized to three dimensions and noncollinear configurations in [68]. Similar to the Zhang-Levy-Fert model introduced in Section 5.3, the Valet-Fert model defines charge and spin currents J e and J s as well as an electric potential φ and a spin potential φ s which corresponds to the spin accumulation s. However, in contrast to the Zhang-Levy-Fert model, the spin current and spin potential are assumed to be collinear to the magnetization in the ferromagnetic regions Ω m With these simplified assumptions, the Valet-Fert model defines the currents in the magnetic region Ω m as with ρ˚being connected to the electric conductivity and β 2 being a dimensionless polarization parameter. In the nonmagnetic region Ω n , the magnetization m as well as the polarization parameter β 2 vanishes and the currents are defined as with (160) being Ohm's law. The source equations for the currents are given as where λ is a spatially varying material parameter that denotes the characteristic spin-flip relaxation length. In contrast to the Zhang-Levy-Fert model, the Valet-Fert model does not require the electric potential φ and the spin potential φ s to be continuous across interfaces. Instead, these potentials are subject to a set of welldefined jump conditions. The charge current j e is continuous everywhere which includes interfaces. The spin current j s is continuous within magnetic/nonmagnetic layers. Furthermore, its longitudinal component is continuous across magnetic/nonmagnetic interfaces n¨Jś " " m¨Js ı¨n where the '´' superscript corresponds to values in the magnetic layer and the '`' superscript corresponds to values in the nonmagnetic layer, see Figure 14. Both, the charge potential φ and the spin potential φ s may have jumps at magnetic-nonmagnetic interfaces defined by where the interface properties rb , γ sf and g Ö denote the resistivity, the spin-flip probability and the spin-mixing conductance of the interface, respectively. While the interface resistivity rb and the splin-flip probability γ sf have bulk counterparts in the Zhang-Levy-Fert model, namely C´1 0 and τ´1 sf , the spin-mixing conductance g Ö is unique to the Valet-Fert model. It describes the interface resistivity for electrons polarized perpendicular to the magnetization in the ferromagnetic layer and contributes significantly to the overall resistivity of the interface. Moreover, the spin-mixing conductance g Ö plays an important role for the description of spin torque in the context of the Valet-Fert model. Since the spin potential φ s is collinear in the magnetic regions by definition of the model, see (156), it cannot exert a torque on the magnetization m. Thus, torque can only be generated at the interface between magnetic and nonmagnetic regions, where the spin potential φ s can have components perpendicular to the magnetization. In general, the spinmixing conductance is assumed to be complex-valued with the real and imaginary part describing the strength of the dampinglike and fieldlike torque, respectively. This said, the Valet-Fert model does not predict the ratio of these torque contributions in contrast to the spin-diffusion model introduced in Section 5.3. As for the simplified model by Slonczewski, this ratio is an input parameter to the method.

Connecting the spintronics models
Various micromagnetic models for the description of spin torque and other spintronics effects are described in the preceding section. While the spin-diffusion model introduced in Section 5.3 covers a multitude of spintronics effects, specialized models like the model by Slonczewski, see Section 5.1, were developed for very specific purposes.

Slonczewski model
The Slonczewski model describes the spin torque in multilayers in terms of a macrospin approach. The characteristic properties of this model are the angular dependencies η damp pϑq and η field pϑq of the torques defined in (127) and (128). While the original angular dependency proposed by Slonczewski (132) is predicted to be valid for structures with two similar magnetic layers, the more general form (133) is expected to work also for asymmetric structures. In order to compare the model of Slonczewski with the spin-diffusion model, the torque for a homogeneously magnetized free layer with varying tilting angle ϑ is computed with the spin-diffusion model. The angular dependence of this torque is extracted and compared to the Slonczewski model in Figure 15. The asymmetric system used for this comparison has two magnetic layers with thicknesses 3 nm and 5 nm and typical material parameters as used in [69]. The symmetric system has similar material parameters, but uses the same layer thickness of 3 nm for both magnetic layers. The symmetric structure is well fitted by the original expression for the angular expression, see Figure 15a. For the asymmetric structure, i.e. a structure with different free-layer and pinned-layer thicknesses, the more general expression (133) is required in order obtain an accurate fit, see Figure 15b. This proves agreement of the two models in the application scope of the Slonczewski model and consequently the superiority of the spin-diffusion model which is able to accurately describe further spin-transport effects and devices.

Zhang-Li model
Like the model of Slonczewski, the model of Zhang and Li introduced in Section 5.2 can be perceived as a special case of the spin-diffusion model. Setting D 0 " 0 in (141) and inserting in (142) βµ B e pj e¨∇ qm`s τ sf´J mˆs " 0.
Multiplying with m and mˆm respectively and eliminating mˆpmˆsq terms results in the torque which has the form as the model of Zhang and Li (134) with the degree of nonadiabacity being defined as This means, that the model of Zhang and Li exactly reproduces the spin-diffusion model for vanishing diffusion D 0 . Neglecting the spin diffusion restricts the model of Zhang and Li to the description of local torque phenomena, i.e. torque due to local magnetization gradients such as domain walls.

Valet-Fert model
In contrast to the Slonczewski and Zhang-Li models, the Valet-Fert model is closely linked to the spin-diffusion model introduced in Section 5.3. Within the nonmagnetic layers the models are completely similar and in the magnetic regions the equations have common terms. A major difference of the models is the role of interfaces that have distinct properties such as the resistivity rb , the spinflip probability γ sf and the spin-mixing conductance g Ö in the Valet-Fert model. In contrast, the Zhang-Levy-Fert model introduced in Section 5.3 solely relies on bulk material properties. For the bulk properties, there is a straightforward mapping of the solution variables and material parameters. Using (156) and (157) and assuming a constant magnetization in the ferromagnetic layers ∇m " 0 yields for the Valet-Fert model. With spatially resolved material parameters ρ and β 2 , these current definitions can be used for the complete sample region Ω. Assuming the following relations, (172) and (173) can be identified with the definitions of the Zhang-Levy-Fert model (140) and (139).
where it should be noted that instead of the two polarization parameters β and β 1 of the Zhang-Levy-Fert model, the Valet-Fert model introduces only a single parameter β 2 . Comparison of the source equations for the spin current in the different models (142) and (163) yields the following additional parameter mappings With these parameter mappings, the Zhang-Levy-Fert model can be used to solve the Valet-Fert model in both the magnetic regions Ω m as well as the nonmagnetic regions Ω n . While both models perfectly agree in the bulk, the essential differences of the models are the continuity conditions of the potentials. While the Zhang-Levy-Fert model assumes continuous potentials u and s, the Valet-Fert model allows jumps which are defined by interface properties, namely the interface resistivity rb , the spin-flip probability γ sf and the spin-mixing conductance g Ö . However, these jumps can be mimicked in the Zhang-Levy-Fert model by introducing thin layers with effective material properties at the positions of the respective interfaces. Approximation of the charge current in the Zhang-Levy-Fert model (140) within the effective-interface layer Ω e by means of finite differences and multiplication with the boundary normal n yields with d being the thickness of the effective-interface layer.
Considering the potential mappings (176), (177) and the current mapping (174), this translates to the following jump condition across the effective-interface layer with d being the thickness of the layer. Comparison with the jump condition (165) of the Valet-Fert model results in the parameter mappings for the effective-interface layer. Applying the same procedure to the spin current of the Zhang-Levy-Fert model (139) yields which translates to and leads to the additional mappings According to these relations, the spin-mixing conductance g Ö depends implicitly on rb and γ sf which contradicts the Valet-Fert model where g Ö is an independent interface property. However, while the charge current j e is approximately constant throughout the effective-interface layer Ω e due to the continuity equation (143), the spin currentj s has sources in Ω e according to (142) which renders the finite-difference approximation (187) inaccurate.
While an appropriate choice of τ sf and J in the effectiveinterface layer can approximately reproduce the behavior of the Valet-Fert model, the dependency of these material parameters from the spin-mixing conductance g Ö is nontrivial and is best resolved by a fitting procedure. For the special case of a collinear magnetization configuration in the magnetic layers, the parameter mapping between the Zhang-Levy-Fert model and the Valet-Fert model is exact. In this case, the spin-accumulation in both models is also collinear to the magnetization. As a consequence, the spin-mixing-conductance term vanishes which reflects the fact, that a collinear magnetization configuration results in a vanishing spin torque. That said, applying the bulk mappings (174)-(181), and adding effective-interface layers with thickness d and material parameters (184)-(185) in order to simulate the interface properties of the Valet-Fert model in the Zhang-Levy-Fert model, results in the same spin accumulation and electric potential for both models in the limit of small d.
The Valet-Fert model is very popular in the experimental community where the interface properties rb , λ and g Ö are discussed and determined for various material systems. However, the Zhang-Levy-Fert model provides a more general approach for the description of spin transport and spin torque. In particular, the bidirectional description of spin torque that accounts for both the torque exerted from the current polarization on the magnetization and vice versa leads to a better representation of the physical processes. Interface properties as defined by the Valet-Fert model can be modeled with additional thin layers.

Beyond the spin-diffusion model
While the spin-diffusion model introduced in Section 5.3 has been shown to incorporate various models for the description of spin torque and other spintronics effects, its area of application is restricted to diffusive transport. Some effects that are not included in the equations presented in Section 5.3, such as inplane GMR, spin pumping and anomalous Hall effect might be added by means of additional terms to the diffusion model since they are in principle compatible with the diffusive transport assumption [70]. An important class of spintronics devices though is making use of magnetic tunnel junctions in order to exploit tunnel-magneto resistance (TMR) and spin torque. The spin transport in tunnel junctions, however, is not diffusive. Various ab initio models have been developed in order to accurately describe magnetic tunnel junctions [71][72][73]. However, ab initio models are computationally challenging and do not integrate well with the semiclassical micromagnetic model. The development of a suitable model that integrates well with the spin-diffusion model and micromagnetics is still subject to ongoing research.

Discretization
The micromagnetic model as introduced in the preceding sections defines a set of nonlinear partial differential equations in space and time, which can only be solved analytically for simple edge cases. In general, the solution of both static and dynamic micromagnetics calls for numerical methods. However, the development of efficient numerical methods is challenging due to the following properties of the micromagnetic equations.
1. The demagnetization field, that describes the dipoledipole interaction in the magnetic material, is a long-range interaction. A naive implementation of such an interaction has a computational complexity of Opn 2 q with n being the number of simulation cells. Various methods have been proposed to reduce this complexity to Opn log nq [74] or even Opnq [75]. 2. The exchange interaction adds a local coupling with high stiffness due to its second order in space. While the competition of the long-range demagnetization field with the local exchange field is crucial for the generation of magnetic domains, it also poses high demands on the numerical time-integration methods. 3. The nonlinear nature of most of the energy contributions leads to a complex energy landscape, which makes it difficult to efficiently seek for energy minima. In the context of quasistatic hysteresis computation, the nearest local energy minimum to a given magnetization configuration has to be found. A complex energy landscape increases the risk to miss a local minimum, and calls for a thoughtful choice of minimization algorithm.
A large variety of tailored numerical methods for the solution of micromagnetic equations has been proposed. Typically, these methods introduce distinct discretizations for space and time. Among the spatial discretizations, the most popular methods applied in micromagnetics are the FDM and the FEM. For both methods the magnetic region is subdivided into simulation cells resulting in a mesh of cells. However, the requirements for the mesh differ significantly for both methods. While the FDM usually requires a regular cuboid mesh, the FEM typically works on irregular tetrahedral meshes, see Figure 16. While FDM or FEM are used for spatial discretization in order to compute the effective field H eff or the respective energy contributions E, another class of algorithms is required in order to either minimize the total energy with respect to the magnetization configuration m or to compute the time evolution of m according to the LLG. Independent from the discretization method, the cell size has to be chosen sufficiently small in order to accurately resolve the structure of domain walls. The characteristic length for the domain-wall width is the so-called exchange length which is defined by where K eff is the effective anisotropy constant that includes contributions from the crystalline anisotropy as well as the shape anisotropy which is introduced by the demagnetization field [13]. Note that for both energy minimization and magnetization dynamics, the effective field needs to be computed in the magnetic domain only. In the following sections, the spatial discretization with FDM and FEM is discussed in detail. In further sections, numerical methods for efficient integration of the LLG, energy minimization and energy barrier calculations will be discussed.

The finite-difference method
The FDM is a very popular numerical tool for the solution of micromagnetic equations. While the restrictions to regular meshes renders this method inapt for certain problems involving complex geometries, this restriction allows the application of very fast algorithms.

Demagnetization field
Among the effective-field contributions introduced in Section 2, the demagnetization field holds the special role as the only long-range interaction. Since long-range interactions are computationally costly, the value of a spatial discretization strategy is significantly influenced by its demagnetization-field algorithm. The FDM solves partial differential equations by approximation of the differential operators with finite differences. In case of the demagnetization-field problem (16), which has the form of Poisson's equation, this would require the approximation of the Laplacian. However, the application of this classical finite-difference procedure is complicated by the open boundary condition (15), which prevents the restriction of the computational domain to the magnetic region Ω m . Hence, in finite-difference micromagnetics instead of discretizing Poisson's equation, the demagnetization field is usually solved by direct integration of (21). Consider a cellwise constant normalized magnetization with the cells Ω i being an arbitrary partitioning of the magnetic domain Ω m Inserting the discretization (192) into the integral formulation of the demagnetization field (21) yields In order to compute the demagnetization field H dem with the same discretization as the magnetization m, the field is averaged over each cell Ω i which results in where V i is the volume of the simulation cell i. Here, A denotes the linear demagnetization-field operator, which is represented by a dense 3nˆ3n matrix with n being the number of simulation cells. While this method could be used for the numerical demagnetization-field computation, it scales with Opn 2 q for both storage requirements and computational complexity which is unfeasible for large problems. A better scaling can be accomplished by exploiting the convolutional structure of the integral equation (21). In order to preserve this structure on the discrete level, a regular spatial discretization is required, i.e. all simulation cells Ω i must be of the same shape Ω ref  Figure 17. Consequently, the offset from a simulation cell Ω i to another simulation cell Ω j is given by a multiple of the cell spacing ∆x in every spatial dimension Using (198) and (199), the integration in (197) can be carried out over the reference cell Ω ref which has the form of a discrete convolution since it only depends on the difference of the multiindices i and j The discrete demagnetization tensorñ i´j has entries for every possible cell distance which amounts to ś k p2n k1 q « 8n, where n k denotes the number of cells in spatial dimension k. Compared to the direct demagnetizationfield operator A, the demagnetization tensor reduces the storage requirements from Opn 2 q to Opnq. The computational complexity, however, still amounts to Opn 2 q when implementing (201) literally.
In order to reduce the computational complexity, the discrete convolution (18) is computed in Fourier space where it reduces to a cell-wise multiplication according to the convolution theorem where the Fourier transform is applied componentwise.
Since the cell-wise multiplication has a low complexity of Opnq, the overall complexity of the demagnetizationfield computation is governed by the complexity of the Fourier-transform computation. In case of the fast Fourier transform this amounts to Opn log nq.
Note that this fast-convolution algorithm requires the discrete demagnetization tensorÑ i´j to be of the same size as the discrete magnetization m i in order to perform cell-wise multiplication in Fourier space. However, while the magnetization is discretized with ś i n i cells, the discrete demagnetization tensor has a size of ś i 2n i´1 . Hence, the discrete magnetization has to be expanded in order to match the size of the demagnetization tensor. Due to the cyclic nature of the discrete Fourier transform all entries of the demagnetization tensorÑ i´j are considered for every field evaluation H dem i . For instance, the negative-distance entryÑ´1 ,´1 is considered for the computation of the field H dem at position p0, 0q, although the p0, 0q cell has no neighbors at negative distances. In order to neglect these unphysical distances, the only reasonable choice for the expansion of the magnetization is by adding zero entries, which is often referred to as zeropadding, see [76]. The complete convolution algorithm for the demagnetization-field computation is visualized in Figure 18, where m andÑ are reduced to two spatial dimensions for the sake of simplicity. The result of the convolution algorithm is of the size ś i 2n i´1 like the demagnetization tensor. Physical meaningful values of the computed field H dem , however, are found only in the first ś i n i entries. The remaining entries are algorithmic byproducts and can be neglected.
Note that the only requirement for the application of the fast convolution is a mesh regularity as described by (198) and (199). However, the evaluation of the demagnetization tensor (202) might be unfeasible for complicated reference cells such as illustrated in Figure 17a. For threedimensional cuboid cells, an analytical formula for the demagnetization tensorÑ i´j was derived by Newell et al. [77]. According to this work the diagonal element N 1,1 of the tensor is given by where the function f is defined as The elements N 2,2 and N 3,3 are obtained by circular permutation of the coordinates N 2,2 px, ∆xq " N 1,1 rpx 2 , x 3 , x 1 q, p∆x 2 , ∆x 3 , ∆x 1 qs (207) N 3,3 pr, ∆rq " N 1,1 rpx 3 , x 1 , x 2 q, p∆x 3 , ∆x 1 , ∆x 2 qs. (208) The off-diagonal element N 1,2 is given by where the function g is defined as See equation (210) above.
Like the continuous tensorÑ px´x 1 q the discrete tensor N i´j is symmetric Hence, the above definitions of N 1,2 , N 1,3 and N 2,3 can be used to obtain the remaining off-diagonal elements. While these analytical expressions are exact, their numerical evaluation can lead to inaccuracies due to floating-point errors. These errors especially occur for large cell distances and degenerated cells and can be avoided by numerical integration as shown by Lebecki et al. [78] and Krüger et al. [79]. The fast-convolution algorithm can be further optimized by considering the specific properties of the demagnetization-field problem, i.e. Fourier transforms of zero values and evaluation of unneeded data can be omitted [80] and symmetries in the demagnetization tensor can be exploited [8] to further speed up computations.
Instead of computing the demagnetization field directly from a convolution of the magnetization m with the demagnetization tensorÑ , the field can also be computed as negative gradient of the scalar potential u [81,82]. The scalar potential u itself is the result of a convolution, see (20), and can be computed by the fast-convolution algorithm. The gradient can be approximated with finite differences. Compared to the direct computation, the advantage of this approach is the reduction of Fouriertransform operations. However, the additional gradient computation compromises the overall accuracy of the computation.

Local field contributions
Except for the demagnetization field, all energy contributions introduced in Section 2 have either local or short-range character in the sense that they depend on derivatives of the magnetization. Local contributions to the effective field, such as the anisotropy fields (68) and (70), are simply evaluated cellwise. Differential operators in short-range contributions such as the exchange field (62) are approximated with finite differences. The secondorder finite-difference approximations for the first and second spatial derivative of the discretized magnetization m i readˆB which results in the following discretization of the exchange field H ex (216) While higher-order finite-differences might lead to better convergence properties for some applications, they also increase the computational effort by involving more neighboring cells in the computation. Since second-order approximation is usually considered sufficiently accurate, most finite-difference codes stick with this approximation.
The boundary condition B " 0, as derived in Section 3, is usually implemented by introducing virtual cells surrounding the magnetic region Ω m , see Figure 19, and computing the magnetization values m i in these cells accordingly. In case of the exchange interaction this boundary condition is defined by Bm{Bn " 0 as derived in Section 3.2. For the boundary at x 1 " 0 this leads to the following equation

Bm
Bx i˙p 0,j,kq " m p1,j,kq´mp´1,j,kq 2∆x i " 0 (217) which needs to hold for all 0 ď j ď n 2´1 and 0 ď k ď n 3´1 and determines the values of the virtual cells m p´1,j,kq as m p´1,j,kq " m p1,j,kq . After computing the values of the virtual cells, the derivatives required for the distinct field computations can be evaluated according to (214) and (215) without special treatment for the boundary cells. Note that the boundary conditions, and thus the computation of the virtual cells, differ depending on the choice of field contributions, see Section 3.6.

Spintronics
The implementation of the spin-torque effects of Slonczewski as well as Zhang and Li as introduced in Sections 5.1 and 5.2 is straightforward as the torque contributions depend on the local magnetization and their derivatives only. In order to solve the spin-diffusion model, a second-order partial differential equation in space has to be solved in order to compute the spin-accumulation s and the electric potential u. The discretization of equations (142) and (143) with finite differences is straightforward. However, care has to be taken in order to properly account for discontinuous material parameters when dealing with multilayer structures. For instance, the spin-currentj s , as defined by (139), may be discontinuous across material interfaces due to discontinuities of the material constants C 0 and D 0 . Hence, the discretization of the divergence in (142) must be carefully chosen to account for these jumps in a distributional sense. A possible finite-difference discretization of the time dependent spin-diffusion model with prescribed electric current (138) and (141) is presented by García-Cervera et al. [83].

Existing software packages
Various software packages implementing the FDM with FFT accelerated demagnetization-field computation have been developed. Probably the most popular open-source finite-difference micromagnetic software is OOMMF [84]. OOMMF is a multiplatform code running on central processing units (CPUs). Other finite-difference CPU codes include the open-source software Fidimag [85] and the commercial package MicroMagus [86]. A very simple CPU implementation of the finite-difference algorithms with the Python library NumPy is presented in [87]. The recent advent of general-purpose graphicsprocessing units (GPGPUs) allowed for the significant acceleration of scientific software. A popular open-source package for finite-difference micromagnetics on GPGPUs is MuMax3 [88,89]. Other GPGPU codes include magnum.fd [90] and the recently developed GPGPU extension to OOMMF [82].

The finite-element method
The FEM is a powerful numerical tool for the solution of partial differential equations. In contrast to the FDM where the differential operators are discretized directly, in finite-elements the original problem is transformed into a variational problem before discretization. Consider Poisson's equation´∆ u " f in Ω. (218) Multiplying both sides with a test function v from a suitable function space V and integrating the original problem turns the original problem into a variational problem. The solution u P V is required to satisfý Applying integration by parts yields 220) where ∇u¨n in the boundary integral was replaced by u N in order to implement Neumann boundary conditions ∇u¨n " u N for x P BΩ. This formulation is referred to as weak form of the original problem as it weakens the and additionally restricting the test functions to V D given by so that the solution u is not tested on the Dirichlet boundary Γ D . By appropriate choice of the function space V , the variational problem (220), can be shown to have a unique solution [91]. Discretization of the continuous problem (220) is achieved by choice of a discrete function space V h Ă V . While the FEM is very general concerning the choice of discrete function spaces [91], the most common choice is that of piecewise affine, globally continuous functions. A suitable function basis is constructed using a tetrahedral mesh as depicted in Figure 16b. Each mesh node x i is associated with a basis function φ i with node values which is affine within each cell, see Figure 20. In order to discretize the weak formulation (220), both the solution function u and the test function v are expressed in terms of the basis functions Inserting into (220) and neglecting the boundary term for the sake of simplicity yields ż Ω f φ j dx @ j P r1, ns.
(226) Instead of testing with all possible test functions v h , the test functions are reduced to the individual basis functions φ j . Since both, the left-hand side and the right-hand side of (220) are linear in the test function v, the equality (226) holds also true for any test function v h . The discretized solution u i is given by a linear system of equations where the matrix A ij , which is referred to as stiffness matrix, depends only on the geometry of the mesh. Furthermore, A ij is sparsely populated due to the choice of basis functions that result in nonzero contributions only for neighboring nodes. The sparsity is an important aspect from a computational point of view, since it reduces the storage requirements from Opn 2 q to Opnq. Furthermore, this property can be exploited by the use of iterative methods for the solution of linear systems. These methods avoid the direct inversion of the matrix and require only the computation of matrix-vector multiplications which leads to a superior scaling compared to direct methods [92]. The procedure of turning a partial differential equation into a variational problem by multiplication with test functions, which is referred to as Galerkin method, can be applied to a large variety of problems.

Demagnetization field
As shown in Section 2.3, the demagnetization-field potential is given by Poisson's equation (16). However, in contrast to the example given in the preceding section, the demagnetization-field problem is subject to open boundary conditions (15). Hence, Poisson's equation has to be solved in the complete space R 3 which is not trivially possible with the FEM that applies only to finite regions.
A simple approach to approximate the required boundary conditions with finite elements is the so-called truncation approach. In order to compute the demagnetization field H dem of a finite magnetization region Ω m , the outer region R 3 zΩ m is approximated with a finite external region Ω e . Since the magnetic scalar potential is known to which solves the homogeneous Dirichlet problem in the domain Ω " Ω m Y Ω e . Note that the right-hand side of (232) is integrated by parts in order to loosen the continuity requirements on the magnetization m. This is essential in order to accurately account for the discontinuity of the magnetization at the boundary of the magnetization region BΩ m . If the magnetization m is discretized with the usual piecewise affine, globally continuous functions, this discontinuity can be described by restricting the integration domain of the right-hand side of (232) to the magnetic domain Ω m . The accuracy of the truncation solution is significantly influenced by the size of the external region Ω e . However, a larger size of the external region usually increases the number of mesh nodes and thus the computational effort.
Choosing Ω e five times the size of Ω m has been found to be a reasonable trade-off between accuracy and computational complexity [93]. In order to further reduce the computational costs, the mesh in the external region can be gradually coarsened to the outside, see Figure 21. This leads to a reduction of degrees of freedom and has no significant impact on the accuracy of the solution, since the scalar potential is known to smoothly decay to the outside.
A related approach to the truncation method is the socalled shell-transformation method. Instead of describing the exterior space R 3 zΩ m with a relatively large region Ω e , this method uses a rather small shell region Ω s and employs a transformation T that maps the shell region onto the complete exterior space Specifically, the transformation T maps the inner boundary of the shell region onto itself and the outer outer shell boundary onto infinity For a spherical shell, the transformation is illustrated in Figure 22. By substitution of variables, integrals over the exterior space R 3 zΩ can be expressed by integrals over the shell region Ω s . For the left-hand side of the weak formulation (232) this substitution yields where g is the so-called metric tensor defined by with J " ∇T being the Jacobian of the transformation T . This leads to the weak formulation for the solution of the scalar potential u. Various shell geometries and transformations have been proposed to solve (238) [94][95][96][97]. Like for the truncation method, application of the shell-transformation method results in sparse linear systems. The advantage of the transformation method over the truncation approach is a finite representation of the infinite exterior region. However, the metric tensor g is by definition singular at the outer shell boundary which makes the computation of the discrete entries for the system matrix difficult. Furthermore, this singularity leads to a bad matrix condition, which has a negative effect on the numerical solution of the linear system.
In order to further reduce the degrees of freedom, it is desirable to consider only the magnetic region Ω m for spatial discretization. This can be achieved by a hybrid finiteelement-boundary-element method proposed by Fredkin and Koehler [98]. From (14) and the divergence theorem the following jump conditions for the scalar potential at the boundary of the magnetic region BΩ m apply where u´denotes the value in the magnetic region and u`denotes the corresponding value outside. Consider the following splitting of the potential u where u 1 is solved by ∆u 1 "´∇¨pM s mq (242) in the magnetic region Ω m and zero in the exterior region R 3 zΩ m . While u 1 satisfies the jump condition (240), it violates the continuity condition (239). Thus, u 2 has to be chosen to fix (239) while preserving (240) and being a solution to the Laplace equation ∆u 2 " 0, i.e.
u2´u2 " u1 These requirements are fulfilled by the double-layer potential defined as While this integral expression can be used to evaluate u 2 in Ω m , this procedure has a high computational complexity of Opn 2 q. Hence, (247) is only evaluated at the boundary BΩ m using the boundary-element method. The result of this computation is then used as Dirichlet boundary condition for a Laplace problem that is solved with finite-elements. This means that a computation of the scalar potential u requires the finite-element solution of Poisson's equation with Neumann boundary conditions, followed by a boundary-element evaluation of the doublelayer potential and the finite-element solution of a Laplace problem with Dirichlet conditions. Compared to the pure finite-element presented in this section, the advantage of this hybrid method is that only the magnetic region Ω m has to be discretized. This is even true if Ω m consists of various separated domains. The disadvantage, however, is the loss of sparsity. The boundary-element operator for the double-layer computation is a dense n BˆnB matrix with n B being the number of boundary nodes. A common procedure to reduce the storage requirements of this matrix is the application of hierarchical matrices [99,100] which reduces the storage requirements as well as the computational complexity for the potential evaluation from Opn 2 B q to Opn B log n B q. An alternative hybrid method using the single-layer potential was proposed by García-Cervera and Roma [101].
The methods presented in this section are devoted to the calculation of the magnetic scalar potential u rather than the demagnetization field H dem which is defined as the negative gradient of the potential H dem "´∇u.
Since the discrete solution of u is a piecewise affine function, the gradient of u can be easily computed. However, the gradient of a piecewise affine function is a piecewise constant, discontinuous function, see Figure 23. Often it is desirable to compute the demagnetization field with the same discretization as the magnetization. This is required for nodewise operations such as the nodewise evaluation of the right-hand side of the LLG equation (117). The approximation of the demagnetization field as a piecewise affine, globally continuous function can be achieved by solving the weak form with H dem P V . This procedure is referred to as projection. Discretization of (248) with the usual basis functions The system matrix resulting from this weak form is called mass matrix M Like the stiffness matrix arising from discretization of the Laplacian, the mass matrix is sparse since only basis functions of neighboring mesh nodes give nonzero contributions. In order to avoid solving another linear system for the computation of the demagnetization field, the mass matrix can be approximated by a diagonal matrix Mg iven by with 1 " p1, 1, 1q. This approximation is referred to as mass lumping since all off-diagonal mass entries are lumped in the diagonal This approximation preserves the integral of the solution variable H dem which justifies its application Since the lumped mass matrix is diagonal, its inverse is trivially given by the reciprocal diagonal entries. Hence, instead of solving a linear system, the projection can be expressed as a single sparse matrix-vector multiplication The illustrative description of mass lumping is that of an weighted average with the weight being the cell sizes as described in [1]. The result of a gradient projection is depicted in Figure 23 along with the original function.

Local field contributions
Local field contributions are usually computed with a similar procedure as applied for the gradient computation in Section 6.2.1 For the exchange field defined by (62) this procedure yields where integration by parts was applied in order to avoid second spatial derivatives. Due to the boundary condition (63), the boundary integral on the right-hand side vanishes which results in While this weak form poses no further difficulties for single-phase magnets with a constant saturation magnetization M s , the gradient on the test function diverges at material interfaces with discontinuous M s . In order to avoid this problem, the original problem (256) can be multiplied with´µ 0 M s before integrating by parts resulting iń with δE ex being the exchange differential defined in (60). This method can be used to compute any local field contribution introduced in Section 2. In order to avoid the solution of a linear system for the retrieval of field contributions, the mass-lumping procedure introduced in Section 6.2.1 can be applied. For energy contributions quadratic in the magnetization m, such as the exchange field, the field computation can be reduced to a single sparse matrix-vector multiplication given by Note that this method accounts for the boundary conditions introduced by the exchange and antisymmetric exchange field in a generic fashion since it considers the differential of the energy δE rather than the functional derivative δE{δm. Since the boundary conditions are embedded in the differential, they do not have to be explicitly applied. While some effective-field contributions like the uniaxial-anisotropy field (68) can also be computed nodewise for single-phase magnets, the nodewise computation fails for discontinuous material parameters. Due to the integral formulation, the calculation according to (260) and (261) also works in these cases.

Interface contributions
Some energy contributions introduced in Section 2 are interface effects rather than bulk effects. This means that the energy connected to these effects is the result of an interface integral. These interface contributions usually enter the micromagnetic formalism through boundary conditions, see e.g. Section 3.4. However, in the framework of dynamical micromagnetics, it is often desirable to be able to express all energy contribution in terms of an effective field instead of a boundary condition in order to add up the different contributions for their use in the LLG.
In order to compute a discretized effective field that accounts for interface energy contributions, we require the energy generated by the effective field to equal the interface energy contribution. Consider the interface exchange energy given by (38). Since this energy contribution is quadratic in m, the equality of energies reads (262) with p " AmrP pxqs. Discretization of the field H eff yields which has to hold for arbitrary magnetizations m. Hence, the magnetization can be replaced with basis functions and the equality has to hold for any φ j . Application of mass lumping yields the system which has the exact same form as any local bulk-field contribution as shown in Section 6.2.2.

Spintronics
As for the FDM, the simplified spin-torque models by Slonczewski as well as Zhang and Li can be incorporated in finite-element micromagnetics along the lines of the local field contributions shown in Section 6.2.2. In order to account for the Slonczewski spin-torque as an interface effect, the consideration of Section 6.2.3 can be applied.
A more challenging task is the discretization of the spin-diffusion model introduced in Section 5.3. The spindiffusion model requires the computation of the spin accumulation s which is the solution to a partial differential equation. This equation is turned into a weak form using the Galerkin method. In order to solve the equilibrium spin accumulation s for a given magnetization m and a prescribed charge current j e as defined by (141) and (142), the following weak form applies Note that integration by parts was not only applied in order to eliminate second derivatives of the spin accumulation s, but also in order to eliminate any derivative on the magnetization m or material parameters since these variables may have discontinuities in the problem domain Ω. While material parameters may have arbitrary discontinuities at material interfaces, the magnetization is continuous within the magnetic domain Ω m . However, since the magnetization vanishes in the nonmagnetic domain ΩzΩ m , it is by definition discontinuous across magnetic-nonmagnetic interfaces. For the discrete formulation, these properties are take into account by appropriate choices of function spaces and integration domains. Namely, all material parameters are discretized with piecewise constant functions. Since the magnetization is continuous within the magnetic region, it is discretized with the usual piecewise affine, globally continuous functions. In order to account for the discontinuity at magnetic-nonmagnetic interfaces, any integral including the magnetization m is restricted to the magnetic domain Ω m which is equivalent to a rapid drop of the magnetization to zero outside the magnetic domain. The boundary integrals arising from partial integration of ∆s vanish due to the homogeneous Neumann boundary condition on s, see Section 5.4. Since the spin accumulation s is assumed to be continuous even across material interfaces, both s as well as the test functions are discretized with componentwise piecewise affine, globally continuous functions s, v P V . For the solution of the self-consistent spin-diffusion model given by the current definitions (140), (139) and their respective source equations (143), (143), the function space for both the solution and the test functions has to be extended. The solution of the self-consistent model comprises both the spin accumulation s and the electric potential u. Both the components of s and the scalar field u are discretized with piecewise affine, globally continuous functions ts, uu P VˆV . Applying the Galerkin method and performing integration by parts in order to avoid diverging derivatives yield two coupled weak formulations. The weak formulation of (140) and (143) reads and the weak formulation of (139) and (142) is given by Here, Γ N P BΩ denotes all external interfaces that act as contacts with prescribed charge-current inflow j 0 e¨n and Γ D P BΩ denotes contacts with prescribed electric potential u 0 . While the current inflow is implemented as natural Neumann boundary condition, the prescribed potential is implemented as Dirichlet boundary condition, see Section 5.4. In addition to the boundary conditions on the potential u, homogeneous Neumann conditions are applied to the spin accumulation s. Discretization of (269) and (270) yields a single sparse system of the size 4nˆ4n with n being the number of mesh nodes. Incorporating spin-orbit interactions given by (151) and (152) or spin dephasing given by (155) is straightforward and can be done along the lines of the weak forms (268)-(270).

Existing software packages
The implementation of finite-element solvers is a challenging task, since it involves the nontrivial generation of tetrahedral meshes, the numerical computation of integrals for the system-matrix assembly and the solution of large linear systems. Various software packages and libraries have been developed in order to solve one or more of these tasks. Popular open-source packages for the mesh generation are Gmsh [102] and NetGen [103]. With ONELAB [104] and NGSolve [105], these mesh generators also act as full stack finite-element libraries. Other open-source libraries for the formulation and solution of finite-element problems are Escript [106] and MFEM [107]. A very comprehensive and fast, yet easy to use finite-element library is FEniCS [108] which, by default, uses PETSc [109] as linear-algebra backend. Libraries for the boundary-element method as required by the hybrid demagnetization-field method introduced in Section 6.2.1 include BEM++ [110] and H2Lib [111]. Both libraries are open source and provide routines for the assembly of boundary-element matrices and their compression with hierarchical matrices.
Besides these multipurpose libraries, a number of specialized micromagnetic finite-element packages are available. The open-source library FinMag [112] and the closed-source library magnum.fe [113] are micromagnetic simulators based on FEniCS. While FinMag concentrates on classical micromagnetics, magnum.fe also implements the spin-diffusion model [114][115][116]. Another closed-source finite-element code that solves the micromagnetic equations coupled to the spin-diffusion model is FEELLGOOD [117,118]. Other finite-element codes include the opensource packages Magpar [119] and NMag [120] as well as the closed-source package FEMME [121]. The closedsource packages Tetramag [122] and Fastmag [123] make use of graphics processing units (GPUs) to speed up computations.

Other spatial discretization methods
While the FDM and the finite-element are by far the most common discretization methods used in the micromagnetic community, other methods have been proposed to solve parts of the micromagnetic model. Especially the computation of the demagnetization field is an ongoing matter of research.
A method that aims to combine the speed of the FFT based convolution with the flexibility of irregular meshes is the non-uniform FFT [124,125]. Both the FFT accelerated convolution and the finite-element demagnetization-field computation exhibit at least a computational complexity of Opn log nq. A well known method for the interaction of particle clouds which scales with n is the fast multipole method which has been shown to be also applicable to the continuous demagnetization-field problem [126,127]. Another class of methods that scales below Opn log nq employs low-rank tensor approximations [128,129].

Time integration
Numerical integration of the LLG equation (117) poses several challenges on the applied method. One of these challenges is the high stiffness that is introduced by the exchange interaction [130]. Another challenge is the micromagnetic unit-sphere constraint |m| " 1 that is required to be preserved by a time-integration scheme. Several methods, that address these difficulties, have been proposed in order to efficiently integrate the LLG.
Most numerical time-integration methods used in micromagnetics completely separate the spatial discretization from the time discretization. For these methods, both the magnetization m and its time derivative B t m are represented by vectors and the differential equation in space and time is transformed into n coupled ordinary differential equation In order to evaluate the time derivative f i , the effectivefield contributions are computed according to the methods introduced in the preceding sections, and the righthand side of the LLG (117) or (118) is evaluated cellwise/nodewise. Due to its first order in time, the LLG is an initial value problem. We denote the initial value of the magnetization at time t 0 as A discrete time-integration scheme approximates the magnetization dynamics as a series of magnetization snapshots at times t i that we denote by In the following, the most common time integrators in micromagnetics will be discussed in detail. A more comprehensive overview of methods can be found in [41].

Explicit Runge-Kutta methods
A reliable and efficient class of numerical integrators for initial value problems are the explicit Runge-Kutta methods. According to the most general form of the Runge-Kutta method, the magnetization configuration m i is obtained from a known magnetization configuration m i´1 by where a set of coefficients b j , c j and a ij defines a particular Runge-Kutta method. The auxiliary results k j are computed one after another. Due to the restriction j ă k in the summation on the right-hand side, each k j depends only on previously computed k k and the previous magnetization m i´1 which makes this method explicit. Explicit methods are usually computationally cheap, since they are linear in the solution variable m i . However, explicit methods lack stability which might be a disadvantage especially in the case of stiff problems. The simplest Runge-Kutta method is the explicit Euler method which is obtained by setting b 1 " 1 and c 1 " 0 This first-order method is not a good choice in terms of accuracy and efficiency. However, it is well suited to investigate the preservation of the unit-sphere constraint.
Inserting the LLG (117) into (276) results in Multiplying with m i`mi´1 and reinserting (276) yields and thus |m i | ď |m i´1 | where the preservation of norm |m 1 | " |m i´1 | only holds for a vanishing right-hand side of the LLG B t m " 0. In order to enforce the norm preservation, the magnetization is usually renormalized after each Runge-Kutta step, i.e. the magnetization m i is replaced by m i 1 given by This renormalization is also performed for higher order methods. While these methods reduce the violation of the unit-sphere constraint, they do not guarantee its preservation. Higher-order schemes are obtained by the choice of appropriate parameters b i , c j and a ij . The classical Runge-Kutta method requires the computation of four auxiliary results k j and is of fourth order. Since the evaluation of each k j comes at the price of an effective-field evaluation, the computation of a single integration step is more expensive for higher-order methods than for lowerorder methods. However, usually this disadvantage is more than compensated by the size of the time step which may by significantly larger for higher-order methods without compromising accuracy. While the classical fourth-order Runge-Kutta method is very efficient, the choice of time step is difficult and may even change in the course of integration. This problem can be solved by application of Runge-Kutta methods with adaptive stepsize control. These methods derive differentorder approximations from a shared pool of auxiliary results k j and estimate the integration error by comparison of these approximations. Prominent candidates which have proven valuable for micromagnetics are the Runge-Kutta-Fehlberg method [131] and the Dormand-Prince method [132] that both use fourth/fifth-order approximations for the stepsize control.
The usage of explicit methods is usually not advised for stiff problems because of their lack of numerical stability [133]. In micromagnetics, a strong exchange coupling can lead to a very high stiffness of the differential equation. However, since the time-step size is coupled to the size of the smallest mesh cell, the use of a regular grid can mitigate this problem [130]. Hence, Runge-Kutta methods are the standard choice in finite-difference micromagnetics [8].

Implicit midpoint scheme
An implicit integration scheme that specifically accounts for the micromagnetic unit-sphere constraint is the implicit midpoint rule [40]. According to this method, the magnetization at m i is obtained from the previous magnetization snapshot m i´1 by Inserting the time derivative of m according to the LLG (117) and setting B t m " pm i´mi´1 q{∆t on the righthand side yields Multiplying both sides with m i`mi´1 immediately yields |m i | " |m i´1 |. Hence, the midpoint rule exactly preserves the magnetization norm for arbitrary time-step sizes. Moreover, it can be shown by a similar procedure, that the midpoint scheme preserves the energy of a magnetic system with quadratic energy contributions in m for vanishing damping α [40]. However, the implicit nature of this method comes at the price of nonlinearity in the solution variable m i . Since the nonlinear system (281) cannot be solved analytically, it requires the application of an iterative procedure such as Newton's method for the computation of m i . Each Newton iteration requires the evaluation of the effective field, which results in a high computational effort.

Tangent-plane integration
Another class of integrators especially suited for the application in the framework of finite elements was introduced by Alouges [134]. This method relies on an alternative formulation of the LLG which is obtained by crossmultiplying the LLG (117) with m αB t m`mˆB t m " γH eff´γ pm¨H eff qm.
This formulation is equivalent to the original LLG since all terms in (117) are perpendicular to m. In order to solve for the time derivative B t m, (282) is reformulated in a weak form. However, instead of seeking for the solution to w " B t m in the complete solution space V : R 3 Ñ R 3 , the solution space is restricted to the tangent space of the magnetization V T " tv : v¨m " 0u. This also allows for the restriction of the test space to the same space V T which simplifies the weak formulation of (282) to Instead of the original LLG, the right-hand side of this form is of the same order in the magnetization m as the effective field H eff . This feature can be exploited in order to construct an implicit integration scheme for field terms linear in m without losing the linearity of the weak form in w. The time derivative w at t i´1 is obtained by setting m " m i´1`θ ∆tw (284) with 0 ď θ ď 1 where θ " 0 leads to an explicit scheme and θ " 1 leads to an implicit scheme. Employing the considerations concerning discontinuous material parameters in Section 6.2.2 yields the weak form ż Ωm µ 0 M s pαw`m i´1ˆw q¨v dx"γδEpm i´1`θ ∆tw, vq.
(285) Considering only the exchange field, the weak formulation reads Each field contribution can be treated with an individual θ. While the exchange field adds a high measure of stiffness to the problem which calls for an implicit integration scheme, other field terms may well be treated explicitly [135]. This usually applies to the demagnetization field. Despite its linearity in m the demagnetization field cannot be treated implicitly in the same manner as the exchange field, since it is the solution to another partial differential equation, see Section 6.2.1. After computing the time derivative w " B t m, the actual time step is performed by computing where the right-hand side is evaluated nodewise. Various techniques have been proposed in order to implement the tangent-plane function space V T within a finite-element formulation. A possible solution is the application of Lagrange multipliers [113,135] that enforce the orthogonality condition on w. An alternative approach uses a local mapping of two-dimensional vector fields R 3 Ñ R 2 onto the tangent plane V T [136]. A combination of this integration scheme with a time marching scheme for the spin accumulation yields a method for the coupled solution of micromagnetics with the dynamic spin-diffusion equation (138) [114,137].
Variants of this method include higher-order methods [117,124]. However, by increasing the order of the method, the linearity in the solution variable w is lost, which leads to a higher computational effort.

Backward differentiation formula
All previously introduced methods are so-called singlestep methods, that require the knowledge of a single magnetization snapshot m i´1 in order to compute the subsequent magnetization m i . In contrast, the backward differentiation formula (BDF) is a multi-step method, i.e. the magnetization m i is computed from a series of preceding snapshots m i´j with 0 ď j ď s and s being the order of the method. In the most general form, the BDF method is given by where a i and β can be chosen such that the method is of order s. Normalizing the parameter a i and β so that a 0 "´1 yields which is a nonlinear equation in m i that can be solved with Newton's method. Starting from an initial configuration m i,0 that is usually approximated by extrapolation of previous snapshots m i´j , the Newton iteration reads where the variational derivative δG{δm i is given by In order to avoid the numerically expensive inversion of δG{δm, the discretized version of the linear equation (290) is usually solved iteratively for ∆m i,k " m i,k`1´mi,k . However, in order to reduce the number of required iterations, it is advised to employ a preconditioning procedure to (290). Instead of solving (290) directly, the preconditioned system is solved for P ∆m i,k . The preconditioner P is chosen to approximate the original problem while being easily invertible. A good choice for this purpose is a simplification of (291) where the linearization of the LLG δw{δm only includes local and linear effective-field contributions like the exchange field. Following this procedure, the computation of a single time step is computationally expensive since for every Newton iteration, the preconditioning system has to be solved. However, this method has proven to be an superior choice for a large range of problems in the framework of finite-element micromagnetics [130].

Existing software packages
The implementation of explicit integration schemes, such as the Runge-Kutta scheme introduced in Section 6.4.1, is straightforward and does not necessarily require the use of external libraries. However, implicit methods call for the solution of nonlinear systems. Linear algebra libraries such as PETSc [109] provide useful functionality for the implementation of suitable Newton methods. The opensource library SUNDIALS [138] includes a very efficient implementation of an adaptive-time-step BDF integration scheme as introduced in Section 6.4.4.

Energy minimization and barrier computation
While dynamical micromagnetic simulations are useful to gain insight into the mechanism of fast magnetization processes like magnetization switching, they are not feasible for investigations in the MHz regime and below. A typical example for quasi static micromagnetics is the computation of hysteresis curves. Experimental measurements of hysteresis curves are performed with field sweeps that are orders of magnitude slowers than the response time of magnetic systems which is in the GHz regime. Hence, the hysteresis curve can be computed by energy minimization instead of dynamical simulation. This is done by increasing the external field stepwise end computing the new energy minimum for each step. In order to compute hysteresis properties, the minimization algorithm is required to converge into the nearest local energetic minimum starting from a given magnetization configuration.
A global minimization would yield unique magnetization configurations for a given external field and hence be unsuited for hysteresis computations. Numerical minimization is usually performed iteratively based on gradient evaluations. These gradient methods are particularly useful for hysteresis computation since they start from a given magnetization configuration and converge to local minima. The simplest gradient based minimization algorithm is the steepest-descent method with a single iteration given by with τ being the step size. In order to account for the micromagnetic unit-sphere constraint |m| " 1, the descent direction is usually projected onto the tangent space of the magnetization Inserting into (293) and considering the definition of the effective field (116) yields which is equivalent to an integration step of the LLGdamping-term with an explicit Euler method. While this method still violates the unit-sphere constraint as shown in Section 6.4.1, it represents a significant improvement compared to the original steepest-descent method (293). Instead of using an explicit Euler method, any of the integration methods introduced in Section 6.4 may be used to progress the steepest descent. However, time-integration methods are optimized to accurately solve for the complete magnetization trajectory, while the only measure of interest for a minimization problem is the final magnetization configuration. Various tailored approaches for micromagnetic energy minimization including optimized steepest descend methods and variants of the conjugate gradient method have been proposed in order to achieve high performance and reduce the risk to miss local minima [28,[139][140][141].
Another class of minimization methods addresses the computation of energy barriers between two given local minima. Energy barriers are an important measure for the stability of magnetic devices. In magnetic storage devices, binary information is usually stored by putting the system in one of two possible energy minima, e.g. the up or down configuration in an anisotropic magnetic layer. The energy barrier between those minima defines the characteristic lifetime of the information by the Néel-Arrhenius equation with E being the energy barrier, T being the temperature and τ 0 being the attempt time.
A robust yet simple method for the numerical calculation of the energy barrier between two magnetic states is the string method that yields the minimum energy path between two states [142]. The string method is an iterative method that evolves an initial magnetization path mpϕq toward a minimum energy path defined bý where ϕ denotes the position on the path and the K subscript denotes magnetization variations that are orthogonal to the path mpϕq. The magnetization path is discretized by a finite number of magnetization images mpϕq Ñ m i with i P t0, 1, . . . , nu.
As illustrated in Figure 24, these images are evolved individually. For a single string-method step, each magnetization image is relaxed toward an energetic minimum by a certain amount. Using the projected steepest descent method (295), the updated images are given as with ζ being the step size of the descent method. Evolving each magnetization image toward the nearest local minimum would eventually relax every image into either the first minimum m 0 or the second minimum m n . In this case all information about the transition path, and thus also the barrier information, would be lost. In order to prevent the images to separate on the transition path, each evolution step (299) is followed by a reparametrization of the path. The position ϕ i of each magnetization image m i on the transition path is determined by the pairwise distance norm After determining the position φ i , the magnetization images m i are interpolated onto the regular grid by cubic spline interpolation. Various variants of this method have been proposed in order to optimize the accuracy and convergence properties. An effective improvement is the use of an energy weighted norm for the reparametrization in order to increase the image density in the barrier regime. A popular alternative to the string method is the nudged elastic band method that introduces a spring force between the images in order to preserve a homogeneous discretization of the transition path [143].

Applications
This section is dedicated to applications of micromagnetic simulations. While the first application deals with a classical micromagnetic problem and discusses the performance differences of finite-difference and finite-element tools, the remaining examples focus on the simulation of spintronics effects.

Standard problem #4
A well known dynamical micromagnetic problem is the standard problem #4 [144] that was developed by the µMAG group with the aim to serve as a benchmark problem for micromagnetic simulation tools. The problem considers a cuboid shaped magnetic system with dimensions 500 nmˆ125 nmˆ3 nm and material parameters similar to permalloy M s " 8ˆ10 5 A{m, A " 1.3ˆ10´1 1 J{m, Fig. 25. Time evolution of the average magnetization components for the switching process of a permalloy thin film according to the first part of the µMAG standard problem #4 computed with a finite-difference solver (FDM) and a finite-element solver (FEM).
K u1 " 0 and α " 0.02. The system is prepared in a socalled s-state e.g. by relaxing the initial magnetization m 0 "`cosr0.1s, sinr0.1s, 0˘into an energetic minimum. After relaxation, an external field with a magnitude of 25 mT is applied directed 170˝counterclockwise from the positive x-axis. This field results in the switching of the magnetization in the thin film. The dynamics of the switching process should be resolved by numerical integration of the LLG. An alternative field which is defined in the original problem specification in not considered in this work. Figure 25 shows the time evolution of the averaged magnetization components computed with the finitedifference solver magnum.fd [90] and the finite-element solver magnum.fe [113] that exhibit good agreement. The finite-difference solver employs the FFT-accelerated demagnetization-field method presented in Section 6.1.1 and an adaptive Runge-Kutta-Fehlberg method of 4/5 order for the time integration, see Section 6.4.1. The FEM employs the hybrid FEM/BEM method presented in Section 6.2.1 and an adaptive preconditioned BDF scheme for time integration, see Section 6.4.4. For the FDM we choose a discretization of 200ˆ50ˆ1 " 10 000 cells and for the FEM we use a regular tetrahedral grid based on a 100ˆ24ˆ2 cuboid grid which leads to 7878 mesh nodes. Both grids are chosen considering the exchange length of permalloy in order obtain accurate results.
Due to the different algorithms for demagnetizationfield computation and time integration, the finitedifference solver and the finite-element solver show significantly deviating computation times for the problem, see Figure 26. While the finite-element solver spends more than half of the computation time on the assembly and solution of the preconditioner of the time integration scheme, the finite-difference solver is completely dominated by the demagnetization-field computation.
Another interesting comparison of the two methods is illustrated in Figure 27 where the required number of LLG right-hand-side evaluations and the total simulation time is compared for different problem sizes. Starting from the above discretization, we perform mesh refinement  on both the finite-difference grid and the finite-element mesh, which leads to an increased stiffness of the problem. For the implicit time integration scheme used in the finite-element code, the number of right-hand-side evaluations remains almost the same for a given integration accuracy. In contrast, the number of right-hand-side evaluations of the finite-difference solver increase faster than linear with the number of simulation cells. However, despite the much larger number of right-hand-side evaluations of the FDM, it still beats the FEM with respect to the overall simulation time. One of the reasons for the superiority of the FDM for this example is the thin shape of the problem domain. While this reduces the demagnetization-field computation to a two-dimensional convolution in the case of the finite-difference solver, the demagnetization-field algorithm of the finite-element code is dominated by the boundary-element method that has inferior scaling properties.

Standard problem #5
The first and only µMAG standard problem dealing with spintronics effects is the standard problem #5 [145] which is derived from a work of Najafi et al. [146]. It describes the precessional motion of a magnetic vortex core due to spin-torque. A cuboid magnetic thin film, with dimensions 100 nmˆ100 nmˆ10 nm and the material parameters of permalloy, see Section 7.1, is initialized in a magnetic vortex state. The application of a homogeneous charge current leads to the precessional motion of the vortex core which eventually relaxes in a shifted equilibrium position. The problem definition suggests the application of the model of Zhang and Li as introduced in Section 5.2 for the calculation of the magnetization dynamics. The current j e is driven through the sample in x-direction. Setting the product of the coupling constant and the current density bj e " 72.17 m{s and choosing the degree of nonadiabaticity ξ " 0.05 and the damping α " 0.1 yields the damped precessional motion depicted in Figure 28.
As shown in Section 5.9.2, the spin-diffusion model is equivalent to the model of Zhang and Li in the case of a vanishing diffusion constant D 0 . The equality of the two models can also be achieved with a finite D 0 by considering the limit of vanshing λ sf and λ J while preserving the ratio λ 2 J λ 2 sf " ξ " 0.05.
In this limit, the terms linear in D 0 become neglectable and the resulting torque is the same as for vanishing D 0 . In order to investigate the influence of diffusion effects on the magnetization dynamics we consider a finite diffusion constant of D 0 " 10´3 A{m which is a reasonable choice for magnetic materials [147]. Furthermore, we choose β 1 " 0.8 and β " 0.9 which yields a current density of j e " 1.15ˆ10 12 A{m 2 according to the required definition of bj e . We perform dynamic simulations for different λ sf while always choosing λ J to fulfill (302). The resulting equilibrium magnetizations are summarized in Figure 29a. While the simulations with small diffusion lengths λ sf and λ J show a perfect agreement with the results obtained from the model of Zhang and Li, there are significant deviations for larger diffusion lengths. These deviations are most significant in the x-shift of the vortex core which is characterized by the y-component of the averaged magnetization. The same comparison for a different degree of nonadiabacity ξ " 0.5 yields significant deviations in the y-shift of the vortex core, see Figure 29b. These simulations demonstrate that the model of Zhang and Li is not accurate in the presence of spin diffusion. While the simplifications of the Zhang-Li model make it very attractive for computational micromagnetics, the results obtained from this model have to be evaluated with care. Simulation results for materials with large diffusion lengths might be useful for qualitative investigations. However, for an accurate description, the solution of the full diffusion model should be considered.

Spin-torque oscillator
Another spin-torque driven device with various potential applications is the spin-torque oscillator (STO) [48,49]. A simple STO consists of two magnetic layers separated by a nonmagnetic layer. One of the magnetic layers is designed to be very hard magnetic in order to act as stable spin polarizing layer. The other magnetic layer, referred to as free layer, is soft magnetic and usually stabilized by an external field. By applying a current perpendicular to the layer system, the free layer is subject to spin torque which drives its magnetization out of the externalfield direction. By suitable choice of current strength, the free-layer magnetization performs a precessional motion due to the external field while the damping contribution of the external field is exactly compensated by the spin torque, see Figure 30. Numerous variants of STOs have been proposed and simulated with the spin-torque model of Slonczewski [54,[148][149][150]. The model of Slonczewski is in principle perfectly suited for the simulation of STOs. However, the input parameters to this model are the angular dependencies η damp and η field that cannot be trivially Fig. 30. Schematic illustration of a field-stabilized spin-torque oscillator. The pinned layer acts as spin polarizer and the spin-torque compensates the damping contribution from the external field. derived from the geometry and material parameters of the system.
The spin-diffusion model renders very useful for the investigation of geometry dependent properties of such devices. We simulate a cylindrical system with a radius of R " 30 nm, a free-layer thickness of d free " 3 nm and a spacer-layer thickness of d spacer " 1.5 nm. The material parameters of the magnetic layers are chosen as α " 0.1, A " 2.8ˆ10´1 1 J{m, D 0 " 10´3 A{m, β " 0.8, β 1 " 0.9, λ sf " 10 nm, λ J " 2.24 nm. For the free layer we further choose µ 0 M s " 1 T and K u1 " 0 and for the spin-polarizing layer we chose µ 0 M s " 1.24 T and K u1 " 10 6 J{m 3 with a perpendicular anisotropy axis. The nonmagnetic spacer layer is simulated with material parameters D 0 " 5ˆ10´3 A{m and λ sf " 10 nm. We set the external field µ 0 H " 0.6 T in z-direction and a current density j e " 4ˆ10 11 A{m 2 in negative z-direction as shown in Figure 30. This system is simulated with varying polarization-layer thicknesses with the spin-diffusion model with prescribed current density. The resulting oscillation frequencies and tilting angles are shown in Figure 31. The qualitative dependence of the frequency and tilting angle from the thickness of the polarizing layer does not come as a surprise, since a thicker polarizing layer will obviously lead to higher spin polarization of the electrons in the free layer. However, the spin-diffusion model does not only account for geometry changes, but also for the changes of material parameters in distinct layers. In this respect it outperforms the simple model of Slonczewski that, on the other hand, has a much lower computational complexity which makes it a fast alternative to the spin-diffusion model for many problem settings.

Spin-orbit torque MRAM
In a final numerical experiment, we demonstrate the capabilities of the self-consistent spin-diffusion model with spin-orbit interactions. We aim to simulate the write and read process of a perpendicular spin-orbit torque magnetoresistive random-access memory (SOT MRAM) [151,152]. Consider a circular magnetic multilayer with a radius of R " 10 nm consisting of 4 layers with thicknesses 1 nm, 0.5 nm, 1 nm and 4 nm from bottom to top, see Figure 32. From these 4 layers, only the bottom layer (free layer) and the third layer from below (pinned layer)  are magnetic while all layers are conducting. The circular stack is centered on top of a rectangular conducting underlayer with dimensions 50 nmˆ50 nmˆ10 nm and the complete structure is meshed with prescribed mesh size of 2 nm. We define the two faces of the underlayer lying in the yz-plane as contact 1 and contact 2 respectively and the top interface of the circular stack as contact 3. Material parameters in the conducting underlayer are chosen as D 0 " 10´3 m 2 {s, C 0 " 6ˆ10 6 A{Vm, τ sf " 2ˆ10´1 5 s, θ " 0.3 which is typical for heavy metals such as Ta that give rise to the spin Hall effect. For the magnetic free layer we choose material parameters typical for perpendicular MRAM, namely M s " 0.796ˆ10 6 A{m, α " 0.02, A " 16ˆ10´1 2 J{m, K u1 " 0.4ˆ10 6 J{m 3 , D 0 " 10´3 m 2 {s, C 0 " 10 6 A{Vm, β " 0.9, β 1 " 0.8, τ sf " 5ˆ10´1 4 s, J " 2.1ˆ10´1 7 J and θ " 0 with the anisotropy axis pointing in z-direction. For the pinned layer, we choose the same parameters, but with a higher anisotropy K u1 " 10 6 J{m 3 . The nonmagnetic spacer and cap layers are simulated with parameters similar to Ag D 0 " 5ˆ10´3 m 2 {s, C 0 " 6ˆ10 6 A{Vm, τ sf " 1.225ˆ10´1 3 s and θ " 0.0.
For the simulation of the write process, the magnetic layers are initialized in positive z-direction. By applying an in-plane electric current in the underlayer, the spin Hall effect gives rise to a spin current in the z-direction. This leads to spin-accumulation in the multilayer stack and thus to spin torque. The torque on the free layer is expected to be much larger than on the pinned layer due its direct neighborhood with the underlayer. Moreover, the pinned layer has a much higher K u1 than the free layer. Hence, the free layer is expected to switch easily while the pinned layer is expected to preserve its magnetization configuration even at high currents.
In a first numerical experiment, we determine the current dependent effective-field contributions, namely the spin accumulation s and the Oersted field H c , for the initial magnetization configuration. We prescribe a constant electric potential at contact 1 (u " 0) and a constant current outflow at contact 2 (j e " 10 12 A{m 2 ). All remaining interfaces are treated with homogeneous Neumann conditions (Bj e {Bn " 0). The resulting fields are plotted in Figure 33. Both the spin accumulation and the Oersted field exhibit a curl-like structure but with opposite sign. In contrast to the Oersted field, the spin accumulation is much smaller in the stack compared to the underlayer, This is due to the interplay of the spin accumulation with the magnetization in the magnetic stack layers.
In the next experiment, we perform time integration in order to resolve the switching of the free layer. In order to enable the spin-orbit torque driven switching we apply an additional external field H zee " p´31830, 0, 0qA{m. Since an electric current in x-direction generates a spin current with polarization y in the z-direction, this additional field is required for the perpendicular switching of the magnetization [152]. Otherwise the spin torque would draw the magnetization toward the xy-plane and the magnetization would return to its initial configuration once the current is turned off.
In a first step, we perform time integration of the LLG (117) including the external field H zee , the exchange field H ex , the demagnetization field H dem and the anisotropy field H ani starting from the initial configuration in order to find an energetic minimum for the system without electric current. After relaxation, we apply a constant current  pulse of 10 12 A{m 2 for the first 0.5 ns that linearly decays to 0 within another 0.5 ns. This pulse is applied in the same fashion as for the field computations, i.e. u " 0 at contact 1 and j e set as Neumann boundary condition on contact 2. In addition to the effective-field contributions considered for the relaxation process, the spin torque due the spin accumulation as well as the Oersted field is included in the simulation. The resulting magnetization dynamics are shown in Figure 34. While the pinned-layer magnetization remains completely fixed in z-direction, the free-layer magnetization performs a fast switch during the current pulse and then relaxes into the´z-direction.
As a final experiment, we simulate the read process of the SOT MRAM. In order to read the magnetization of the free layer, the magnetization dependent resistance of the magnetic multilayer is exploited. When applying a current through the multilayer stack, the resistance of the structure changes either due to giant magnetoresistance (GMR) in the case of a conducting spacer or due to the TMR in the case of an insulating spacer. In both cases, the multilayer is expected to have a lower resistance in case of parallel alignment of the free layer and the pinned layer and a high resistance in case of antiparallel alignment.
For the readout we apply constant potential boundary conditions u " 0 at both contact 1 and contact 2 and set a constant current outflow j e " 10 A{m 2 at contact 3. The pinned-layer magnetization is homogeneously set in z-direction while the free-layer magnetization is homogeneously set to m free " p0, sin θ, cos θq. The coupled system for the spin accumulation and the electric potential is solved for various angles and the potential difference between contact 3 and contact 1/2 is evaluated. Figure 35 shows the simulation results. As expected, the potential, and thus the resistance, of the stack is higher for antiparallel configurations. The difference of the resistance for parallel R p and antiparallel R a configurations is significant pR a´Rp q{R a « 23%. This proves that the presented model includes all crucial physical effects for the self-consistent simulation of both the write process and the read process of an SOT MRAM device.

Conclusion
By bridging the gap between experiment and simple analytical models, computational micromagnetics has proven to be an invaluable tool for the development of magnetic devices. With the rise of spintronics, many extensions to the classical micromagnetic model have been proposed in order to account for the bidirectional coupling of spin-polarized currents and the magnetization configuration. Among the existing models, the spin-diffusion model by Zhang, Levy and Fert is one of the most complete approaches. However, besides other shortcomings, the spin-diffusion model is only valid for diffusive transport which renders this model useless for the accurate description of tunnel barriers. Since the spin-transport through tunnel barriers is dominated by quantummechnical effects, existing models mostly rely on ab initio techniques. The efficient integration of such methods with the micromagnetic model is a challenging yet important task for future research.
The micromagnetic spintronics models introduced in this work already cover a lot of applications. However, detailed knowledge of the applied model and its limitations is crucial for the successful application of micromagnetic simulations. While a complex model might be required in order to understand certain details of magnetization dynamics, for other applications a simpler model might provide sufficient detail and allow to quickly perform a large number of simulations. Moreover, the choice of discretization method has a significant impact on the accuracy and computation time and should be chosen carefully depending on the problem at hand.