Hysteresis loop of a magnetic nanoparticle in an AC field

What is computed

A single-domain nanoparticle with uniaxial anisotropy sits in an alternating field H(t) = H₀ sin(2πf t) applied along its easy axis. The particle is immobilised (no mechanical or Brownian rotation of the particle body) and does not interact with neighbours, so the only degree of freedom is the direction of its magnetic moment μ = M_sV, which fluctuates thermally between the two easy directions and is driven back and forth by the field. The applet solves Brown’s Fokker–Planck equation for the probability density of that direction, live in your browser (nothing is precomputed, nothing runs on a server), and integrates the dynamic hysteresis loop of the ensemble average m(t) = ⟨cos θ⟩ over one field period. The area of the loop is the energy dissipated per cycle; from it follow the specific absorption rate SAR and the efficiency ε relative to the bound SAR_max.

The energy landscape

With θ the angle between the moment and the easy axis (which is also the field axis), the energy is the sum of anisotropy and Zeeman terms,

E(θ) = K₁ V sin²θ − μ₀ M_s V H(t) cos θ

In units of k_BT, with x = cos θ, the reduced field h = H/H_K, H_K = 2K₁/(μ₀M_s) and the reduced barrier σ = K₁V/k_BT:

E/k_BT = σ sin²θ − 2σ h cos θ = σ(1 − x²) − 2σ h x

For |h| < 1 there are two minima, θ = 0 (x = +1, energy −2σh) and θ = π (x = −1, energy +2σh), separated by a maximum at x = −h. Seen from the well aligned with the field the barrier is σ(1 + h)², seen from the well opposed to the field it is σ(1 − h)²: the field lowers one barrier and raises the other, which is what drives the moment across. At h ≥ 1 the opposed minimum disappears altogether and the moment switches deterministically — the Stoner–Wohlfarth limit, which the loop approaches for large σ (large particles, low temperature, high frequency).

The unknown: W(x, t)

W(x, t) is a probability density in one orientation variable, x = cos θ ∈ [−1, 1], and time: W(x,t) dx is the probability that the projection cos θ of the moment onto the axis lies in [x, x+dx] at time t. Because the field is parallel to the easy axis the problem is axially symmetric: the azimuth φ of the moment carries no information and has been integrated out, so W is a function of just two variables and is normalised, ∫₋₁⁺¹ W dx = 1. W is not the magnetisation; the magnetisation is its first moment m(t) = ∫ x W(x,t) dx = ⟨cos θ⟩ = M(t)/M_s. The heat map in the panel “Distribution W(x, t)” shows W itself; the loop is m plotted against h.

The equation

∂W/∂t = (1 / 2τ_N) ∂/∂x [ (1 − x²) ( ∂W/∂x + W ∂(E/k_BT)/∂x ) ]

This is a drift–diffusion (Smoluchowski) equation on the unit sphere. The first term in the inner bracket is diffusion: thermal agitation spreading the orientation. The second is drift: the torque −∂E/∂θ pulling the moment towards the energy minima. The factor (1 − x²) = sin²θ is the metric of the sphere in the variable x; it makes the flux vanish at the poles, so probability is conserved. The single time constant is the free rotational diffusion time of the moment,

τ_N = (1 + α²) M_s V / (2 α γ k_BT) = σ τ₀,   τ₀ = (1 + α²) / (α γ B_ani),   B_ani = 2K₁/M_s,   γ = 1.760859×10¹¹ rad s⁻¹ T⁻¹

Brown (1963) derived this equation as the exact orientation-space form of the stochastic Landau–Lifshitz–Gilbert equation for this geometry — it is the same physics as a thermal LLG simulation of one macrospin, but for the full ensemble at once and without noise. The smallest non-zero eigenvalue of the operator at fixed field is the Néel relaxation rate 1/τ; the familiar two-level (Arrhenius) master equation, τ ∝ τ₀ exp[σ(1 ∓ h)²], is its high-barrier limit. Here no such approximation is made: the full distribution is propagated for any σ, h and frequency, including the superparamagnetic and the Stoner–Wohlfarth regimes and everything in between.

Symbols

Beyond the axially symmetric case (sphere modes)

If the easy axis is tilted by ψ from the field, or the anisotropy is cubic, the azimuth no longer drops out and W depends on both angles. The mode selector then switches to Brown’s equation on the full unit sphere,

∂W/∂t = −∇_s·(v W) + D ∇_s²W,  v = −γ' m×B − αγ' m×(m×B),  D = αγk_BT/((1+α²)M_sV) = 1/(2τ_N),  B = −(1/M_s)∂F/∂m + μ₀H(t) n

The first term of v is the precession about the effective field. In the axially symmetric case it only rotates W about the axis and drops out; with tilted or cubic minima it couples the two angles and is essential (it is the reason why a tilted particle is not simply a projected 1-D one). Numerically, damping drift and diffusion use Scharfetter–Gummel fluxes, the precession uses central fluxes (upwinding it would add a large artificial diffusion when the minima lie away from the grid poles). Central fluxes are not monotone: when an energy well becomes narrower than about two cells (q = Δθ√σ > 0.6) W can turn negative, and the status line flags such results.

