Equations

Atomistic spin model

The basic assumption of the atomistic spin model is that each lattice site is associated with a magnetic moment ${\mu }_s$. For metal systems with quenched orbital moments, the magnetic moment is mainly related to its spin angular momentum

$$\mu =-g{\mu }_B\mathbf{S}=-\hslash \gamma \mathbf{S}$$

where ${\mu }_B=e\hslash /(2m)$ is the Bohr magneton, $e(>0)$ is the electron charge, $\gamma =g{\mu }_B/\hslash (>0)$ is the gyromagnetic ratio, $g=2$ is the g-factor. The LLG equation governs the dynamics of the magnetic moment, which reads

$$\frac{\partial \mathbf{m}}{\partial t}=-\gamma \mathbf{m}\times {\mathbf{H}}_{\mathrm{eff}}+\alpha \mathbf{m}\times \frac{\partial \mathbf{m}}{\partial t}$$

where $\mathbf{m}$ is the unit vector of the magnetic moment, ${\mathbf{H}}_{\mathrm{eff}}$ is the effective field.

Unlike the continuous micromagnetic model, the atomistic spin model is discrete, and thus, the effective field ${\mathbf{H}}_{\mathrm{eff}}$ is defined as the partial derivative of the total Hamiltonian with respect to $\mathbf{m}$, i.e.,

$${\mathbf{H}}_{\mathrm{eff}}=-\frac{1}{{\mu }_s}\frac{\partial \mathcal{H}}{\partial \mathbf{m}}$$

where $\mathcal{H}$ is the total Hamiltonian including the exchange interaction, Dzyaloshinskii-Moriya interaction, dipolar interaction, anisotropy interaction, Zeeman interaction and so on.

The exchange interaction is given by

$${\mathcal{H}}_{\mathrm{ex}}=-J\sum _{⟨i,j⟩}{\mathbf{m}}_i\cdot {\mathbf{m}}_j$$

where $J$ denotes the exchange constant and $⟨i,j⟩$ represents a unique pair between lattice sites $i$ and $j$ and we assume that the summation is taken only once for each pair. A positive $J$ results in the ferromagnetic state while a negative $J$ leads to the antiferromagnetic state, which can be seen by minimizing the exchange Hamiltonian. To solve the LLG equation numerically, we need the effective field at each site, which is given

$${\mathbf{H}}_{\mathrm{ex},\mathrm{i}}=\frac{J}{{\mu }_s}\sum _{⟨i,j⟩}{\mathbf{m}}_j.$$

Similarly, the Dzyaloshinskii-Moriya interaction reads

$${\mathcal{H}}_{\mathrm{dmi}}=\sum _{⟨i,j⟩}{\mathbf{D}}_{ij}\cdot ({\mathbf{m}}_i\times {\mathbf{m}}_j)$$

where ${\mathbf{D}}_{ij}$ is the DM vector. Typically, there are two situations: (1) ${\mathbf{D}}_{ij}=D{\hat{r}}_{ij}$, which corresponds to the Bulk DMI. (2) ${\mathbf{D}}_{ij}=D{\hat{r}}_{ij}\times \hat{z}$, which corresponds to the interfacial DMI. Unlike the exchange interaction, the DM Hamiltonian is minimal when the two adjacent magnetic moments are perpendicular. The corresponding effective field can be computed by

$${\mathbf{H}}_{\mathrm{dmi},\mathrm{i}}=\frac{1}{{\mu }_s}\sum _{⟨i,j⟩}{\mathbf{D}}_{ij}\times {\mathbf{m}}_j.$$

The other two common interactions are the anisotropy and the Zeeman, which are given by

$$\begin{gathered} {\mathcal{H}}_{\mathrm{an}}=-K\sum _i{({\mathbf{m}}_i\cdot \hat{u})}^2 \\ {\mathcal{H}}_{\mathrm{ze}}=-{\mu }_s{\mathbf{m}}_i\cdot \mathbf{H} \end{gathered}$$

The effective field of the anisotropy is

$${\mathbf{H}}_{\mathrm{an},\mathrm{i}}=\frac{2K}{{\mu }_s}({\mathbf{m}}_i\cdot \hat{u})\hat{u}.$$

Multiferroic Insulators

