A three-dimensional spin-diffusion model for micromagnetics

We solve a time-dependent three-dimensional spin-diffusion model coupled to the Landau-Lifshitz-Gilbert equation numerically. The presented model is validated by comparison to two established spin-torque models: The model of Slonzewski that describes spin-torque in multi-layer structures in the presence of a fixed layer and the model of Zhang and Li that describes current driven domain-wall motion. It is shown that both models are incorporated by the spin-diffusion description, i.e., the nonlocal effects of the Slonzewski model are captured as well as the spin-accumulation due to magnetization gradients as described by the model of Zhang and Li. Moreover, the presented method is able to resolve the time dependency of the spin-accumulation.

An electric or spin current of spin angular momentum acting on a ferromagnet exerts a torque on the magnetization, by doing that driving it out of equilibrium 1,2 . Due to its non-conservative nature this spin torque can act either as effective negative or positive magnetic damping and thereby excite magnetization self-oscillations 3,4 . Spin-torque oscillators (STO) have been realized in nanoscale spin valve systems [5][6][7] ranging from point contact magnetic multilayer structures 8,9 to nanoscale magnetic tunnel junction systems [10][11][12] . Recently, a new type of nano spin-torque oscillator based on current-induced spin orbit torques in a permalloy(Py)/platinum(Pt) bilayer systems was introduced 13,14 . Such spin orbit torques originate from the spin Hall effect in Pt [15][16][17] and the Rashba effect at the Pt/Py interface 18,19 . Applications for STOs include various nanoscale communication technologies. In the upcoming field of microwave assisted recording, a STO is used to generate fields in the GHz regime in order to locally assist the reversal process in hard disks 20 . Another important application of spin torque is to switch magnetic layers in nano structures 21 .
In order to obtain deeper understanding of the applications utilizing spin torque and to guide the design of novel devices, various numerical models have been developed: A model for the description of magnetic multi-layers was first introduced by Slonczewski 1,22 and Berger 23 . The model describes a current flow perpendicular to a layer system consisting of a fixed magnetic layer, a nonmagnetic layer and a free magnetic layer. The current is assumed to pick up its spin polarization in the fixed layer and to exert a torque on the magnetization in the free layer. While this model is successfully applied to MRAM like structures, it is not suited for the description of current driven domain-wall motion, where no magnetic fixed layer is involved. A more general model to spin polarized currents introduces a spin-accumulation field that is bidirectionally coupled to the magnetization 24,25 . The explicit computation of the spin-accumulation can be avoided by neglecting spin-diffusion effects and treating the spin-accumulation in an adiabatic fashion 26 . This simplified model is able to describe domain-wall motion, but fails in the description of multi-layer structures.
In this work, we consider a spin-diffusion model that represents a more general approach to describe the interaction between spin current and magnetization. While similar models were implemented for two-dimensional problems 27 , we consider a three-dimensional space and time regime in this work. While several extensions to this model were proposed, see e.g. 28 , we stick to the original model in this work.
The model is compared to the results of the simplified models of Slonczewski 1,22 and Zhang and Li 25 . For the discretization, we implement a finite-element scheme that was introduced and mathematically analyzed in our previous work 29 . In numerical simulations, we underline that our model generalizes and combines those of 1,22,[24][25][26] in the sense that the same results are also obtained with the more general model, demonstrating the capability to describe and predict more complex features that depend on an inhomogenous spin current distribution.