Cubic anisotropy: F = K₁(m_x²m_y² + m_y²m_z² + m_z²m_x²) + K₂m_x²m_y²m_z². For K₁ < 0 (magnetite, K₁ ≈ −1.1×10⁴ J/m³) the eight ⟨111⟩ directions are minima and the barrier between neighbouring minima (over ⟨110⟩) is only |K₁|V/12, compared with KV for a uniaxial particle — a cubic magnetite particle of the same volume is therefore much more superparamagnetic. The field direction in the crystal frame can be chosen ([111], [100], [110]) or averaged over six quasi-uniform (Fibonacci) directions. M_s = 400 kA/m is the paper’s value; bulk magnetite has about 480 kA/m. For tilted uniaxial particles see I. S. Poperechny, Yu. L. Raikher, V. I. Stepanov, Phys. Rev. B 82, 174423 (2010).

Numerics in the browser: finite volumes on an N_θ×2N_θ grid (default 32×64), implicit Euler with 400 steps per period; each step is solved exactly with a band LU (bandwidth 2N_θ; the unresolved precession, ~10 rad per step, defeats Krylov solvers), the factors of the S/2+1 distinct step matrices are kept in float32 (~110 MB at 32×64) and every solve is finished by float64 iterative refinement; the periodic state comes from GMRES on (I − R M_half)W = 0 with R the inversion m → −m (typically two iterations). A loop takes 3–6 s (×6 for the orientation average).

Numerics — what “live” means

The interval x ∈ [−1, 1] is divided into N finite-volume cells uniform in θ (fine near the poles, where the population sits at large σ). Fluxes between cells use the Scharfetter–Gummel (exponentially fitted) scheme, so the discrete equilibrium is exactly the Boltzmann distribution for any σ. Time stepping is Crank–Nicolson with S steps per period, one tridiagonal solve per step. The periodic steady state is not reached by cycling until transients die out: the propagator over half a period, M_half, is assembled column by column, and the symmetry of the equation under (x, h) → (−x, −h) turns the periodicity condition into W* = R M_half W* with R the reflection x → −x. One dense N×N linear solve, with one row replaced by the normalisation, gives W*. From W* one period is propagated to record m(t), h(t) and W(x, t), whose normalisation is checked at every step (see the status line). With the defaults (N = 120, S = 600) this takes about 60 ms in a Web Worker, so every input change recomputes immediately. Crank–Nicolson is A-stable but not L-stable: when a large field (h₀ > 1) sweeps the barrier of the opposed well from ≈25 k_BT to zero within only a few time steps, the switching front is unresolved in time and the solution rings (negative W, |m| > 1). The applet therefore counts the steps in that sweep window, n = S √(25/σ) / (2π √(h₀²−1)), and raises S automatically to the smallest multiple of 600 (at most 24 000) that gives n ≥ 20; the status line says so when it happens.

Hysteresis loops: stochastic LLG (solid) vs Fokker–Planck (dashed)

Each panel shows the steady-state loop of the ensemble-averaged magnetisation along the field; the legend gives the loop area / 4B₀ computed from the plotted curve. M_s = 400 kA/m (particles) or 800 kA/m (disks), α = 0.1, T = 300 K throughout.

Cubic points in numbers

The full-sphere Fokker–Planck mode was checked against stochastic Landau–Lifshitz–Gilbert simulations with jaxmag (matrix-free cubic anisotropy operator make_cubic_anisotropy, ensembles of 8192 independent macrospins, Heun scheme, Δt = 1 ps, 4 field periods with the first discarded, run on the A100 GPUs of gto4). Material: cubic magnetite K₁ = −1.1×10⁴ J/m³, K₂ = 0, M_s = 4×10⁵ A/m, α = 0.1, T = 300 K, field on the discomfort limit H₀f = C = 5×10⁹ A m⁻¹s⁻¹. Fokker–Planck on a 96 × 192 grid, 600 steps per period.

(a) Symmetry test: a field along x and along z must give the same loop. The grid prefers z (its poles); the difference shrinks from 3.5 % on 64 × 128 to 2 % on 96 × 192.

(b) Resolution: on the coarse browser grids the 60 nm points are not resolved (narrow wells away from the grid poles; the applet flags them red). Use the default browser grid for D ≤ 40 nm in cubic mode and the finer grid option for larger particles.

(c) Physics: pure cubic magnetite loses far less than a uniaxial particle with K = |K₁| (25 nm: 0.004–0.006 vs 0.88; 40 nm: 0.07–0.26 vs 0.95), because the barrier between neighbouring ⟨111⟩ minima is only |K₁|V/12. This agrees with N. A. Usov, M. S. Nesmeyanov, E. M. Gubanova, N. B. Epshtein, Beilstein J. Nanotechnol. 10, 305 (2019), doi:10.3762/bjnano.10.29, who find appreciable SAR for cubic magnetite only at 45–55 nm; real particles that heat well have additional shape or surface anisotropy (combined anisotropy: N. A. Usov, J. M. Barandiarán, J. Appl. Phys. 112, 053915 (2012), doi:10.1063/1.4749799).

Code: jaxmag, gitlab.com/dietersuess/jaxmag, branch feature/md, examples/cubic_macrospin_loop.py.