For multiferroics such as Cu${}_2$OSeO${}_3$, the noncollinear spin texture induces a local electric polarization via the the $d-p$ hybridization mechanism. The local electric dipole moment $\mathbf{P}$ depends on the direction of the applied static magnetic field relative to the crystallographic axes.[PRL 128, 037201 (2022)]

$H_0//001$:

$${\mathbf{P}}_i=\lambda (-m_{i,z}{m}_{i,x},m_{i,y}{m}_{i,z},\frac{-m_{i,x}^2+m_{i,y}^2}{2})$$

$H_0//110$:

$${\mathbf{P}}_i=\lambda (-m_{i,x}{m}_{i,y},\frac{-m_{i,x}^2+m_{i,z}^2}{2},m_{i,y}{m}_{i,z})$$

$H_0//111$:

$$\begin{gathered} {\mathbf{P}}_{i,x}=-\lambda \frac{m_{i,x}(\sqrt{2}{m}_{i,y}+m_{i,z})}{\sqrt{3}} \\ {\mathbf{P}}_{i,y}=\lambda \frac{-m_{i,x}^2+m_{i,y}(m_{i,y}-\sqrt{2}{m}_{i,z})}{\sqrt{6}} \\ {\mathbf{P}}_{i,z}=-\lambda \frac{m_{i,x}^2+m_{i,y}^2-2m_{i,z}^2}{2\sqrt{3}} \end{gathered}$$

The Hamiltonian related to the high-frequency lasers is given by

$${\mathcal{H}}_{\mathrm{laser}}=-\sum _i{\mu }_s{\mathbf{m}}_i\cdot \mathbf{H}(t)-\sum _i{\mathbf{P}}_i\cdot \mathbf{E}(t)$$

The corresponding effective fields associated with the electric field are

$H_0//001$:

$$\begin{gathered} H_x=\frac{\lambda }{{\mu }_s}(-E_z{m}_x-E_x{m}_z) \\ {H}_y=\frac{\lambda }{{\mu }_s}(E_z{m}_y+E_y{m}_z) \\ {H}_z=\frac{\lambda }{{\mu }_s}(-E_x{m}_x+E_y{m}_y) \end{gathered}$$

$H_0//110$:

$$\begin{gathered} H_x=\frac{\lambda }{{\mu }_s}(-E_y{m}_x-E_x{m}_y) \\ {H}_y=\frac{\lambda }{{\mu }_s}(-E_x{m}_x+E_z{m}_z) \\ {H}_z=\frac{\lambda }{{\mu }_s}(E_z{m}_y+E_y{m}_z) \end{gathered}$$

$H_0//111$:

$$\begin{gathered} H_x=-\frac{\lambda }{\sqrt{3}{\mu }_s}(\sqrt{2}(E_y{m}_x+E_x{m}_y)+E_z{m}_x+E_x{m}_z) \\ {H}_y=-\frac{\lambda }{\sqrt{3}{\mu }_s}(\sqrt{2}(E_x{m}_x-E_y{m}_y)+E_z{m}_y+E_y{m}_z) \\ {H}_z=-\frac{\lambda }{\sqrt{3}{\mu }_s}(E_x{m}_x+E_y{m}_y-2E_z{m}_z) \end{gathered}$$

Micromagnetic model

In micromagnetics, the effective field can be computed from the total micromagnetic energy

$${\mathbf{H}}_{\mathrm{eff}}=-\frac{1}{{\mu }_0{M}_s}\frac{\delta E}{\delta \mathbf{m}}$$

The typical energy terms are

• Exchange energy

$$E_{\mathrm{ex}}={\int }_{V}A(\mathrm{\nabla }\mathbf{m}{)}^2\mathrm{d}V$$

where $(\mathrm{\nabla }\mathbf{m}{)}^2={\left(\mathrm{\nabla }{m}_x\right)}^2+{\left(\mathrm{\nabla }{m}_y\right)}^2+{\left(\mathrm{\nabla }{m}_z\right)}^2$. So the corresponding effective field is

$${\mathbf{H}}_{\mathrm{ex}}=\frac{2A}{{\mu }_0{M}_s}{\mathrm{\nabla }}^2\mathbf{m}$$

• Zeeman energy

$$E_{\mathrm{ex}}=-{\mu }_0{\int }_V\mathbf{H}\cdot \mathbf{M}\mathrm{d}V$$

as expected, the effective field is $\mathbf{H}$.

• Anisotropy The uniaxial anisotropy energy is given by

