Finite element modeling of spin – orbit torques

Using a spin drift – diffusion model for coupled spin and charge transport, including the spin Hall effect and the inverse spin Hall effect, spin – orbit torques acting on the magnetization in heavy metal/ferromagnet bilayers are evaluated employing a finite element approach. The behavior of the resulting spin accumulation and torques is analyzed for different magnetization orientations of the ferromagnetic layer and varying thicknesses of the heavy metal layer. A strong damping-like torque component, consistent with experimental results, is observed. In addition, sub-ns switching simulations based on solving the Landau – Lifshitz – Gilbert equation are demonstrated.


Introduction
Spin-orbit torques (SOTs) allow for efficient manipulation of magnetization in emerging spintronic devices such as nano-oscillators and nonvolatile magnetoresistive random access memories (MRAM). To facilitate the research and development of such devices, it is important to have software capable of simulating the spin and magnetization dynamics involved in their operation. To be able to simulate realistic devices the flexibility of the software with regards to geometry is critical. To address these needs, a solver has been developed with the capabilities to model spin-transfer torques (STT) in STT-MRAM cells. The solver uses the finite element method to model magnetization and spin dynamics on a mesh by solving the Landau-Lifshitz-Gilbert equation (LLG) and a drift-diffusion model for spin accumulation [1], respectively. The capabilities of this solver have been extended to include the spin Hall effect (SHE), the main driving force behind SOT. In this work we present simulation results for the spin accumulation and spin torques in a heavy metal/ferromagnet (HM/FM) bilayer.

Theory & method
The spin torques acting on the magnetization of the ferromagnetic layer are given by [2] where S is the spin accumulation, m is the normalized magnetization of the ferromagnetic layer, D e is the electron diffusion coefficient, and λ J and λ φ are the exchange and spin dephasing lengths, respectively. Because the spin dynamics occur at a much shorter time scale than magnetization dynamics, we can obtain the spin accumulation by solving its equation of motion for a steady state, which is given by [1][2][3] ∂S ∂t λ sf is the spin-flip length and (∇⋅J S ) i = ∂ j (J S ) ij , where ∂ j is the partial derivative with respect to coordinate j and (J S ) i,j is the rank-2 spin current tensor describing the flow of spin component i in direction j. The spin current is given by [2,4] where the last term describes the spin current generated by the SHE.
Here μ B is the Bohr magneton, e is the elementary charge, σ is the conductivity, and β σ and β D are the charge current spin polarisation and diffusion spin polarisation, respectively. In the second and third term (∇S) ij = ∂ j S i . The ε indicates the Levi-Civita tensor and ⊗ the outer product. The strength of the SHE is described by the dimensionless spin Hall angle θ SHA . E = − ∇V is the electric field and V is the electrostatic potential obtained from solving With the assumption that the spin and charge currents vanish at boundaries not containing electrodes, one obtains the two following Neumann boundary conditions [2] ∇V⋅n = θ SHA D e σ e μ B (∇ × S)⋅n (5) and where n is the surface normal. The condition imposed by Eq. (5) describes the inverse SHE, through which a spin current can induce a charge current. In this work the spin and charge currents are kept continuous at the HM/FM boundary. Therefore, the rapid absorption of transverse spin components due to interfacial effects is modeled through the use of effective exchange and spin dephasing lengths. By deriving the weak formulation of Eqs. (2) and (4) and applying the appropriate Neumann boundary conditions given by Eqs. (5) and (6), the spin accumulation and spin torques can be calculated on a mesh with the desired geometry using the finite element method [3]. The weak formulation for the spin accumulation takes the same form as the one presented in the appendix of [1] with the exception of an additional SHE term in the linear form. Here Ω denotes a domain with non-zero θ SHA , v is the vector test function and A : B = ∑ i,j A ij B ij is the Frobenius inner product between two matrices. Furthermore, the weak formulation for the electrical potential in [1] is also extended by adding to the linear form the term corresponding to the inverse SHE ∫ where v is the scalar test function. The resulting weak formulations for the electric potential and spin accumulation are coupled since they rely on each other's solutions. Therefore, both have to be solved consecutively, until the solution for the spin accumulation converges. With the solution for S, the torques can be calculated with Eq. (1). To simulate magnetization dynamics with the LLG equation, the spin torques have to be re-calculated at every time step.
To study the dependence of the torque on several parameters it is useful to study the average torque acting on the magnetic layer, given by where V ω is the volume of the magnetic domain ω.

