A self-consistent spin-diffusion model for micromagnetics

We propose a three-dimensional micromagnetic model that dynamically solves the Landau-Lifshitz-Gilbert equation coupled to the full spin-diffusion equation. In contrast to previous methods, we solve for the magnetization dynamics and the electric potential in a self-consistent fashion. This treatment allows for an accurate description of magnetization dependent resistance changes. Moreover, the presented algorithm describes both spin accumulation due to smooth magnetization transitions and due to material interfaces as in multilayer structures. The model and its finite-element implementation are validated by current driven motion of a magnetic vortex structure. In a second experiment, the resistivity of a magnetic multilayer structure in dependence of the tilting angle of the magnetization in the different layers is investigated. Both examples show good agreement with reference simulations and experiments respectively.

eff s  where m is the normalized magnetization, γ is the gyromagnetic ratio, α is the Gilbert damping, and h eff is the effective field that usually contains the demagnetization field, the exchange field, as well as other contributions depending on the problem setting. The effective field is complemented by a contribution from the spin accumulation s with J being the exchange strength between itinerant and localized spins, ħ being the reduced Planck constant, and M s being the saturation magnetization. The spin accumulation itself is defined in the conducting region Ω and satisfies the equation of motion where τ sf is the spin-flip relaxation time, and j s is the matrix-valued spin current. According to refs 9 and 11, the spin current j s and the electric current j e are defined by where E is the electric field, D 0 is a diffusion constant, C 0 is related to electric resistivity ρ by C 0 = 1/2ρ, and β and β′ are dimensionless polarization parameters. Solving (3) for E and inserting this into (4) yields which, combined with (2), gives the simplified diffusion model with prescribed electric current j e used in refs 8, 12 and 13. However, instead of prescribing the electric current, it is possible to solve the coupled system (3) and (4). For this purpose, we consider the following simplifications: As proposed in ref. 13, we assume the spin accumulation to be in equilibrium at all times, i.e., This simplification is justified by the fact that the characteristic time scale of the spin accumulation dynamics is two orders of magnitude smaller than that of the magnetization dynamics. This was predicted theoretically in ref. 7 and shown by time resolved solution of (2) for typical material parameters in ref. 8. Since sample sizes are usually very small, eddy currents can be neglected 14 . Therefore, the electric field is curl free and thus given as the gradient of a scalar potential Moreover, it is assumed that the conducting region Ω, that is considered for the solution of the system (3) and (4), does not contain any sources of electric currents, i.e., Inserting these assumptions into (2)-(4) yields the overall system The system (9)-(10) introduces the need for further boundary conditions for the electric potential u and the spin accumulation s. A set of mixed boundary conditions is used to prescribe the electric potential and current inflow at the boundary of the conducting region ∂Ω. The Dirichlet condition is applied directly to the potential u A typical choice of these boundary conditions is depicted in Fig. 1(b). In an MRAM like structure, the top and bottom surfaces Γ 1 and Γ 2 are electric contacts. In order to prescribe the current flow through the sample, like it is usually done in experiments, one might set the potential to zero at one contact u = 0 on Γ 1 . On the other contact Γ 2 a finite current inflow is prescribed j e · n = g. The rest of the sample boundary ∂Ω\(Γ 1 ∪Γ 2 ) is treated with homogeneous Neumann conditions j e · n = 0.
The boundary conditions are completed with homogeneous Neumann conditions for the spin accumulation ∇ ⋅ = ∂Ω. s n 0 on (15) This choice can be justified by considering the boundary flux of the spin current. Multiplying (5) with the boundary normal n and inserting the Neumann condition yields This expression is nonzero only at boundaries with both nonvanishing current in-/outflow and nonvanishing magnetization. Hence the homogeneous Neumann condition (15) is equivalent to the noflux condition j s · n = 0 for systems as depicted in Fig. 1, where the electric contacts are part of the nonmagnetic region. The noflux condition itself is a reasonable choice if the thickness of the electrodes is large against the diffusion length. In this case the spin accumulation and hence also the spin current is expected to be approximately zero at the contacts.
For systems where the magnetic region is contacted directly, the homogeneous Neumann condition leads to unphysical behaviour since the spin accumulation that is generated at the contact interface is neglected. This accumulation strongly depends on the material of the leads that is not known when directly contacting the magnet. However, the choice of homogeneous Neumann conditions leads to a good agreement with the predictions of the model by Zhang and Li 7 that also neglects surface effects at the contacts.