$$E_{\mathrm{anis}}=-{\int }_V{K}_u(\mathbf{m}\cdot \hat{u}{)}^{2}dV$$

from which the effective field can be computed as

$${\mathbf{H}}_{\mathrm{an}}=\frac{2K_u}{{\mu }_0{M}_s}(\mathbf{m}\cdot \hat{u})\hat{u}$$

• Cubic Anisotropy The cubic anisotropy energy is given by

$$E_{\mathrm{cubic}}=-{\int }_V{K}_c(m_x^4+m_y^4+m_z^4)dV$$

and thus the corresponding effective field reads

$${\mathbf{H}}_{\mathrm{cubic}}=\frac{4K_c}{{\mu }_0{M}_s}(m_x^3{\mathbf{e}}_x+m_y^3{\mathbf{e}}_y+m_z^3{\mathbf{e}}_z)$$

• Hexagonal Anisotropy The energy density of the hexagonal anisotropy is given by:

$$E=K_1{\sin }^2\theta +K_2{\sin }^4\theta +K_3{\sin }^6\theta \cos 6\varphi$$

Here, $\theta$ is the angle between the magnetization vector $\mathbf{m}$ and the $c$-axis ($z$-axis). $\varphi$ is the angle of the projection of $\mathbf{m}$ on the hexagonal plane, measured with respect to the $x$-axis.

Using the identity $\cos 6x=-{\sin }^{6}x+15{\cos }^{2}x{\sin }^{4}x-15{\cos }^{4}x{\sin }^{2}x+{\cos }^{6}x$, the energy density can be rewritten in terms of $m_x$, $m_y$, and $m_z$ as:

$$E=K_1(1-m_z^2)+K_2(1-m_z^2{)}^2+K_3(m_x^6-15m_x^4{m}_y^2+15m_x^2{m}_y^4-m_y^6)$$

Therefore, The corresponding effective field is given by:

$${\mathbf{H}}_{\mathrm{eff}}=-\frac{6K_3}{{\mu }_0{M}_s}(m_x^5-10m_x^3{m}_y^2+5m_x{m}_y^4){\mathbf{e}}_x-\frac{6K_3}{{\mu }_0{M}_s}(-5m_x^4{m}_y+10m_x^2{m}_y^3-m_y^5){\mathbf{e}}_y+\frac{2m_z}{{\mu }_0{M}_s}[K_1+2K_2(1-m_z^2)]{\mathbf{e}}_z$$

Dzyaloshinskii-Moriya Energy

In the continuum limit, the DMI energy density $w_{\mathrm{dmi}}$ is associated with the so-called Lifshitz invariants, which are terms in the form

$$L_{ij}^{(k)}=m_i\frac{\partial m_j}{\partial x_k}-m_j\frac{\partial m_i}{\partial x_k}.$$

The form of DMI energy density varies depending on the symmetry class. For bulk DMI, corresponding to symmetry class $T$ or $O$, the expression is given by:

$$w_{\mathrm{dmi}}=D(L_{yx}^{(z)}+L_{xz}^{(y)}+L_{zy}^{(x)})=D\mathbf{m}\cdot (\mathrm{\nabla }\times \mathbf{m}).$$

The associated effective field is

$${\mathbf{H}}_{\mathrm{dmi}}=-\frac{2D}{{\mu }_0{M}_s}(\mathrm{\nabla }\times \mathbf{m}).$$

For a thin film with interfacial DMI or a crystal with symmetry class $C_{nv}$, the energy density is

$$w_{\mathrm{dmi}}=D(L_{xz}^{(x)}+L_{yz}^{(y)})=D(\mathbf{m}\cdot \mathbf{\nabla }{m}_z-m_z\mathbf{\nabla }\cdot \mathbf{m}),$$

and the effective field is

$${\mathbf{H}}_{\mathrm{dmi}}=-\frac{2D}{{\mu }_0{M}_s}({\mathbf{e}}_y\times \frac{\partial \mathbf{m}}{\partial x}-{\mathbf{e}}_x\times \frac{\partial \mathbf{m}}{\partial y}).$$

For a crystal with symmetry class $D_{2d}$, the DMI energy density is given by $w_{\mathrm{dmi}}=D(L_{xz}^{(y)}+L_{yz}^{(x)})$, resulting in the effective field