Let
 ω ⊂ 3 be some ferromagnetic material. The magnetization dynamics is described by the Landau-Lifshitz-Gilbert equation subject to initial and boundary conditions where m is the normalized magnetization, γ is the gyromagnetic ratio, and α is the Gilbert damping. The effective field h eff is given by the negative variational derivative of the total free energy U where μ 0 is the magnetic constant and M s is the saturation magnetization. Contributions to the effective field usually include the exchange field given by with the material dependent exchange constant A and the demagnetization field = −∇ h u demag generated by a magnetic region ω with the potential u defined by with jump and boundary conditions where n is the outer unit normal of ω and, e.g., ⋅ m n is the Euclidean scalar product. Further possible contributions are the external field h ext and the anisotropy field h aniso whose analytical expression depends on the lattice structure of the magnetic material. The vector field s in (1) is the spin-accumulation generated by a spin current J s and c is the corresponding coupling constant. Let  Ω ⊂ 3 be some conducting region with ω ⊂ Ω. According to 30 , the spin-accumulation obeys the equation of motion where D 0 is the diffusion constant, λ sf is related to the spin-flip relaxation time τ sf by τ λ = D 2 sf 0 sf and λ = / hD J 2 J 0 ex with h being Planck's constant and J ex being the exchange integral. For some given current density J e , the matrix-valued spin current J s is given by where μ B is the Bohr magneton, e is the elementary charge and β and β′ are dimensionless polarization parameters. Here,  ( ) 3 3 is the Jacobian of s and  ⊗ = ∈ × a b ab T 3 3 is the outer product. Figure 1 shows a typical multi-layer structure for the investigation of spin-torque effects. A magnetic fixed layer ω fixed is separated from a magnetic free layer ω free by a nonmagnetic layer ω nm . Additional nonmagnetic layers above and beneath the magnetic layers serve as electrical contacts. The Landau-Lifshitz-Gilbert equation (1)

Algorithm
The model introduced in the preceding section is discretized by an algorithm, we proposed and analyzed mathematically in 29 . The vector fields m, s, = ∂ /∂ v m t, and test functions w, ζ are discretized with piecewise affine, globally continuous functions constructed on a tetrahedral tessellation of the domain Ω . Material constants that differ from subdomain to subdomain, e.g., the diffusion constant D 0 , are discretized by piecewise constant functions. A splitting scheme is used to decouple the Landau-Lifshitz-Gilbert equation (1) from the spin-accumulation equation (10).
Suppose that the magnetization m k at the k-th timestep is known. Then, we use the tangent-plane scheme of 31 to compute m k+1 at the next time-step, where the essential idea is to discretize ∂ /∂ m t at the k-th timestep by some v k which mimics the continuous orthogonality ⋅ (∂ /∂ ) = m m t 0. Following 32-34 , all terms but the exchange term are treated explicitly. This also includes the spin-accumulation s. Based on m k+1 and s k , we then use the implicit Euler scheme to compute the spin-accumulation s k+1 at the next time-step. The complete algorithm for a timestep size τ > 0 reads as follows, where we suppose that m k and s k are already known: 1. Consider the discrete tangent space  m k to the known magnetization m k defined by is the Frobenius inner product of two matrices  , ∈ × A B 3 3 . While the exchange field ist integrated implicitly, all other effective-field contributions are summarized in h explicit and integrated explicitly. 2. Define the piecewise affine and globally continuous approximation m k+1 by ∪ ω ω Ω = nm . A magnetic fixed layer ω fixed is separated from a magnetic free layer ω free by a nonmagnetic spacer layer. The system is completed by two leads made of nonmagnetic material.
where n is the outer normal of ω and Update the timestep and iterate.
While our preceding work 29 did not include any numerical simulations with the proposed time-marching scheme, we proved that the coupled system of LLG with spin-diffusion admits weak solutions and that the proposed integrator converges to these in an appropriate sense. With the present work, we aim to provide empirical evidence for the physical relevance of the proposed model and the numerical results that are obtained by the described discretization.

Implementation
The presented algorithm is implemented within the finite-element micromagnetic code magnum. fe 35,36 . magnum.fe uses FEniCS 37 for the assembly and solution of the finite-element systems. The linear systems arising from discretization of (14) and (16) are solved iteratively by a GMRES solver. The demagnetization-field problem is solved by a hybrid finite-element boundary-element method 38 , where our implementation couples FEniCS to the boundary-element library BEM+ + 39 . Each timestep requires the solution of two linear systems of the dimension 3N × 3N, where N is the number of mesh nodes. The finite-element discretization leads to sparse matrices, whose storage requirements scale linearly with the number of mesh nodes. Hence systems with 100000 nodes and more may be handled in reasonable time by magnum.fe, which corresponds to system sizes up to the m regime depending on the spatial discretization.

