The tangent-plane scheme (TPS), introduced by Alouges [135] and further developed by Abert et al. [90] and Praetorius et al. [160], reformulates the LLG equation to
exploit its geometric structure. Rather than solving directly for \(\mathbf {m}_{n+1}\), the scheme solves for the velocity \(\mathbf {v} = \partial _t \mathbf {m}\) constrained to lie in the tangent plane at the current
magnetization. This reformulation offers two key advantages: it transforms the nonlinear problem into a linear one and enforces the unit-sphere constraint to \(\mathcal {O}(\Delta t^2)\) via the tangent-plane projection, with
exact restoration achieved through post-step renormalization.
where \(\mathbf {H}_\text {eff}\) is the effective field defined in Equation (3.19), can be rewritten by taking the cross product with
\(\mathbf {m}\) and using the vector triple product identity. With \(\mathbf {v} = \partial _t \mathbf {m}\), one obtains:
The RHS is the projection of \(\mathbf {H}_\text {eff}\) onto the tangent plane at \(\mathbf {m}\). Crucially, if \(\mathbf {v}\) satisfies the constraint \(\mathbf {v} \cdot \mathbf {m} = 0\), then the RHS
automatically lies in the tangent plane, and the equation is consistent.
5.3.2 Weak Formulation
The finite-element discretization proceeds by multiplying Equation (5.22) by a test function \(\mathbf {w}\) and integrating
over the domain \(\Omega \). In weak form:
where \(\nabla \mathbf {m}: \nabla \mathbf {w} = \sum _{i,j} \partial _j m_i \partial _j w_i\) denotes the Frobenius inner product, and homogeneous Neumann boundary conditions have been assumed. Other
field contributions (external, anisotropy, demagnetization) enter as volume integrals against \(\mathbf {w}\).
5.3.3 Complete Weak Form with Spin-Transfer Torque
For simulations of current-induced magnetization dynamics, spin-transfer torque contributions must be incorporated. The spin-accumulation \(\mathbf {S}\), computed from the drift-diffusion equations, enters via the torque
expression derived in Section 6.1. The complete weak form reads:
\(\seteqnumber{0}{5.}{24}\)
\begin{align}
\int _{\Omega _m} \bigl ( \alpha \mathbf {v} + \mathbf {m}^n \times \mathbf {v} \bigr ) \cdot \mathbf {w} \, \mathrm {d}\mathbf {x} &= \frac {2A\gamma _0}{\mu _0 M_s} \int _{\Omega _m}
\nabla \mathbf {m}^n : \nabla \mathbf {w} \, \mathrm {d}\mathbf {x} \nonumber \\ &\quad - \gamma _0 \int _{\Omega _m} \bigl ( \mathbf {H}_\text {ext} + \mathbf {H}_\text {ani}(\mathbf
{m}^n) + \mathbf {H}_\text {dem}(\mathbf {m}^n) \bigr ) \cdot \mathbf {w} \, \mathrm {d}\mathbf {x} \nonumber \\ &\quad + \frac {D_e}{M_s} \int _{\Omega _m} \left ( \frac {\mathbf
{S}^n}{\lambda _J^2} + \frac {\mathbf {m}^n \times \mathbf {S}^n}{\lambda _\varphi ^2} \right ) \cdot \mathbf {w} \, \mathrm {d}\mathbf {x}, \label {eq::TI_full_weak_STT}
\end{align}
where \(D_e\) is the electron diffusion coefficient, and \(\lambda _J\), \(\lambda _\varphi \) are the exchange and dephasing lengths characterizing the spin-transfer torque. The integrals are evaluated only over the magnetic
subdomain \(\Omega _m\).
5.3.4 Time Discretization: The \(\theta \)-Method
A \(\theta \)-method time discretization treats the effective field at a weighted average of the current and next time levels:
For the exchange field, which depends linearly on \(\mathbf {m}\), this introduces implicit dependence on \(\mathbf {m}^{n+1} = \mathbf {m}^n + \Delta t \mathbf {v}\). The key observation is that:
In the full weak formulation (Equation (5.25)), the first term is known from \(\mathbf {m}^n\) and moves to the right-hand side,
while the second term depends on the unknown \(\mathbf {v}\) and contributes to the left-hand side.
5.3.5 The Saddle-Point System
Enforcement of the tangent-plane constraint \(\mathbf {v} \cdot \mathbf {m} = 0\) via Lagrange multipliers leads to a block-structured linear system. Let \(\mathbf {v}_h\) and \(\lambda _h\) denote the finite-element
approximations to the velocity and Lagrange multiplier, expanded in the basis \(\{\boldsymbol {\phi }_i\}\) defined in Section 4.2.2:
where \(\phi _i\) are scalar basis functions (with \(\boldsymbol {\phi }_i\) the corresponding vector-valued functions). The constraint is imposed weakly:
(i) \(A\) is the \((1,1)\)-block with entries \(A_{ij} = \int _\Omega (\alpha \boldsymbol {\phi }_j + \mathbf {m}^n \times \boldsymbol {\phi }_j) \cdot \boldsymbol {\phi
}_i \, d\mathbf {x} + \theta \Delta t K_{ij}\),
(ii) \(K_{ij} = \frac {2A\gamma _0}{\mu _0 M_s} \int _\Omega \nabla \boldsymbol {\phi }_j : \nabla \boldsymbol {\phi }_i \, d\mathbf {x}\) is the exchange stiffness
matrix,
(iii) \(B\) is the constraint matrix with entries \(B_{ij} = \int _\Omega (\mathbf {m}^n \cdot \boldsymbol {\phi }_j) \phi _i \, d\mathbf {x}\),
(iv) \(\mathbf {f}\) is the RHS incorporating explicit field contributions.
The zero block in the \((2,2)\) position renders the overall system indefinite, giving rise to a saddle-point problem [163]: the solution corresponds to a minimum with respect to \(\mathbf {v}\) and a maximum with
respect to \(\boldsymbol {\lambda }\).
The \((1,1)\)-block \(A\) has the structure:
\(\seteqnumber{0}{5.}{31}\)
\begin{equation}
A = \alpha M + [\mathbf {m}^n]_\times + \theta \Delta t K, \label {eq::TI_TPS_A_structure}
\end{equation}
where \(M\) is the mass matrix, \([\mathbf {m}^n]_\times \) is the skew-symmetric matrix representing the cross product with \(\mathbf {m}^n\), and \(K\) is the stiffness matrix (see Figure 5.3). In addition to the indefiniteness of the overall saddle-point structure, the skew-symmetric contribution renders \(A\) itself non-symmetric, necessitating
GMRES [164] or other Krylov methods for iterative solution.
The linear algebra routines are provided by the MFEM library [144, 145].
5.3.6 Constraint Enforcement and Magnetization Update
The tangent-plane constraint \(\mathbf {v} \cdot \mathbf {m}^n = 0\) is enforced via Lagrange multipliers by solving the full saddle-point system (5.31). This yields a velocity \(\mathbf {v}\) that satisfies the constraint in the finite-element sense, ensuring that the magnetization update remains tangent to the
unit sphere.
After solving for the velocity \(\mathbf {v}\), the magnetization is updated:
Even with exact tangent-plane enforcement, this first-order update does not preserve \(|\mathbf {m}| = 1\) exactly and therefore the deviation is \(\mathcal {O}(\Delta t^2)\) (as shown in Section 5.1.2).
Figure 5.3: Block structure of the tangent-plane scheme saddle-point system corresponding to (5.31). The \((1,1)\)-block \(A\)
combines the damping, precession, and implicit exchange terms. The off-diagonal blocks \(B\) and \(B^T\) enforce the tangent-plane constraint \(\mathbf {v} \cdot \mathbf {m} = 0\), while the solution vector contains both
the velocity \(\mathbf {v}\) and Lagrange multipliers \(\boldsymbol {\lambda }\).