$${\mathbf{H}}_{\mathrm{dmi}}=-\frac{2D}{{\mu }_0{M}_s}({\mathbf{e}}_y\times \frac{\partial \mathbf{m}}{\partial y}-{\mathbf{e}}_x\times \frac{\partial \mathbf{m}}{\partial x}).$$

Although the effective fields for different symmetries differ, the numerical implementation can be unified as follows

$${\mathbf{H}}_{\mathrm{dmi},\mathrm{i}}=-\frac{1}{{\mu }_0{M}_s}\sum _{j\in N_i}{D}_{ij}\frac{{\mathbf{e}}_{ij}\times {\mathbf{m}}_j}{{\mathrm{\Delta }}_{ij}},$$

where $D_{ij}$ represents the effective DMI constant and ${\mathbf{e}}_{ij}$ denotes the DMI vectors. For bulk DMI, ${\mathbf{e}}_{ij}={\hat{\mathbf{r}}}_{ij}$ where ${\hat{\mathbf{r}}}_{ij}$ is the unit vector between cell $i$ and cell $j$. For interfacial DMI,  ${\mathbf{e}}_{ij}={\mathbf{e}}_z\times {\hat{\mathbf{r}}}_{ij}$, i.e., ${\mathbf{e}}_{ij}=\{-{\mathbf{e}}_y,{\mathbf{e}}_y,{\mathbf{e}}_x,-{\mathbf{e}}_x,0,0\}$ for the 6 neighbors $N_i=\{-x,+x,-y,+y,-z,+z\}$. For the symmetry class $D_{2d}$ one has ${\mathbf{e}}_{ij}=\{{\mathbf{e}}_x,-{\mathbf{e}}_x,-{\mathbf{e}}_y,{\mathbf{e}}_y,0,0\}$. If the cell-based DMI is provided, the effective DMI constant can be computed as

$$D_{ij}=\frac{2D_i{D}_j}{D_i+D_j}.$$

• Bulk DMI energy The Bulk DMI energy reads

$$E_{\mathrm{dmi}}={\int }_{V}D\mathbf{m}\cdot (\mathrm{\nabla }\times \mathbf{m})\mathrm{d}V$$

so the effective field is

$${\mathbf{H}}_{\mathrm{D}}=-\frac{2D}{{\mu }_0{M}_s}(\mathrm{\nabla }\times \mathbf{m})$$

• Magnetostatic energy

$$E_{\mathrm{d}}=-\frac{{\mu }_0}{2}{\int }_V{\mathbf{H}}_{\mathrm{d}}(\mathbf{r})\cdot \mathbf{M}(\mathbf{r})dV$$$${\mathbf{H}}_{\mathrm{d}}(\mathbf{r})=\frac{1}{4\pi }({\int }_V{\rho }_m\left({\mathbf{r}}^{\prime }\right)\frac{\mathbf{r}-{\mathbf{r}}^{\prime }}{{|\mathbf{r}-{\mathbf{r}}^{\prime }|}^3}{\mathrm{d}}^3{r}^{\prime }+{\int }_S{\sigma }_m\left({\mathbf{r}}^{\prime }\right)\frac{\mathbf{r}-{\mathbf{r}}^{\prime }}{{|\mathbf{r}-{\mathbf{r}}^{\prime }|}^3}{\mathrm{d}}^2{r}^{\prime })$$

LLG equation

The driver LLG solves the standard LLG equation, which can be written as

$$\frac{\partial \mathbf{m}}{\partial t}=-\gamma \mathbf{m}\times \mathbf{H}+\alpha \mathbf{m}\times \frac{\partial \mathbf{m}}{\partial t}$$

and the corresponding LL form is given by

$$(1+{\alpha }^2)\frac{\partial \mathbf{m}}{\partial t}=-\gamma \mathbf{m}\times \mathbf{H}-\alpha \gamma \mathbf{m}\times (\mathbf{m}\times \mathbf{H})$$

Inertial LLG Equation

The Inertial LLG equation describes the dynamics of magnetization with inertial effects, extending the classical LLG equation by including a second-order time derivative term:

$$\frac{d\mathbf{m}}{dt}=-\gamma \mathbf{m}\times {\mathbf{H}}_{\mathrm{eff}}+\alpha \mathbf{m}\times \frac{d\mathbf{m}}{dt}+\eta \mathbf{m}\times \frac{d^2\mathbf{m}}{d{t}^2}$$