Numerical Experiments
As a first test the time development of the spin diffusion for a fixed magnetization configuration is computed. An ellipsoidal multi-layer structure as pictured in Fig. 1 with axis lengths 130 × 70 nm and layer thicknesses (5, 2, 3, 2, 5) nm is considered. In the following the coordinate system is chosen such that the long and short axes of the ellipse align with the x and y-axis respectively. The diffusion constant in the magnetic layers ω free and ω fixed is set to D 0 = 1 × 10 −3 m 2 /s. In the nonmagnetic region the diffusion constant is chosen as D 0 = 5 × 10 −3 m 2 /s. The remaining constants are λ sf = 10 nm, λ J = 2 nm, β = 0.9, β′ = . . The magnetization m is initialized homogeneously parallel to the x-axis in the fixed layer and parallel to y-axis in the free layer. Equation (16)  , a time step τ = 1 fs and an initial spin-accumulation of s = 0. Figure 2(a) shows the simulation results. The system relaxes within approximately 70 fs. This time scale is 2 orders of magnitude below the typical reaction time of the magnetization m. Hence the magnetization dynamics may be accurately described by assuming the spin-accumulation to be in equilibrium at any times 24 . However, the presented integration scheme for the spin-accumulation is able to handle very large time steps due to its implicit nature and the linearity of the problem. In fact it shows that even for an initial spin-accumulation of s = 0 two integration steps with τ = 1 ps are sufficient to obtain the equilibrium state accurately. When moving from one magnetization configuration to the next during time integration of the Landau-Lifshitz-Gilbert equation, a single integration step of the spin-accumulation yields the equilibrium state. The presented algorithm is thus well suited for investigating the time evolution of both the spin diffusion using time steps τ ≈ 1fs as well as the magnetization using time steps τ ≈ 1ps. Figure 2(b) shows the equilibrium spin-accumulation s on the z-axis in the middle of the multi-layer structure. From (10) it is clear, that the source of the spin-accumulation is the spatial change of the magnetization m. Since the magnetization is homogeneous within the magnetic regions, the spin-accumulation where q, A, and B are constants that depend on the geometry and the material parameters of the system. Consider the previously described ellipsoidal multi-layer structure with the magnetization being homogeneous and in-plane in both the fixed layer and the free layer. In this case the additional torque term in the free layer (18) is expected to be in-plane with a magnitude given by Furthermore the in-plane magnitude is fitted with the magnitude as predicted by Slonczewski (20). Since the system is overdefined for fitting, q is set to 1 and A and B are chosen as fitting parameters. As seen in Fig. 3(a) the fit is very accurate. The relative fitting error for both A and B is well below 10 −6 . Not only the in-plane magnitude, but also the in-plane angle of the torque term shows excellent accordance to the spin-torque model of Slonczewski as depicted in Fig. 3(b). However, the spin-accumulation model predicts a non-neglible out-of-plane torque that cannot be explained by the simplified model. As pointed out in 40 the expression for the angular coefficient (19) only holds for the symmetric case of equal height fixed and free layer. For the general case the angular coefficient reads  a z-component of the torque that is not described by the spin-torque model of Slonczewski. The ratio of in-plane torque to out-of-plane torque for a given structure is largely influenced by the choice of λ J .
Choosing large values for λ J leads to high out-of-plane components of the torque, see Fig. 4(b). The coefficients A, B, q, q + , and q − depend, among other things, on the geometry of the system. In the following the dependence of these constants on the thickness of the nonmagnetic layer and the fixed layer is investigated in detail. A physically motivated choice for the right-hand-side coefficients in (21) is introduced in 40 where Λ L , Λ R , P L , and P R describe the properties of the fixed and the free layer respectively. For symmetric structures, i.e. Λ = Λ L R and P L = P R , (21) reduces to  .
The dependence of Λ L and Λ R on the thickness of the fixed layer in an asymmetric structure is depicted in Fig. 5(b). Both coefficients are well fitted by a shifted reciprocal square function with X, Y, and Δ being fit parameters. The model of Slonczewski provides a straightforward way to deal with multi-layer structures. However, it is obviously not suited for the description of spin-torque effects in smoothly variing magnetization configurations since it requires a separated fixed layer to generate the spin polarization of the current.  relative L 2 error norm below 5 × 10 −3 . The assumption of a vanishing gradient of the spin-accumulation ∇s, made in the model by Zhang and Li, is valid within the magnetic regions as discussed in 26 . However, at the interfaces of magnetic multi-layers this assumption is violated. Hence the model by Zhang and Li is not suited for the description of such structures.