Validation
The presented model is implemented within the finite-element code magnum.fe 16 . The discretization is explained in detail in Appendix A. For validation purposes, the standard problem #5 proposed by the μMAG group 17 is computed with the self-consistent model and compared to results obtained with the model of Zhang and Li 7 and the simplified diffusion model used in ref. 8. While this problem does not require particular features of the proposed self-consistent model, it serves as an excellent experiment for the validation of the proposed algorithm. The standard problem #5 describes the motion of a magnetic vortex in a thin square of size 100 nm × 100 nm × 10 nm under the influence of a DC current defined by βj e = (10 12 , 0, 0)A/m 2 . For our simulations we choose β = 1 and thus j e = (10 12 , 0, 0) A/m 2 . The material parameters are chosen similar to those of permalloy, namely M s = 8 × 10 5 A/m, A = 1.3 × 10 −11 J/m, and α = 0.1. In the original problem definition, it is proposed to apply the model of Zhang and Li that extends the LLG (1) by current dependent terms

Resistance of multilayer stack with perpendicular current
The key advantage of the presented method over the simplified diffusion model introduced in ref. 8 is the self-consistent treatment of the electric potential u. The potential is computed considering not only Ohm's law j e = 2C 0 E but also magnetization dependent contributions. This dependence is exploited for instance in sensor applications 1 . Consider a magnetic multilayer stack as shown in Fig. 1 with two homogeneously magnetized layers ω 1 and ω 2 separated by a conducting layer. The resistivity of this structure heavily depends on the tilting angle θ of the magnetization in ω 1 and ω 2 . Namely an antiparallel configuration is known to result in a high resistivity while a parallel configuration leads to a low resistivity. This effect is referred to as giant magnetoresistance (GMR). The anisotropic magnetoresistance (AMR), which describes the change of electrical resistance as a function of the angle between electric current and the magnetization, is not considered in the presented model. In order to calculate the electric resistivity with the diffusion model, the potential difference between the contacts Γ 1 and Γ 2 is computed for a given current inflow. Note, that due to the homogeneous magnetization configuration within the magnetic layers, the simulation is quasi one-dimensional and the results do not depend on the lateral dimension of the stack.
For numerical experiments two magnetic layers with 5 nm thickness separated by a conducting layer with 1.5 nm thickness are considered. The system is contacted with 100 nm thick leads in order to justify the homogeneous boundary conditions on the spin accumulation as described in the model section. The potential is set to u = 0 at the bottom contact Γ 1 and the current inflow is set to j e · n = 10 12 A/m 2 on the top contact Γ 2 . Note that the cross section of the system does not have any influence on the potential computation as long as it is constant throughout the stack. Figure 3 shows the computed potential difference for different tilting angles of the magnetization in the two layers ω 1 and ω 2 . The material parameters in the magnetic regions are chosen similiar to those in the preceding section. In the conducting region Ω\ω, a different diffusion constant of D 0 = 5 × 10 −3 m 2 /s and a conductivity of C 0 = 6.0−10 6 A/(Vm) is applied. Moreover, in Fig. 3(a), exchange strength J is varied in the whole stack Ω. In Fig. 3(b), the polarization parameter β′ is varied. The resulting potential is compared to a sine parameterization a + b sin 2 (θ/2) that is often used to describe the GMR effect in such a stack 18 in the presence of some in-plane current. The presented simulations, however, suggest that the potential and thus the resistivity of the stack in the presence of out-of-plane currents is not well described by a sine, but has a much narrower peak for certain choices of material parameters. Specifically the sine approximation is accurate only in the case of small β′ and J as shown in Fig. 4(a), where the parameters are chosen as β′ = 0.1 and J = 0.013 eV.
The deviation of the angular dependence of the resistivity from the sine approximation for some out-of-plane current was already observed in experiment 19 . The authors of this work suggest to include a higher order sine term in order to fit the resistivity curve. Figure 4(b) shows the result of the fit with a + b sin 2 (θ/2) + c sin 4 (θ/2) for J = 0.082 eV and β′ = 0.8 which shows a good agreement with the simulated curve. A similar effect is also predicted in ref. 6, where the authors suggest the following expression for the angular dependence of the resistivity The parameter χ depends on geometry and material of the involved layers and leads to a steeper peek if positive, which, according to the reference, is the case for all investigated systems up to then. Using χ as a fitting parameter, the results for J = 0.082 eV and β′ = 0.8 can be described with very high precision, see Fig. 4(b).
The same effect was also predicted theoretically in ref. 20 with a model of Valet and Fert 21 and in ref. 22 with a two-dimensional diffusion model. However, these papers do not discuss the influence of material parameters onto the narrowing of the sine response in detail.
In the context of the diffusion model, the deviation from the simple sine approximation sin 2 (θ/2) has its origin in two different effects. First, the cross product term Js × m/ħ in (10) describes the torque that is exerted from the magnetization on the spin polarization of the itinerant electrons. This torque is zero for parallel and antiparallel alignment of the magnetic layers and reduces the angle of the polarization of the itinerant electrons to the magnetization for any other alignment. Hence, for large J this contribution leads to a lowered resistance for all alignments other than parallel and antiparallel, which results in the narrow peak observed in the simulation.
The second effect is a bit more subtle. For vanishing β′ the potential u from (9) varies linearly. Small values of β′ lead to small perturbations of this linear solution. While these perturbations have a clear effect on the overall potential difference, their effect on the spin accumulation s due to (10) is negligible, leading to a clean sinusoidal response of the system as shown in Fig. 4(a). With increasing β′ the perturbations of u gain influence on the solution of s which results in a distorsion of the sinusoidal response as seen in Figs 3 and 4.

