A self-consistent spin-diffusion model for micromagnetics
Claas Abert, Michele Ruggeri, Florian Bruckner, Christoph Vogler, Aurelien Manchon, Dirk Praetorius, Dieter Suess
AA self-consistent spin-diffusion model for micromagnetics
Claas Abert ∗ , Michele Ruggeri , Florian Bruckner , Christoph Vogler , AurelienManchon , Dirk Praetorius , and Dieter Suess Christian Doppler Laboratory of Advanced Magnetic Sensing and Materials, Institute of Solid StatePhysics, TU Wien, Austria Institute for Analysis and Scientific Computing, TU Wien, Austria Institute of Solid State Physics, TU Wien, Austria Physical Science and Engineering Division, King Abdullah University of Science and Technology(KAUST), Thuwal 23955-6900, Kingdom of Saudi Arabia
October 13, 2018
Abstract
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 pre-vious methods, we solve for the magnetization dynamics and the electric potential in aself-consistent fashion. This treatment allows for an accurate description of magnetizationdependent resistance changes. Moreover, the presented algorithm describes both spin ac-cumulation due to smooth magnetization transitions and due to material interfaces as inmultilayer structures. The model and its finite-element implementation are validated bycurrent driven motion of a magnetic vortex structure. In a second experiment, the resistiv-ity of a magnetic multilayer structure in dependence of the tilting angle of the magnetizationin the different layers is investigated. Both examples show good agreement with referencesimulations and experiments respectively. ∗ [email protected] a r X i v : . [ phy s i c s . c o m p - ph ] S e p Introduction
Spin-tronic devices are versatile candidates for a variety of applications including sensors [1, 2],storage devices [3], and frequency generators [4, 5]. Different quantum mechanical mechanismscontribute to the coupling of the electrical current to the magnetization. Simulations of spin-tronic systems usually apply the micromagnetic model extended by simplified coupling termssuch as the model by Slonczeswki [6] and the model by Zhang and Li [7]. In [8] it was shownthat a simplified spin diffusion model incorporates these models. In all these approaches, theeffects of the magnetic state of the system on the electronic transport are neglected. Indeed,the electric current density is assumed to be fixed, so that the models can be only used toinvestigate how the electronic transport affects the magnetization dynamics, but not vice versa.In this work, we present a three-dimensional finite-element method for the solution of the fullspin diffusion model from [9] which includes a bidirectional coupling of the magnetization tothe electric current. In [10] the model from [9] is applied to a single phase magnetic nanostruc-ture in order to predict domain wall motion. In this work we extend this model to compositestructures consisting of magnetic and nonmagnetic materials which enables us to compute themagnetization-dependent resistivity caused by the GMR effect.
According to the micromagnetic model, the magnetization dynamics in a three-dimensionalmagnetic domain ω is described by the Landau-Lifshitz-Gilbert equation (LLG) ∂ m ∂t = − γ m × (cid:18) h eff + J (cid:126) γM s s (cid:19) + α m × ∂ m ∂t in ω, (1)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 exchangefield, as well as other contributions depending on the problem setting. The effective field iscomplemented by a contribution from the spin accumulation s with J being the exchangestrength between itinerant and localized spins, (cid:126) being the reduced Planck constant, and M s being the saturation magnetization. The spin accumulation itself is defined in the conductingregion Ω and satisfies the equation of motion ∂ s ∂t = − ∇ · j s − s τ sf − J s × m (cid:126) in Ω , (2) (a) ω ω ω Ω (b) Γ Γ Figure 1: Typical magnetic/nonmagnetic material stack as used for MRAM devices. (a) Theregion ω = ω ∪ ω denotes the magnetic material, while Ω denotes the complete sample. (b)Electrical contacts Γ and Γ . 2here τ sf is the spin-flip relaxation time, and j s is the matrix-valued spin current. Accordingto [9, 11], the spin current j s and the electric current j e are defined by j e = 2 C E − β (cid:48) D eµ B ( ∇ s ) T m , (3) j s = 2 βC µ B e m ⊗ E − D ∇ s , (4)where E is the electric field, D is a diffusion constant, C is related to electric resistivity ρ by C = 1 / ρ , and β and β (cid:48) are dimensionless polarization parameters. Solving (3) for E andinserting this into (4) yields j s = β µ B e m ⊗ j e − D (cid:2) ∇ s − ββ (cid:48) m ⊗ (( ∇ s ) T m ) (cid:3) , (5)which, combined with (2), gives the simplified diffusion model with prescribed electric current j e used in [12, 13, 8].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 [13], weassume the spin accumulation to be in equilibrium at all times, i.e., ∂ s ∂t = 0 . (6)This simplification is justified by the fact that the characteristic time scale of the spin accumula-tion dynamics is two orders of magnitude smaller than that of the magnetization dynamics, see[7]. 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 E = − ∇ u. (7)Moreover, it is assumed that the conducting region Ω, that is considered for the solution of thesystem (3) and (4), does not contain any sources of electric currents, i.e., ∇ · j e = 0 . (8)Inserting these assumptions into (2) – (4) yields the overall system − ∇ · (cid:20) C ∇ u + β (cid:48) D eµ B ( ∇ s ) T m (cid:21) = 0 , (9) ∇ · (cid:104) βC µ B e m ⊗ ∇ u + 2 D ∇ s (cid:105) − s τ sf − J s × m (cid:126) = 0 , (10)that determines the electric potential u and the spin accumulation s .A number of boundary conditions are required in order to close the LLG (1) coupled to thespin-diffusion system for the magnetization m ( t ). The LLG itself is an initial value problemand requires the initial magnetization m ( t ) = m . (11)If the exchange field h exchange = 2 A/ ( µ M s )∆ m , with the exchange constant A and the sat-uration magnetization M s , is included in the list of effective field contributions, an additionalboundary condition for the magnetization is required. If the domain ω for the solution of the3LG coincides with the magnetic region, it was shown in [15] that homogeneous Neumannboundary conditions are the right choice in a physical sense ∂ m ∂ n = 0 on ∂ω. (12)The system (9) – (10) introduces the need for further boundary conditions for the electricpotential u and the spin accumulation s . A set of mixed boundary conditions is used to prescribethe electric potential and current inflow at the boundary of the conducting region ∂ Ω. TheDirichlet condition is applied directly to the potential uu = u on Γ D ⊆ ∂ Ω , (13)while the Neumann condition is applied to the electric current j e · n = − (cid:20) C ∇ u + β (cid:48) D eµ B (cid:2) ( ∇ s ) T m (cid:3)(cid:21) · n = g on Γ N = ∂ Ω \ Γ D . (14)A typical choice of these boundary conditions is depicted in Fig. 1(b). In an MRAM likestructure, the top and bottom surfaces Γ and Γ are electric contacts. In order to prescribethe current flow through the sample, like it is usually done in experiments, one might set thepotential to zero at one contact u = 0 on Γ . On the other contact Γ a finite current inflowis prescribed j e · n = g . The rest of the sample boundary ∂ Ω \ (Γ ∪ Γ ) is treated withhomogeneous Neumann conditions j e · n = 0.The boundary conditions are completed with homogeneous Neumann conditions for the spinaccumulation ∇ 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 j s · n = β µ B e m ( j e · n ) − D (cid:2) ∇ s · n − ββ (cid:48) m (( ∇ s · n ) · m ) (cid:3) (16)= β µ B e m ( j e · n ) (17)This expression is nonzero only at boundaries with both nonvanishing current in-/outflow andnonvanishing magnetization. Hence the homogeneous Neumann condition (15) is equivalent tothe noflux condition j s · n = 0 for systems as depicted in Fig. 1, where the electric contacts arepart of the nonmagnetic region. The noflux condition itself is a reasonable choice if the thicknessof the electrodes is large against the diffusion length. In this case the spin accumulation andhence also the spin current is expected to be approximately zero at the contacts.For systems where the magnetic region is contacted directly, the homogeneous Neumanncondition leads to unphysical behaviour since the spin accumulation that is generated at thecontact interface is neglected. This accumulation strongly depends on the material of the leadsthat is not known when directly contacting the magnet. However, the choice of homogeneousNeumann conditions leads to a good agreement with the predictions of the model by Zhang andLi [7] that also neglects surface effects at the contacts. The presented model is implemented within the finite-element code magnum.fe [16]. The dis-cretization is explained in detail in Appendix A. For validation purposes, the standard prob-lem µ MAG group [17] is computed with the self-consistent model and4 a) − . − . − . t [ns] h m x i Zhang and Lispin diffusionself consistent (b) − . . . t [ns] h m y i Zhang and Lispin diffusionself consistentFigure 2: Results for the standard problem x -component (b) y -component.compared to results obtained with the model of Zhang and Li [7] and the simplified diffu-sion model used in [8]. While this problem does not require particular features of the pro-posed self-consistent model, it serves as an excellent experiment for the validation of the pro-posed algorithm. The standard problem ×
100 nm ×
10 nm under the influence of a DC current defined by β j e = (10 , ,
0) A/m . For our simulations we choose β = 1 and thus j e = (10 , ,
0) A/m .The material parameters are chosen similar to those of permalloy, namely M s = 8 × A/m, A = 1 . × − J/m, and α = 0 .
1. In the original problem definition, it is proposed to applythe model of Zhang and Li that extends the LLG (1) by current dependent terms ∂ m ∂t = − γ m × h eff + α m × ∂ m ∂t − b m × [ m × ( j e · ∇ ) m ] − ξb m × ( j e · ∇ ) m , (18)where b = 72 . × − m /(As) is a coupling constant and ξ = 0 .
05 the degree of nonadi-abacity. As shown in [8] an equivalent set of material parameters for the diffusion model canbe obtained by perceiving the Zhang and Li model as diffusion model in the limit of vanishingdiffusion. For the diffusion model we choose D = 10 − m /s, β (cid:48) = 0 .
8, and τ sf = 5 × − s.The remaining constant J = 0 .
263 eV is then uniquely defined by the relations given in [8]. Inorder to solve this problem with the self-consistent model the additional conductivity constant C = 1 . × A/(Vm) is introduced. Furthermore, instead of prescribing a constant currentwithin the magnetic material, the current is applied in terms of boundary conditions. On theleft side of the sample Γ N = { r | r x = −
50 nm } current inflow is set to j e · n = 10 A/m andon the right side of the sample Γ D = { r | r x = 50 nm } the potential is set to 0. The remainingboundary is treated with homogeneous Neumann conditions in order to simulate current in-and outflow only through the contacts.The results for the computation of standard problem x -component of the magnetization arein very good agreement, the results for the y -component show a notable offset. The offset ofthe results of the Zhang and Li model to the remaining models is caused by the neglecteddiffusion. The offset of the self-consistent model to the simple diffusion model is caused by theinhomogeneous current distribution resulting from the self-consistent treatment.5 a) θ [ ° ] ∆ u [ m V ] sin (cid:16) θ (cid:17) J = 1 .
316 eV J = 0 .
329 eV J = 0 .
013 eV (b) θ [ ° ] ∆ u [ a . u . ] sin (cid:16) θ (cid:17) β = 0 . β = 0 . β = 0 . Figure 3: Potential difference ∆ u between top and bottom contact required to generate anaverage current density of j e = 10 A/m depending on the tilting angle of the magnetizationin the free and fixed layer of a magnetic multilayer structure. (a) variation of J for β (cid:48) = 0 . β (cid:48) for J = 0 .
263 eV. The results are renormalized in order to facilitate thecomparison to the sine parameterization.
The key advantage of the presented method over the simplified diffusion model introduced in[8] is the self-consistent treatment of the electric potential u . The potential is computed con-sidering not only Ohm’s law j e = 2 C E but also magnetization dependent contributions. Thisdependence is exploited for instance in sensor applications [1]. Consider a magnetic multilayerstack as shown in Fig. 1 with two homogeneously magnetized layers ω and ω separated by aconducting layer. The resistivity of this structure heavily depends on the tilting angle θ of themagnetization in ω and ω . Namely an antiparallel configuration is known to result in a highresistivity while a parallel configuration leads to a low resistivity. This effect is referred to asgiant magnetoresistance (GMR). In order to calculate the electric resistivity with the presentedmodel, the potential difference between the contacts Γ and Γ is computed for a given currentinflow.For numerical experiments two magnetic layers with 5 nm thickness separated by a conduct-ing layer with 1 . u = 0 at the bottom contact Γ and the current inflow is setto j e · n = 10 A/m on the top contact Γ . Note that the cross section of the system does nothave any influence on the potential computation as long as it is constant throughout the stack.Fig. 3 shows the computed potential difference for different tilting angles of the magne-tization in the two layers ω and ω . The material parameters in the magnetic regions arechosen similiar to those in Sec. 3. In the conducting region Ω \ ω , a different diffusion constantof D = 5 × − m /s and a conductivity of C = 6 . × A/(Vm) is applied. Moreover,in Fig. 3(a), exchange strength J is varied in the whole stack Ω. In Fig. 3(b), the polariza-tion parameter β (cid:48) is varied. The resulting potential is compared to a sine parameterization a + b sin ( θ/
2) that is often used to describe the GMR effect in such a stack [18] in the presenceof some in-plane current. The presented simulations however suggest that the potential and thusthe resistivity of the stack out-of-plane currents is not well described by a sine, but has a muchnarrower peak for certain choices of material parameters. Specifically the sine approximation6 a) . . θ [ ° ] ∆ u [ m V ] sin (cid:16) θ (cid:17) simulation (b) θ [ ° ] ∆ u [ m V ] sin (cid:16) θ (cid:17) sin , sin χ cos modelsimulation Figure 4: Potential difference ∆ u between top and bottom contact required to generate anaverage current density of j e = 10 A/m depending on the tilting angle of the magnetizationin the free and fixed layer of a magnetic multilayer structure along with fits to different modelsfor (a) J = 0 .
013 eV and β (cid:48) = 0 . J = 0 .
082 eV and β (cid:48) = 0 . β (cid:48) and J as shown in Fig. 4(a), where the parameters arechosen as β (cid:48) = 0 . J = 0 .
013 eV.The deviation of the angular dependence of the resistivity from the sine approximation forsome out-of-plane current was already observed in experiment [19]. The authors of this worksuggest to include a higher order sine term in order to fit the resistivity curve. Fig. 4(b) showsthe result of the fit with a + b sin ( θ/
2) + c sin ( θ/
2) for J = 0 .
082 eV and β (cid:48) = 0 . R = 1 − cos ( θ/ χ cos ( θ/ . (19)The parameter χ depends on geometry and material of the involved layers and leads to a steeperpeek if positive, which, according to the reference, is the case for all investigated systems up tothen. Using χ as a fitting parameter, the results for J = 0 .
082 eV and β (cid:48) = 0 . ( θ/
2) has its origin in two different effects. First, the cross product term J s × m / (cid:126) in(10) describes the torque that is exerted from the magnetization on the spin polarization of theitinerant electrons. This torque is zero for parallel and antiparallel alignment of the magneticlayers and reduces the angle of the polarization of the itinerant electrons to the magnetizationfor any other alignment. Hence, for large J this contribution leads to a lowered resistance forall alignments other than parallel and antiparallel, which results in the narrow peak observedin the simulation.The second effect is a bit more subtle. For vanishing β (cid:48) the potential u from (9) varieslinearly. Small values of β (cid:48) lead to small perturbations of this linear solution. While theseperturbations have a clear effect on the overall potential difference, their effect on the spinaccumulation s due to (10) is negligible, leading to a clean sinusoidal response of the system as7hown in Fig. 4(a). With increasing β (cid:48) the perturbations of u gain influence on the solution of s which results in a distorsion of the sinusoidal response as seen in Fig. 3 and 4. We propose a three-dimensional spin-diffusion model that simultaneously solves for the spinaccumulation 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 givencurrent inflow. In order to validate the model and its implementation, we simulate the standardproblem µ MAG group and compare the outcome to results obtained withsimplified models.In a second numerical experiment, we compute the resistivity of a magnetic multilayerstructure in dependence on the tilting angle of the magnetization in the two magnetic layers.In the limit of small polarization parameter β (cid:48) and a small exchange strength J , we show thatthe resistivity is well approximated by a sine. For realistic choices of β (cid:48) and J the angulardependence shows a significantly narrower peak than the simple sine approximation. Whileexisting models already predict the observed behaviour in a macro spin approach, the presentedmodel is able to accurately describe GMR effects for both dynamically and spatially varyingmagnetization configurations. A Discretization
We solve the coupled equations (1) and (9) – (10) numerically by employing the finite-elementmethod for spatial discretization and a preconditioned BDF scheme as described in [23] forthe time integration. The demagnetization field is solved by a hybrid FEM–BEM methodas introduced in [24]. For the discretization of (9) and (10) special care has to be taken inthe 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 andintegration by parts is applied to avoid second derivatives. Furthermore, first derivatives of themagnetization m are eliminated in the same way and the integration domain for terms includingthe magnetization is restricted to the magnetized region ω in order account for discontinuities.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 tetrahedralmesh. The material parameters β , β (cid:48) , C , D , τ sf , and J are discretized with piecewise constantfunctions. For a given magnetization m , the weak version of (9) reads (cid:90) Ω C ∇ u · ∇ v d x + (cid:90) ω β (cid:48) D eµ B ( ∇ s ) T m · ∇ v d x = − (cid:90) Γ N gv d s , (20)where the current in-/outflow is given by g as Neumann condition. The additional Dirichletboundary conditions on Γ D are embedded in the function space of the solution u as usual whenemploying the finite-element method. The weak version of (10) reads − (cid:90) ω βC µ B e m ⊗ ∇ u : ∇ ζ d x + (cid:90) ∂ω ∩ Γ D βC µ B e ( ∇ u · n )( m · ζ ) d s − (cid:90) Ω D ∇ s : ∇ ζ d x − (cid:90) Ω s · ζ τ sf d x − (cid:90) ω J ( s × m ) · ζ (cid:126) d x = (cid:90) ∂ω ∩ Γ N β µ B e g ( m · ζ ) d s . (21)8 cknowledgements The financial support by the Austrian Federal Ministry of Science, Research and Economy andthe National Foundation for Research, Technology and Development as well as the AustrianScience Fund (FWF) under grant W1245 and F4112 SFB ViCoM, the Vienna Science andTechnology Fund (WWTF) under grant MA14-44, the innovative projects initiative of TUWien is gratefully acknowledged. A.M. acknowledges financial support from the King AbdullahUniversity of Science and Technology (KAUST).
References [1] Daughton, J. GMR applications.
J. Magn. Magn. Mater. , 334–342,DOI:10.1016/S0304-8853(98)00376-X (1999).[2] Freitas, P., Ferreira, R., Cardoso, S. & Cardoso, F. Magnetoresistive sensors.
J. Phys.:Condens. Matter , 165221, DOI:10.1088/0953-8984/19/16/165221 (2007).[3] Huai, Y. Spin-transfer torque MRAM (STT-MRAM): Challenges and prospects. AAPPSBulletin , 33–40 (2008).[4] Kiselev, S. I. et al. Microwave oscillations of a nanomagnet driven by a spin-polarizedcurrent.
Nature , 380–383, DOI:10.1038/nature01967 (2003).[5] Mistral, Q. et al.
Current-driven microwave oscillations in current perpendicular-to-planespin-valve nanopillars.
Applied physics letters , 192507, DOI:10.1063/1.2201897 (2006).[6] Slonczewski, J. C. Currents and torques in metallic magnetic multilayers. J. Magn. Magn.Mater. , 324–338, DOI:10.1016/S0304-8853(02)00291-3 (2002).[7] Zhang, S. & Li, Z. Roles of nonequilibrium conduction electrons on the magnetization dy-namics of ferromagnets.
Phys. Rev. Lett. , 127204, DOI:10.1103/PhysRevLett.93.127204 (2004).[8] Abert, C. et al. A three-dimensional spin-diffusion model for micromagnetics.
Sci. Rep. ,DOI:10.1038/srep14855 (2015).[9] Zhang, S., Levy, P. & Fert, A. Mechanisms of spin-polarized current-driven magnetizationswitching. Phys. Rev. Lett. , 236601, DOI:10.1103/PhysRevLett.88.236601 (2002).[10] Sturma, M., Toussaint, J.-C. & Gusakova, D. Geometry effects on magnetization dynamicsin circular cross-section wires. J. Appl. Phys. , 243901, DOI:10.1063/1.4922868 (2015).[11] Imamura, H. & Sato, J. Spin accumulation and mistracking effects on the magne-toresistance of a ferromagnetic nano-contact. vol. 266, (IOP Publishing, 2011).[12] Garc´ıa-Cervera, C. J. & Wang, X.-P. Spin-polarized currents in ferromagnetic multilayers.
Journal of Computational Physics , 699–711, DOI:10.1016/j.jcp.2006.10.029 (2007).[13] Ruggeri, M., Abert, C., Hrkac, G., Suess, D. & Praetorius, D. Coupling of dynamicalmicromagnetism and a stationary spin drift-diffusion equation: A step towards a fully self-consistent spintronics framework. Physica B , DOI:10.1016/j.physb.2015.09.003 (2015).914] Hrkac, G. et al.
Influence of eddy current on magnetization processes insubmicrometer permalloy structures.
IEEE Trans. Magn. , 3097–3099,DOI:10.1109/TMAG.2005.855234 (2005).[15] Rado, G. & Weertman, J. Spin-wave resonance in a ferromagnetic metal. J. Phys. Chem.Solids , 315–333, DOI:10.1016/0022-3697(59)90233-1 (1959).[16] Abert, C., Exl, L., Bruckner, F., Drews, A. & Suess, D. magnum. fe: A micromagneticfinite-element simulation code based on FEniCS. J. Magn. Magn. Mater. , 29–35,DOI:10.1016/j.jmmm.2013.05.051 (2013).[17] muMAG standard problem .Accessed: 2016-03-04.[18] Dieny, B. et al.
Giant magnetoresistive in soft ferromagnetic multilayers.
Phys. Rev. B ,1297, DOI:10.1103/PhysRevB.43.1297 (1991).[19] Dauguet, P. et al. Angular dependence of the perpendicular giant magnetoresistance ofmultilayers.
Phys. Rev. B , 1083, DOI:10.1103/PhysRevB.54.1083 (1996).[20] Stiles, M. D. & Zangwill, A. Noncollinear spin transfer in co/cu/co multilayers (invited). J. Appl. Phys. , 6812–6817, DOI:10.1063/1.1446123 (2002).[21] Valet, T. & Fert, A. Theory of the perpendicular magnetoresistance in magnetic multilay-ers. Phys. Rev. B , 7099–7113, DOI:10.1103/PhysRevB.48.7099 (1993).[22] Strelkov, N. et al. Spin-current vortices in current-perpendicular-to-plane nanoconstrictedspin valves.
Phys. Rev. B , 024416, DOI:10.1103/PhysRevB.84.024416 (2011).[23] Suess, D. et al. Time resolved micromagnetics using a preconditioned time integrationmethod.
J. Magn. Magn. Mater. , 298–311, DOI:10.1016/S0304-8853(02)00341-4 (2002).[24] Fredkin, D. & Koehler, T. Hybrid method for computing demagnetizing fields.
IEEETrans. Magn. , 415–417, DOI:10.1109/20.106342, 415–417, DOI:10.1109/20.106342