Standard Problem #5. A suitable benchmark for the interplay of magnetization and spin-diffusion
dynamics is the micromagnetic standard problem #5 that was recently proposed by the Mag group, see 41 .
A thin square of size 100 × 100 × 10 nm with material parameters similar to permalloy is prepared in a vortex state. Then a constant current J e = 10 12 A/m is applied in x-direction. The vortex core is expected to perform a damped rotation around a new equilibrium position. The spin-torque related costants in the problem definition are tailored to the model by Zhang and Li. Namely the degree of non-adiabacity ξ is given as well as a factor b that can be expressed in terms of spin-diffusion related constants as / bJ m 72 17 s is exactly fulfilled as required. Figure 7 shows the simulation results in comparsion to a reference solution computed with the finite-difference code MicroMagnum 42 that implements the model by Zhang and Li. Additionally the averaged magnetization in the new equilibrium computed with different simulation tools is compared in Table 1. The most notable difference is the shifted y-position of the vortex core in equilibrium. This aberration is a result of the missing ∇s terms in the model by Zhang and Li. Note that the validity of the spin-diffusion model for Permalloy vortex structures is controversial. The derivation of this model assumes a short mean free path compared to the domain-wall width which is violated in this case 43 . However, the purpose of this numerical experiment is to show that the presented model reproduces the results by the model of Zhang and Li and the standard problem #5 is a widely accepted test problem for current driven domain-wall motion.
Switching of a multi-layer Structure. As another test the current induced switching of a magnetic multi-layer structure is simulated. The geometry considered is ellipsoidal with edge lengths 130 × 70 nm and layer thicknesses (3, 10 . Since the size of the ellipsoid is well under the single-domain limit, the energy of the system is minimal if the magnetization is aligned along the long axis of the ellipsoid in the magnetic layers. The thick and thin magnetic layers are referred to as fixed and free layer respectively, since the energy barrier for magnetic switching in the thick layer is higher than in the thin layer. By applying a current perpendicular to the layer structure, the magnetization in the free layer can be switched. Figure 8(a) shows the hysteresis of the free layer for different currents. The shift of the curve is an expected consequence resulting from the strayfield coupling of the magnetic layers. An antiparallel configuration of the fixed and free layer is energetically favored over a parallel configuration. Figure 8(b) shows the magnetization configuration of the free layer during switching. Note that the out-of plane component is small compared to the in-plane component, which is a consequence of the shape anisotropy caused by the small thickness of the layer. The out-of-plane component of the torque as predicted by the spin-diffusion model is compensated by the shape anisotropy as soon as the magnetization is driven slightly out of plane. This may serve as an justification for the missing out-of-plane spin-torque component in the model of Slonczewski.

Conclusion
The presented method for the solution of the Landau-Lifshitz-Gilbert equation including spin-torque effects is able to describe both the spin transport in multi-layer structures as well as current driven domain-wall motion. The spin-torque in multi-layers as predicted by Slonczewski is essentially reproduced except for an out-of-plane component that is not present in the model by Slonczewski. However, the out-of-plane component is shown to play an inferior role in typical multi-layer structures since it is usually almost completely compensated by shape anisotropy effects. The presented method only requires the geometry and material parameters of the multilayer system as input. In contrast the model of Slonczewski requires some global parameters that depend on the system as a whole and that have to be determined by experiment in order to perform meaningful simulations. For the description of current-driven domain-wall motion the presented method is compared to a simplified model by Zhang and Li and shows a good agreement. Small deviations can be explained by diffusion terms that are omitted in the simplified model. While the present method is able to resolve the temporal evolution of the spin-accumulation exactly it can also be used for an adiabatic description since the time integrator for the spin-accumulation behaves very well even for large time steps.