Conclusion
We propose a three-dimensional spin-diffusion model that simultaneously solves for the spin accumulation s and the electric potential u. By coupling this model to the Landau-Lifshitz-Gilbert equation, we are able to self-consistently solve the magnetization dynamics for a given current inflow. In order to validate the model and its implementation, we simulate the standard problem #5 proposed by the μMAG group and compare the outcome to results obtained with simplified models.
In a second numerical experiment, we compute the resistivity of a magnetic multilayer structure in dependence on the tilting angle of the magnetization in the two magnetic layers. In the limit of small polarization parameter β′ and a small exchange strength J, we show that the resistivity is well approximated by a sine. For realistic choices of β′ and J the angular dependence shows a significantly narrower peak than the simple sine approximation. While existing models already predict the observed behaviour in a macro spin approach, the presented model is able to accurately describe GMR effects for both dynamically and spatially varying magnetization configurations.
While the examples in this work consider either spin-torque effects or magnetoresistance, we want to stress that the presented model is also capable of simulating the interplay of the two effects. This interplay is for instance crucial for the description of spin-torque induced noise in GMR sensors that was found to be a significant contribution to the overall noise 23 .

A Discretization
We solve the coupled equations (1) and (9)-(10) numerically by employing the finite-element method for spatial discretization and a preconditioned BDF scheme as described in ref. 24 for the time integration that is applied nodewise. The BDF scheme is an implicit method and hence well suited for the solution of stiff problems. This stiffness of micromagnetic problems usually originates from the exchange interaction and the additional solution of the diffusion equation does not influence the stiffness noticeably. Typical time steps for this method are 10 −13 s to 10 −12 s. The demagnetization field is solved by a hybrid FEM-BEM method as introduced in ref. 25.
For the discretization of (9) and (10) special care has to be taken in the treatment of the discontinuities, e.g., in the magnetization m introduced by magnetic-nonmagnetic interfaces. As usual the original problem is multiplied with test functions and integration by parts is applied to avoid second derivatives. Furthermore, first derivatives of the magnetization m are eliminated in the same way and the integration domain for terms including the magnetization is restricted to the magnetized region ω. This way, the proposed algorithm naturally accounts for the magnetization jumps at the magnetic-nonmagnetic interfaces without the need of discontinuous function spaces for the discretization of m.
The solution variables m, u and s as well as the test functions v and ζ are discretized by (componentwise) piecewise affine, globally continuous functions constructed on a tetrahedral mesh. The material parameters β, β′, C 0 , D 0 , τ sf , and J are discretized with piecewise constant functions. For a given magnetization m, the weak version of (9) reads