Physical test scenario 5 - single spin Stoner-Wohlfarth hysteresis simulation

[AdamsMP92]

Motivation

The quasistatic Stoner-Wohlfarth model provides an analytical reference for single-spin hysteresis simulations.

The test compares an Ubermag field sweep with the metastable equilibrium branch of a uniaxial macrospin. It validates:

• uniaxial anisotropy;
• the Zeeman term;
• the conversion between magnetic flux density and magnetic field;
• continuation of metastable states during a field sweep;
• and the angular dependence of the switching field.

Physical model

Consider a normalised magnetic moment confined to the (xz)-plane,

$$\mathbf{m}(\tau) = \begin{bmatrix} \sin\tau\\\ 0\\\ \cos\tau \end{bmatrix},$$

with a uniaxial easy axis

$$\mathbf{u} = \begin{bmatrix} \sin\beta\\\ 0\\\ \cos\beta \end{bmatrix}.$$

The applied magnetic flux density is oriented along the (z)-axis,

$$\mathbf{B}=B\mathbf{e}_z.$$

The Stoner-Wohlfarth energy density is

$$\mathcal{E}(\tau;B) = -K_\mathrm{u}\cos^2(\tau-\beta) - M_\mathrm{s}B\cos\tau.$$

Stationary states satisfy

$$\frac{\partial\mathcal{E}}{\partial\tau} = K_\mathrm{u}\sin\left[2(\tau-\beta)\right] + M_\mathrm{s}B\sin\tau = 0.$$

Introducing the reduced field

$$b = \frac{M_\mathrm{s}B}{2K_\mathrm{u}},$$

the stationary-state equation becomes

$$\frac{1}{2} \sin\left[2(\tau-\beta)\right] + b\sin\tau = 0.$$

A stationary state is locally stable when

$$\frac{\partial^2\mathcal{E}}{\partial\tau^2} = 2K_\mathrm{u}\cos\left[2(\tau-\beta)\right] + M_\mathrm{s}B\cos\tau > 0,$$

or equivalently,

$$\cos\left[2(\tau-\beta)\right] + b\cos\tau > 0.$$

The longitudinal magnetisation is

$$m_z=\cos\tau.$$

Metastable branch continuation

For a quasistatic field sweep, the physical hysteresis branch is obtained by selecting, at every field value, the stable stationary solution closest to the solution at the preceding field value.

The field sequence should run from positive saturation to negative saturation and back:

$$+B_{\max} \rightarrow -B_{\max} \rightarrow +B_{\max}.$$

This continuation follows a local energy minimum until that minimum disappears at the Stoner-Wohlfarth switching field.

Analytical reference implementation

import numpy as np


def stable_tau_nearest_previous(
    b,
    beta,
    tau_previous,
    tol=1e-10,
):
    """Return the stable stationary angle nearest the previous state.

    Solves

        0.5 * sin(2 * (tau - beta)) + b * sin(tau) = 0

    using the companion-matrix representation.
    """
    companion = np.array(
        [
            [0, 1, 0, 0],
            [0, 0, 1, 0],
            [0, 0, 0, 1],
            [
                np.exp(4j * beta),
                2 * b * np.exp(2j * beta),
                0,
                -2 * b * np.exp(2j * beta),
            ],
        ],
        dtype=complex,
    )

    roots = np.linalg.eigvals(companion)

    unit_circle = np.abs(np.abs(roots) - 1) < tol
    tau = np.angle(roots[unit_circle])

    tau = np.concatenate(
        (
            tau - 4 * np.pi,
            tau - 2 * np.pi,
            tau,
            tau + 2 * np.pi,
            tau + 4 * np.pi,
        )
    )

    stable = (
        np.cos(2 * (tau - beta))
        + b * np.cos(tau)
        > 0
    )
    tau = tau[stable]

    if tau.size == 0:
        raise RuntimeError("No stable stationary state found.")

    return tau[np.argmin((tau - tau_previous) ** 2)]