Results
With the approach outlined in the previous section, SOTs were simulated for a HM/FM bilayer structure using a spin Hall angle of 0.06 in the HM layer. The geometry used in the simulation is depicted in Fig. 1, where the HM layer is shown in blue while the FM layer is shown in red. The FM has a diameter of 30 nm, and the thicknesses of the HM and FM layers are 4 nm and 2 nm, respectively. The magnetization is fixed along −ŷ. A 2 mA current passes through the heavy metal layer along the x-direction, which generates an out-of-plane spin current through the SHE, resulting in a buildup of spin accumulation along the edges of the HM polarized perpendicular to both the current and the surface normals. In Fig. 2 the calculated spin accumulation is shown, with the color bar showing the magnitude and direction of the ycomponent. At the far and near sides of the FM a net spin accumulation parallel and anti-parallel to the magnetization, respectively, can be observed. This difference in accumulation on either side of the FM induces an in-plane spin current across the HM/FM interface. This effect is also present without the SHE and is caused by the spin injection of the current partly flowing through the FM layer [5]. This additional spin accumulation is polarized longitudinally to the magnetization and therefore does not contribute to the torques.
The transverse components of the spin current incident on the FM layer are absorbed, resulting in SOTs acting on the magnetization. The resulting torques are shown in Fig. 3, where the magnitudes of the torques are given by the color bar and depicted by the length of the vectors. Since the accumulation generated from the SHE is polarized parallel to the magnetization along the center of the FM layer, it does not generate any torques. Towards the edges, the polarization is not completely aligned with the magnetization, the transverse components are therefore absorbed, resulting in torques. However, due to the symmetry the  average damping-like torque is zero.
To investigate the dependence of the torque on the orientation of the magnetization, the average torque given by Eq. (9) acting on the FM was calculated for several magnetization directions along paths in the xyand yz-plane. Fig. 4 shows the average torque as a function of the angle θ between the y-axis and m, with m lying in the xy-plane. Fig. 5 shows the average torque as a function of the angle ϕ between ẑ and m with m constrained to the yz-plane. As expected the magnitude of the average torque is at its peak, where the magnetization is perpendicular to the polarization of the incident spin current, and it is lowest where magnetization and polarization are parallel. These two figures demonstrate, why it is possible to switch the magnetization in-plane deterministically in contrast to the magnetization out-of-plane. Fig. 5 shows that for out-of-plane switching of the magnetization from z to − z the torque acts against the switching after the halfway point where the magnetization is pointing along ŷ. To complete the out-of-plane switching the mirror symmetry of the structure with respect to the xzplane must be broken [6]. However, in Fig. 5, which demonstrates the torques during switching from y to − y, the torque is always working towards switching. Fig. 6 shows the dependence of the damping-like and field-like components of the torque on the thickness of the HM layer. The magnetization was kept in the z-direction and the thickness of the FM layer was kept at 10 nm. It is observed that both the damping-like and field-like components increase with increasing thickness and appear to saturate at a finite thickness of the HM layer, in agreement with [7,8]. As the torques are driven by the spin Hall effect, the damping-like component dominates. The Rashba effect at a HM/FM interface is known to give rise to a stronger field-like component. By adjusting the lengths λ J and λ φ the field-like component can be increased or decreased to reflect this behaviour.
Combining the coupled spin and charge transport calculations with a finite element LLG solver [9], enables to perform switching simulations. Fig. 7 shows the result of an in-plane switching simulation performed for the bilayer in Fig. 1, where the components of m are shown as functions of time. A voltage difference of 1 V was applied to compute the current and the spin-Hall angle was set to 0.1. The switching was achieved in    less than 0.5 ns, consistent with experimental results [10], demonstrating one of the many advantages of magnetization switching with SOTs.

Conclusion
By implementing the SHE into a spin drift-diffusion finite element solver, SOTs were modeled for a HM/FM bilayer. The behavior of the resulting torques for different orientations of the magnetization and varying thicknesses of the SHE layer is in agreement with reported experimental results [8]. Furthermore, sub-ns switching simulations were demonstrated by using our approach in conjunction with solving the Landau-Lifshitz-Gilbert equation numerically.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.