where $\mathbf{m}$ is the normalized magnetization vector ($|\mathbf{m}|=1$), $\gamma$ is the gyromagnetic ratio, ${\mathbf{H}}_{\mathrm{eff}}$ is the effective magnetic field, $\alpha$ is the Gilbert damping parameter, $\eta =\alpha \tau$ is the inertial parameter, with $\tau$ being the angular momentum relaxation time.

To facilitate numerical implementation, we introduce the velocity variable:

$$\mathbf{v}=\frac{d\mathbf{m}}{dt}$$

Substituting this definition, the ILLG equation becomes:

$$\mathbf{v}=-\gamma \mathbf{m}\times {\mathbf{H}}_{\mathrm{eff}}+\alpha \mathbf{m}\times \mathbf{v}+\eta \mathbf{m}\times \frac{d\mathbf{v}}{dt}$$

Taking the cross product of both sides with $\mathbf{m}$:

$$\mathbf{m}\times \mathbf{v}=-\gamma \mathbf{m}\times (\mathbf{m}\times {\mathbf{H}}_{\mathrm{eff}})+\alpha \mathbf{m}\times (\mathbf{m}\times \mathbf{v})+\eta \mathbf{m}\times (\mathbf{m}\times \frac{d\mathbf{v}}{dt})$$

Applying the vector triple product identity $\mathbf{a}\times (\mathbf{b}\times \mathbf{c})=(\mathbf{a}\cdot \mathbf{c})\mathbf{b}-(\mathbf{a}\cdot \mathbf{b})\mathbf{c}$ to the inertial term:

$$\eta \mathbf{m}\times (\mathbf{m}\times \frac{d\mathbf{v}}{dt})=\eta (\mathbf{m}\cdot \frac{d\mathbf{v}}{dt})\mathbf{m}-\eta \frac{d\mathbf{v}}{dt}=-\eta v^2\mathbf{m}-\eta \frac{d\mathbf{v}}{dt}$$

The simplification $(\mathbf{m}\cdot \frac{d\mathbf{v}}{dt})=-v^2$ follows from differentiating the constraint $\mathbf{m}\cdot \mathbf{v}=0$, which itself derives from the conservation of magnetization magnitude $|\mathbf{m}|=1$. The ILLG equation can be transformed into the following coupled first-order system suitable for numerical integration:

$$\begin{array}{rl}\frac{d\mathbf{m}}{dt}& =\mathbf{v}\\ \frac{d\mathbf{v}}{dt}& =-\frac{1}{\eta }\mathbf{m}\times \mathbf{v}-\frac{\gamma }{\eta }\mathbf{m}\times (\mathbf{m}\times {\mathbf{H}}_{\mathrm{eff}})+\frac{\alpha }{\eta }\mathbf{m}\times (\mathbf{m}\times \mathbf{v})-v^2\mathbf{m}\end{array}$$

Spin Transfer Torques

Spin transfer torque (STT) is a fundamental phenomenon in spintronics where spin-polarized electric currents exert torques on magnetic moments, enabling the manipulation of magnetization states. This effect is typically modeled through extensions to the Landau-Lifshitz-Gilbert (LLG) equation.

Zhang-Li Model

The Zhang-Li model incorporates STT into the LLG equation. The equation for the magnetization vector m is:

$$\frac{d\mathbf{m}}{dt}=-\gamma \mathbf{m}\times {\mathbf{H}}_{\mathrm{eff}}+\alpha \mathbf{m}\times \frac{d\mathbf{m}}{dt}-b\mathbf{m}\times [\mathbf{m}\times (\mathbf{j}\cdot \mathrm{\nabla })\mathbf{m}]-\xi b\mathbf{m}\times (\mathbf{j}\cdot \mathrm{\nabla })\mathbf{m}$$

where j is the current density vector and $\xi$ is the non-adiabatic parameter. The coefficient $b$ is defined as:

$$b=\frac{P{\mu }_B}{e{M}_s(1+{\xi }^2)}$$

with $P$ being the spin polarization, ${\mu }_B$ the Bohr magneton, and $e>0$ the elementary charge.

Note that $\mathbf{m}\times [\mathbf{m}\times (\mathbf{j}\cdot \mathrm{\nabla })\mathbf{m}]=-(\mathbf{j}\cdot \mathrm{\nabla })\mathbf{m}$, the equation simplifies to:

$$\frac{\partial \mathbf{m}}{\partial t}=-\gamma \mathbf{m}\times {\mathbf{H}}_{\mathrm{eff}}+\alpha \mathbf{m}\times \frac{\partial \mathbf{m}}{\partial t}-(\mathbf{u}\cdot \mathrm{\nabla })\mathbf{m}+\beta [\mathbf{m}\times (\mathbf{u}\cdot \mathrm{\nabla })\mathbf{m}]$$

where

$$\mathbf{u}=-b\mathbf{J},\beta =\xi$$

and the unit of u is m/s.

Alternative Formulations

Other models may define the current strength parameter differently. For comparison, $\mathbf{u}$ is defined as follows:

$$\mathbf{u}=-\frac{Pg{\mu }_B}{2e{M}_s}\mathbf{j}$$

For detailed discussions on model differences, see this reference.

Atomistic Model

In atomistic simulations, the current strength is adapted for discrete spins:

$$\mathbf{u}=-\frac{Pg{\mu }_B{a}^3}{2e{\mu }_s}\mathbf{j}=-\frac{P{a}^3}{2eS}\mathbf{j}$$

where $a$ is the lattice constant, ${\mu }_s$ is the magnetic moment, $S=|\mathbf{S}|$ is the magnitude of the local spin vector.

Slonczewski spin transfer torque

The LLG equation including the Slonczewski torque is

$$\frac{d\mathbf{m}}{dt}=-\gamma \mathbf{m}\times {\mathbf{H}}_{\text{eff}}+\alpha (\mathbf{m}\times \frac{d\mathbf{m}}{dt})+\gamma \beta ϵ[\mathbf{m}\times ({\mathbf{m}}_p\times \mathbf{m})]-\gamma \beta \xi \mathbf{m}\times {\mathbf{m}}_p$$

where $\mathbf{m}=\mathbf{M}/M_s$, $\gamma$ the gyromagnetic ratio, ${\mathbf{m}}_p$ is electron polarization direction, and $\xi$ is the secondary spin-torque parameter. The coefficient $\beta$ is defined as follows:

$$\beta =\frac{\hslash }{{\mu }_{0}e}\frac{J}{t{M}_s}$$

where $e$ is electron charge (C), $J$ is current density (A/m²), $t$ is free layer thickness (m), $M_s$ is the saturation magnetization (A/m). The coefficient $ϵ$ is defined as

$$ϵ=\frac{P{\mathrm{\Lambda }}^2}{({\mathrm{\Lambda }}^2+1)+({\mathrm{\Lambda }}^2-1)(\mathbf{m}\cdot {\mathbf{m}}_p)}$$

where $P$ is the spin polarization, $\mathrm{\Lambda }$ is the Slonczewski parameter.

For constant $ϵ=P/2$ which corresponding to $\mathrm{\Lambda }=1$, the Slonczewski torque reduced to [PRL 102, 067206 (2009)]

$$\frac{\partial \mathbf{m}}{\partial t}=-\gamma \mathbf{m}\times {\mathbf{H}}_{\text{eff}}+\alpha \mathbf{m}\times \frac{\partial \mathbf{m}}{\partial t}-a_J\mathbf{m}\times (\mathbf{m}\times {\mathbf{m}}_p)-b_J\mathbf{m}\times {\mathbf{m}}_p$$

where

$$a_J=\frac{\hslash \gamma }{2{\mu }_{0}e}\frac{PJ}{t{M}_s},b_J=\xi a_J$$

Damping-Like and Field-Like Torques

The Landau-Lifshitz-Gilbert (LLG) equation can be extended to include damping-like and field-like torques, which are essential for modeling spin-orbit torques (SOT) or Slonczewski torque with constant $ϵ(\theta )$. The modified LLG equation is expressed as:

$$\frac{\partial \mathbf{m}}{\partial t}=-\gamma \mathbf{m}\times {\mathbf{H}}_{\text{eff}}+\alpha \mathbf{m}\times \frac{\partial \mathbf{m}}{\partial t}-a_J\mathbf{m}\times (\mathbf{m}\times \mathbf{p})-b_J\mathbf{m}\times \mathbf{p}$$

where $a_J$ and $b_J$ are the coefficients of the damping-like and field-like torques, respectively, and $\mathbf{p}$ is a unit vector. In the context of SOT in particular, $\mathbf{p}$ represents the direction of the spin current generated by spin-orbit coupling.

Effective Field Formulation