def analytical_sw_hysteresis(
    Ku,
    beta,
    Ms,
    Bmax,
    n,
):
    """Return the quasistatic Stoner-Wohlfarth hysteresis loop."""
    descending = np.linspace(Bmax, -Bmax, n)
    ascending = np.linspace(-Bmax, Bmax, n)
    B = np.concatenate((descending, ascending[1:]))

    tau = np.empty_like(B)
    mz = np.empty_like(B)

    tau_previous = 0.0

    for i, Bi in enumerate(B):
        b = Ms * Bi / (2 * Ku)

        tau[i] = stable_tau_nearest_previous(
            b=b,
            beta=beta,
            tau_previous=tau_previous,
        )

        mz[i] = np.cos(tau[i])
        tau_previous = tau[i]

    return B, mz, tau

Proposed Ubermag setup

A single-cell model can be used because the Stoner-Wohlfarth model assumes coherent rotation.

import numpy as np
import discretisedfield as df
import micromagneticmodel as mm
import oommfc as oc

Ms = 1.7e6
Ku = 4.8e4
beta = np.deg2rad(50)

Bmax = 1.0
n = 250

u = (
    np.sin(beta),
    0.0,
    np.cos(beta),
)

mesh = df.Mesh(
    p1=(0, 0, 0),
    p2=(1e-9, 1e-9, 1e-9),
    n=(1, 1, 1),
)

system = mm.System(name="stoner_wohlfarth")
system.m = df.Field(
    mesh,
    nvdim=3,
    value=(0, 0, 1),
    norm=Ms,
)

At each field value (B_i), define

H = (0, 0, Bi / mm.consts.mu0)

system.energy = (
    mm.UniaxialAnisotropy(K=Ku, u=u)
    + mm.Zeeman(H=H)
)

The factor

$$H_i=\frac{B_i}{\mu_0}$$

is required because mm.Zeeman expects the magnetic field in (\mathrm{A,m^{-1}}), whereas the analytical reference uses (B) in tesla.

The system should be minimised sequentially over the complete field sweep. The relaxed state at field (B_i) must be used as the initial state at field (B_{i+1}) so that the numerical simulation follows the same metastable branch as the analytical model.

No demagnetisation term should be included because it is not part of this reference model.

Primary assertion

The numerical longitudinal magnetisation should be compared with the analytical hysteresis branch:

B_ref, mz_ref, tau_ref = analytical_sw_hysteresis(
    Ku=Ku,
    beta=beta,
    Ms=Ms,
    Bmax=Bmax,
    n=n,
)

np.testing.assert_allclose(
    mz_simulated,
    mz_ref,
    rtol=...,
    atol=...,
)

Values immediately adjacent to a switching event may be excluded from the strict pointwise comparison because the analytical and numerical switching fields are sampled on a discrete field grid.

Away from the switching points, the relaxed numerical state should agree directly with the analytical metastable solution.

Switching-field check

The Stoner-Wohlfarth switching-field magnitude is

$$B_\mathrm{sw}(\beta) = \frac{2K_\mathrm{u}}{M_\mathrm{s}} \left( \lvert\cos\beta\rvert^{2/3} + \lvert\sin\beta\rvert^{2/3} \right)^{-3/2}.$$

The numerical switching field can be identified from the discontinuity in (m_z). Its error should be bounded by the field-step size,

$$\left| B_\mathrm{sw}^{\mathrm{sim}} - B_\mathrm{sw}^{\mathrm{ref}} \right| \lesssim \Delta B.$$

A parameterised test over several angles can additionally verify the Stoner-Wohlfarth astroid.

Special symmetry directions such as (\beta=0) and (\beta=\pi/2) may be tested separately because they contain degenerate stationary branches.

Suggested assertions

The test can verify:

1. agreement of the complete metastable branch away from switching points;
2. agreement of the descending and ascending switching fields;
3. inversion symmetry of the hysteresis loop;
4. saturation at sufficiently large positive and negative fields;
5. the angular dependence of the switching field.

For an unbiased uniaxial model, the loop should satisfy approximately

$$m_z(-B) = -m_z(B)$$

when corresponding points on opposite sweep branches are compared.

Expected value

This test provides a quantitative reference for quasistatic single-spin hysteresis.

It would validate the combined behaviour of the uniaxial anisotropy, Zeeman energy, field-unit conversion, minimisation procedure, and metastable-state continuation. It also provides a direct angular test through the analytical Stoner-Wohlfarth switching astroid.