For numerical simulations, spin-transfer torque (STT) can be equivalently expressed as an effective field contribution. The LLG equation with STT represented as an effective field becomes:

$$\frac{\partial \mathbf{m}}{\partial t}=-\gamma \mathbf{m}\times ({\mathbf{H}}_{\text{eff}}+{\mathbf{H}}_{\text{stt}})+\alpha \mathbf{m}\times \frac{\partial \mathbf{m}}{\partial t}$$

In the Zhang-Li model for STT, the effective field is given:

$${\mathbf{H}}_{\text{stt}}=-\frac{1}{\gamma }[\mathbf{m}\times (\mathbf{u}\cdot \mathrm{\nabla })\mathbf{m}+\beta (\mathbf{u}\cdot \mathrm{\nabla })\mathbf{m}]$$

For damping-like and field-like torques, the effective field can be derived as:

$${\mathbf{H}}_{\text{stt}}=\frac{1}{\gamma }(a_J\mathbf{m}\times \mathbf{p}+b_J\mathbf{p})$$

SLLG equation

The SLLG equation, i.e., LLG equation including the stochastic field $\mathbf{b}$, is given by

$$\frac{\partial \mathbf{m}}{\partial t}=-\gamma \mathbf{m}\times ({\mathbf{H}}_{\mathrm{eff}}+\mathbf{b})+\alpha \mathbf{m}\times \frac{\partial \mathbf{m}}{\partial t}$$

The thermal fluctuation is assumed to be a Gaussian white noise, i.e., the thermal noise $\mathbf{b}$ obeys the properties

$$⟨\mathbf{b}⟩=0,⟨{\mathbf{b}}_i^u\cdot {\mathbf{b}}_j^v⟩=2D{\delta }_{ij}{\delta }_{uv}$$

where $i$ and $j$ are Cartesian indices, $u$ and $v$ indicate the magnetization components and $⟨\cdot ,\cdot ⟩$ represents the average taken over different realizations of the fluctuating field. And

$$D=\frac{\alpha k_{B}T}{\gamma {\mu }_s}.$$

For the micromagnetic case, $D$ is given as

$$D=\frac{\alpha k_{B}T}{{\mu }_0{M}_s\gamma \mathrm{\Delta }V}.$$

which is equivalent to a stochastic field

$${\mathbf{b}}^u=\eta \sqrt{\frac{2\alpha k_{B}T}{{\mu }_0{M}_s\gamma \mathrm{\Delta }Vdt}}$$

where $\eta$ is a random number follows the normal distribution.

Steepest descent method

We provide a steepest descent energy minimization method for a complicated system, which is of the form

$$x_{k+1}=x_k+{\alpha }_k{d}_k$$

where

$$d_k=-\mathrm{\nabla }f(x_k)$$

And for the micromagnetics, we have

$${\mathbf{m}}_{k+1}={\mathbf{m}}_k-{\tau }_k{\mathbf{m}}_k\times ({\mathbf{m}}_k\times {\mathbf{H}}_{\mathrm{eff}})$$

In practice, we use the following update rule to keep the magnetization vector normalized.

$${\mathit{m}}_{k+1}={\mathit{m}}_k-{\tau }_k\frac{{\mathit{m}}_k+{\mathit{m}}_{k+1}}{2}\times ({\mathit{m}}_k\times {\mathit{H}}_{\mathrm{eff}}\left({\mathit{m}}_k\right))$$$${\mathit{m}}_{k+1}^2={\mathit{m}}_k^2$$

From the equation we have:

$$(1+\frac{{\tau }_k^2}{4}{\mathit{f}}_k^2){\mathbf{m}}_{k+1}=(1-\frac{{\tau }_k^2}{4}{\mathit{f}}_k^2){\mathbf{m}}_k-{\tau }_k{\mathbf{g}}_k$$

where

$$\begin{array}{rl}{\mathbf{f}}_k& ={\mathbf{m}}_k\times {\mathbf{H}}_{\mathrm{eff}}\\ {\mathit{g}}_k& ={\mathit{m}}_k\times ({\mathit{m}}_k\times {\mathit{H}}_{\mathrm{eff}})\end{array}$$

The step size ${\tau }_k$ can be computed by

$${\tau }_k^1=\frac{\sum _i{\mathit{s}}_{k-1}^i\cdot {\mathit{s}}_{k-1}^i}{\sum _i{\mathit{s}}_{k-1}^i\cdot {\mathit{y}}_{k-1}^i},{\tau }_k^2=\frac{\sum _i{\mathit{s}}_{k-1}^i\cdot {\mathit{y}}_{k-1}^i}{\sum _i{\mathit{y}}_{k-1}^i\cdot {\mathit{y}}_{k-1}^i}$$

where

$$\begin{array}{rl}{\mathit{s}}_{k-1}& ={\mathit{m}}_k-{\mathit{m}}_{k-1}\\ {\mathit{y}}_{k-1}& ={\mathit{g}}_k-{\mathit{g}}_{k-1}\end{array}$$

Monte Carlo Simulation

The implemented energy reads

$$\mathcal{H}=-\sum _{⟨i,j⟩}(J_x{S}_i^x{S}_j^x+J_y{S}_i^y{S}_j^y+J_z{S}_i^z{S}_j^z)+\sum _{⟨i,j⟩}{\mathbf{D}}_{ij}\cdot ({\mathbf{S}}_i\times {\mathbf{S}}_j)-K\sum _i{(\mathbf{u}\cdot {\mathbf{S}}_i)}^2-\sum _i\mathbf{H}\cdot {\mathbf{S}}_i$$

where ${\mathbf{S}}_i$ is unit vector of the classical spin at site i.

For interfacial DMI,

$${\mathbf{D}}_{ij}=D\hat{z}\times {\hat{r}}_{ij}+D_z^j\hat{z}$$

while for Bulk DMI,

$${\mathbf{D}}_{ij}=D{\hat{r}}_{ij}$$

Note that the Monte Carlo only works for triangular and cubic meshes.

NEB (Nudged elastic band)

NEB is a chain method to find the MEP (minimum energy path) between two states. To start, we need to construct a chain including several images (each image is a copy of the magnetization) and then relax the system. Two ends images that corresponding to the initial and final states will be pinned as they are the energy states given by the users. The system contain all free images will be relaxed to reduce the total energy, which is very similar to the case that relaxing the magnetic system using LLG equation if one disables the precession term. One significant difference is that the effective field in LLG equation is the functional derivative of the system energy with respect to magnetization while in NEB the effective field of image n should also contain the influence of its neighbours (i.e., the images n-1 and n+1). This influence is described by the so-called tangents: only the perpendicaular part of the effective field is kept when relaxing the whole system.

Assume that the whole system has N images

$$\mathbf{Z}=[{\mathbf{Y}}_1,{\mathbf{Y}}_2,...,{\mathbf{Y}}_N]$$

where

$${\mathbf{Y}}_i=[m_{1x},m_{1y},m_{1z},...,m_{nx},m_{ny},m_{nz}]$$

each image has n spins. To relax the system, we could solve the equation

$$\frac{{\mathbf{Y}}_i}{\partial t}=-{\mathbf{Y}}_i\times ({\mathbf{Y}}_i\times {\mathbf{G}}_i)$$

where ${\mathbf{G}}_i$ is effective field that can be computed as

$${\mathbf{Y}}_i={\mathbf{H}}_i-({\mathbf{H}}_i\cdot {\mathbf{t}}_i){\mathbf{t}}_i+{\mathbf{F}}_i$$

The ${\mathbf{H}}_i$ is the normal micromagnetic effective field, ${\mathbf{t}}_i$ is the tangent and ${\mathbf{F}}_i$ is a force that can be used to adjust the distance between images.

$${\mathbf{F}}_i=k(|{\mathbf{Y}}_{i+1}-{\mathbf{Y}}_i|-|{\mathbf{Y}}_i-{\mathbf{Y}}_{i-1}|){\mathbf{t}}_i$$

The distance bewteen images ${\mathbf{Y}}_i$ and ${\mathbf{Y}}_j$ is defined as

$$L={[\sum _k(L_k^{i,j}{)}^2]}^{1/2}$$

where $L_k^{i,j}$ is the geodesic distance of point k that can be computed using Vincenty's formula.

The tangents can be computed as follows

$$\begin{gathered} {\mathbf{t}}_i^{+}={\mathbf{Y}}_{i+1}-{\mathbf{Y}}_i \\ {\mathbf{t}}_i^{-}={\mathbf{Y}}_i-{\mathbf{Y}}_{i-1} \end{gathered}$$

The detailed equations can be found @ [Journal of Chemical Physics 113, 22 (2000)] and [Computer Physics Communications 196 (2015) 335–347].