\documentclass[fleqn,10pt]{wlscirep}

\usepackage[utf8]{inputenc} \usepackage[T1]{fontenc} \usepackage{amsmath,amssymb} \usepackage{siunitx} \usepackage{booktabs} \usepackage{xcolor} \usepackage{graphicx} \usepackage{orcidlink} \usepackage{acronym}
\usepackage{xstring}
\let\bareurl\url
\renewcommand{\url}[1]{\IfBeginWith{#1}{10.}{\bareurl{https://doi.org/#1}}{\bareurl{#1}}}

\sisetup{ per-mode=power, inter-unit-product=\ensuremath{{}\cdot{}}, range-units=single, range-phrase=~to~, separate-uncertainty=true, retain-explicit-plus=true, }
\DeclareSIUnit{\sample}{S}

\newcommand{\Ms}{M_\mathrm{s}} \newcommand{\Man}{M_\mathrm{an}} \newcommand{\Mirr}{M_\mathrm{irr}} \newcommand{\He}{H_\mathrm{e}} \newcommand{\Hc}{H_\mathrm{c}} \newcommand{\Br}{B_\mathrm{r}} \newcommand{\Ku}{K_\mathrm{u}} \newcommand{\Kone}{K_1} \newcommand{\tauR}{\tau_\mathrm{R}} \newcommand{\Dcoef}{D} \newcommand{\DTHM}{\textsc{DTHM}} \newcommand{\DTHMS}{\textsc{DTHM Simplified}} \newcommand{\JA}{J--A} \newcommand{\DDD}{\textsc{D-D-D}} \newcommand{\PI}{\textsc{PI}}

\acrodef{FP}{Fokker--Planck}
\acrodef{DTHMS}[DTHM--S]{simplified diffusion-driven thermodynamic hysteresis model}
\acrodef{DTHM}[DTHM]{diffusion-driven thermodynamic hysteresis model}
\acrodef{GO}{grain-oriented}
\acrodef{NGO}{non-grain-oriented}
\acrodef{HGO}{high-grain-oriented}
\acrodef{RD}{rolling direction}
\acrodef{TD}{transverse direction}
\acrodef{JA}[J--A]{Jiles--Atherton}
\acrodef{DDD}{D'Aloia--Di~Francesco--De~Santis}
\acrodef{PI}{Prandtl--Ishlinskii}
\acrodef{ODE}{ordinary differential equation}
\acrodef{FEM}{finite-element method}
\acrodef{CFL}{Courant--Friedrichs--Lewy}
\acrodef{RMSE}{root-mean-square error}
\acrodef{TPE}{tree-structured Parzen estimator}
\acrodef{LLG}{Landau--Lifshitz--Gilbert}
\acrodef{EBSD}{electron backscatter diffraction}
\acrodef{ODF}{orientation distribution function}

\renewcommand{\thefigure}{S\arabic{figure}} \renewcommand{\thetable}{S\arabic{table}} \renewcommand{\theequation}{S\arabic{equation}} \renewcommand{\thesection}{} \setcounter{section}{0}

\title{Supplementary Information for:\\ A unified rate-dependent anisotropic hysteresis model for grain-oriented, non-oriented, and directional electrical steels}

\author[1,*]{Mario Bendra\,\orcidlink{0000-0001-5994-2682}} \author[1]{Florian Poltschak\,\orcidlink{0000-0002-9618-0901}}

\affil[1]{JKU Hoerbiger Research Institute for Smart Actuators, Altenberger Strasse 69, Science Park 1, 4040 Linz, Austria}

\affil[*]{mario.bendra@jku.at}

\begin{abstract}
This document is the Supplementary Information for the manuscript referenced above. It contains tables of fitted parameters for every material studied, per-loop residuals, the moment-closure verification against the Ramesh--Jiles--Bi 1997 anisotropic-hysteresis benchmark, the cross-model overlays for the Li and Luo amplitude sweep, the discrete Boltzmann convergence study, and the parameter-identifiability map. All numerical results are reproducible from the public \texttt{hysterpy} code repository cited in the main text.
\end{abstract}

\begin{document}
\flushbottom
\maketitle
\thispagestyle{empty}

\noindent This Supplementary Information accompanies the main
manuscript and contains: a model-naming legend (\textbf{Note~S1}) clarifying which model is used where, fitted parameter sets for every material (\textbf{Tables~S1--S4}), per-loop residuals and identifiability summaries (\textbf{Tables~S5--S6}), the per-ring cross-model comparison, unification ablation, and cross-validation summary (\textbf{Tables~S7--S9}), supporting figures (\textbf{Figures~S1--S13}: $S^{2}$ geometry, unification ablation, the \qty{1}{\hertz} fit, cross-model overlays, the cost-landscape multimodality study, the cross-model \ac{RMSE} summary, the multistart-basin scatter, the 1D Maxwell skin-depth verification, the Ramesh--Jiles--Bi consistency check, the third-harmonic excitation demonstration, the HEFMAG external-dataset blind cross-validation, and the \ac{LLG}-direction infrastructure validation), three further explanatory notes (\textbf{Notes~S2--S4}) on the MagHyst$^\circledR$ measurement protocol, the lat--lon grid convergence, and the rationale for fitting the \qty{1}{\hertz} quasi-static loop separately, a cross-validation and cost-landscape-identifiability study (\textbf{Note~S5}), and a 1D Maxwell-coupled derivation of the cubic-anisotropy sigmoid corner frequency (\textbf{Note~S6}, \textbf{Fig.~S9}). All numerical results are reproducible from the public code repository (\texttt{hysterpy}, MIT license) cited in the main text Data Availability statement.

\section*{Supplementary Note S1: Model naming and where each model is used}

The paper compares the present \emph{vector \acf{FP}} model against six published baselines. The naming convention used throughout, and the hysterpy class behind each, is listed below. The main result of the paper relies solely on the first entry. The rest appear as cross-model baselines in main-text Tables~3 and~4, in Supplementary Table~\ref{tab:s8_models_per_ring}, and in the cross-model overlays of Figs.~\ref{fig:s_models_hgo} and~\ref{fig:s_models_ring}.

\begin{itemize}
	\item \textbf{vector~FP $S^{2}$} (\texttt{hysterpy.DTHMS2Model}): the present model. Fokker--Planck equation on the unit sphere $S^{2}$, Chang--Cooper well-balanced flux, Zeeman + uniaxial $\Ku$ + cubic $\Kone(f)$ free energy, polycrystalline orientation distribution. Used for: main-text Fig.~3, main-text Table~2, the ``This work'' column of main-text Table~3. The lab-frame versus crystal-frame geometry underlying the Goss-textured fit is illustrated in Fig.~\ref{fig:s_geometry}, and the per-frequency knobs and per-loop residuals of the joint fit are listed in Table~\ref{tab:s3_baghel_S2}.
	\item \textbf{scalar~FP $S^{1}$} (\texttt{hysterpy.DTHMModel}): the strict $S^{2}\!\to\!S^{1}$ reduction of the vector~\ac{FP} model: same Chang--Cooper flux, anisotropy reduced to a uniaxial scalar, no cubic invariant, and used for: the rate-independent isotropic reduction on the seven \ac{NGO} rings (Tables~\ref{tab:s4_rings}, \ref{tab:s8_models_per_ring} and Fig.~\ref{fig:s_cross_model_bars}).
	\item \textbf{DTHM--S}\acused{DTHMS} (\texttt{hysterpy.DTHMSimplifiedModel}): the simplified diffusion-driven thermodynamic hysteresis model, a Jiles--Atherton-family \ac{ODE} solver with separate $\Mirr$ tracking. Used as a baseline only.
	\item \textbf{J--A} (\texttt{hysterpy.JilesAthertonModel}): classical Jiles--Atherton \ac{ODE}~\cite{jiles1986}.
	\item \textbf{Chwastek} (\texttt{hysterpy.ChwastekModel}): modified J--A with mixing-rule $\beta$ replacing $c$~\cite{chwastek2009}.
	\item \textbf{DDD}\acused{DDD} (\texttt{hysterpy.DDDModel}): D'Aloia--Di~Francesco--De~Santis elastic-collision model~\cite{daloia2021}.
	\item \textbf{Preisach} (\texttt{hysterpy.PreisachModel}): classical Preisach with Gaussian hysteron density~\cite{mayergoyz}.
	\item \textbf{PI}\acused{PI} (\texttt{hysterpy.PrandtlIshlinskiiModel}): Prandtl--Ishlinskii play-operator pile~\cite{krasnoselskii}.
\end{itemize}

Throughout this paper, ``this work'' denotes the vector Fokker--Planck $S^{2}$ model with the parameter set of main-text Table~2. The ``isotropic reduction'' on the \ac{NGO} rings is the same model with $\Ku\!=\!\Kone\!=\!0$ enforced, which is mathematically identical to scalar \ac{FP} on $S^{1}$. The joint single-parameter-set fit of the Li and Luo amplitude sweep is tabulated in Table~\ref{tab:s5_liluo}. The per-grain \ac{LLG}-direction infrastructure released alongside the models above (\texttt{hysterpy.llg\_grain}) is validated in Fig.~\ref{fig:llg_infrastructure}.

\section*{Supplementary Note S2: MagHyst$^\circledR$ measurement protocol}

The seven non-oriented ring-core measurements reported in the main text were acquired in-house on the MagHyst$^\circledR$ Modular system of Ilmenauer Mechatronik GmbH~\cite{maghyst}. The closed magnetic circuit consists of a toroidal ring sample with primary and secondary windings. The primary is current-driven by a software-defined reference $i_\mathrm{prim}(t)$ following a triangular waveform of period $T = 1/f$, the secondary picks up the induced voltage $u_\mathrm{ind}(t) = -\mathrm{d}\Psi/\mathrm{d}t$, and the flux linkage is recovered by numerical integration of the open-circuit secondary voltage. Sample geometry parameters (mean radius $r_\mathrm{Fe}$, iron cross-section $A_\mathrm{Fe}$, primary turns $N_\mathrm{prim}$, secondary turns $N_\mathrm{sec}$) translate $(i_\mathrm{prim}, \Psi)$ to material-frame $(H, B) = (N_\mathrm{prim}\,i_\mathrm{prim}/(2\pi r_\mathrm{Fe}), \Psi/(N_\mathrm{sec} A_\mathrm{Fe}))$. The MagHyst$^\circledR$ device exports the measured loop as an XML record containing the full $(t, i_\mathrm{prim}(t), u_\mathrm{ind}(t), B(t), H(t))$ time series at a sampling rate of typically \qty{100}{\kilo\sample\per\second}. The measurement frequency is the inverse of the four-cycle MeasureTime field in the XML. We use the third loop of each four-cycle record for fitting, discarding the first two as warmup. This is consistent with the $n_\mathrm{warmup} = 2$ convention used in all model simulations~\cite{zirka2008}. The dynamic response of the fitted X10CrAlSi7 ring under multi-harmonic distorted excitation is demonstrated in Fig.~\ref{fig:s_harmonic}.

\section*{Supplementary Note S3: Lat--lon grid convergence study}

The Chang--Cooper well-balanced flux discretization is exact in the sense that the discrete Boltzmann distribution is the exact kernel of the discrete Fokker--Planck operator, regardless of grid resolution. Convergence to the continuum limit is therefore a question of grid quadrature accuracy, not of stability or operator bias. The convergence of the discrete $\langle m_x \rangle = M/\Ms$ at $H/a = 1$ to the continuum Langevin function $\mathcal{L}(1) = 0.3130$ as a function of $N_\theta$ at fixed $N_\phi = 2 N_\theta$ is summarized numerically below:
\begin{equation}
	N_\theta \in \{16, 24, 32, 48, 64\} \Rightarrow
	\frac{|\langle m_x\rangle_\mathrm{disc} - \mathcal{L}(1)|}{\mathcal{L}(1)}
	\in \{1.6\!\times\!10^{-3},\,5\!\times\!10^{-4},\,3\!\times\!10^{-4},\,
	1\!\times\!10^{-4},\,5\!\times\!10^{-5}\}.
	\label{eq:convergence}
\end{equation}
The reported main-text $S^2$ results use $(N_\theta, N_\phi) = (32, 64)$, which gives 0.03\,\% accuracy in the Langevin recovery and is sufficient for all materials in our validation set ($H/a < 5$). Sharper-distribution problems (e.g.\ ferromagnets near the Curie point or low-field limit of permanent magnets) may require higher $N_\theta$. The companion code-equivalence check against the Ramesh--Jiles--Bi anisotropic anhysteretic benchmark, computed at the same grid resolution, is reported in Fig.~\ref{fig:angle_resolved}.

\section*{Supplementary Note S4: Quasi-static \qty{1}{\hertz} loop fit}

The \qty{1}{\hertz} Baghel measurement is fitted separately from the \qtyrange{50}{200}{\hertz} dynamic loops. The reason is geometric. At fixed $\tauR \Dcoef = 1$, the canonical Boltzmann limit is reached when $\omega \tauR \ll 1$ (period much longer than relaxation). For our fitted $\tauR$ values ($\tauR \approx 0.02$~s in dynamic fits), $\omega \tauR \approx 0.13$ at \qty{1}{\hertz} versus $\omega \tauR \approx 6.3$ at \qty{50}{\hertz}: the \qty{1}{\hertz} loop is essentially anhysteretic. The same $(\Ms, a, \alpha, \Ku, \Kone)$ shared parameters that fit the dynamic loops are found by an independent fit of the \qty{1}{\hertz} loop with $\tauR \to 0^+$ in the model. This split is consistent with established practice in the dynamic Jiles--Atherton literature~\cite{zirka2008,baghel2015} and with the Kalmykov--Coffey analysis of the orientational Fokker--Planck relaxation~\cite{coffey2004}. The fitted parameters and residuals of this quasi-static fit are reported in Table~\ref{tab:s1_baghel_qs} and overlaid on the measurement in Fig.~\ref{fig:s_qs}. The per-loop residuals of the corresponding direction-resolved baseline across all five frequencies are collected in Table~\ref{tab:s6_per_loop}.

\section*{Supplementary Note S5: Generalization across frequency splits, $K_1(f)$ functional form, and cost-landscape identifiability}

The main result of the main text (Methods subsection ``Joint fit across directions, amplitudes, and frequencies'', Results Fig.~3) shares a base $(\Ms, a, \alpha, \Ku, D)$ plus the three sigmoid constants $(K_1^{\mathrm{low}}, f_\mathrm{c}, n)$ across all eight \ac{RD} and \ac{TD} dynamic loops at \qtyrange{50}{200}{\hertz}. We test three orthogonal robustness questions on the same dataset and the same shared-base architecture: (a) does the fit generalize across frequency splits? (b) is the sigmoid form for $K_1(f)$ defensible against alternative monotone parametrizations? (c) does the 38-dimensional cost surface admit unique identification, or are there competing local minima? The complementary unification ablation on the same dataset is reported in Table~\ref{tab:s9_unification_ablation} and visualized in Fig.~\ref{fig:s_unification}.

\paragraph*{Cross-frequency generalisation.} Two frequency-split runs probe extrapolation in both directions. A \emph{forward-extrapolation} fit jointly fits the \qty{50}{\hertz} \ac{RD}$+$\ac{TD} pair only and predicts the higher frequencies via the $K_1(f)$ sigmoid, the $\tauR(f) = \tauR^{50}(50/f)$ scaling, and $H_\mathrm{pin}^{\mathrm{TD}}(f) = H_\mathrm{pin}^{50}(f/50)$ of the main text. Mean \ac{RMSE} at \qty{50}{\hertz} is \qty{0.043}{\tesla}. A \emph{reverse-extrapolation} fit holds the shared base fixed and refits only the four $K_1(f)$ sigmoid constants $(K_1^{\mathrm{low}}, f_\mathrm{c}, n, \tau_{\mathrm{const}})$ on the \qtyrange{100}{200}{\hertz} loops, then predicts \qty{50}{\hertz}. Mean fit cost at the held-in frequencies is \qty{0.072}{\tesla}. Both splits reproduce the held-out frequencies at residuals comparable to the main joint fit (Table~\ref{tab:cv_summary}), confirming that the $K_1(f)$ sigmoid is a transferable scaling rather than a frequency-by-frequency refit in disguise.

\paragraph*{Alternative $K_1(f)$ functional forms.} The same reverse-extrapolation split was re-run with the sigmoid replaced by a monotonically \emph{increasing} Hill function $K_1(f) = K_1^{\mathrm{max}} f^n / (f_0^n + f^n)$. The Hill form lands at a fit cost \qty{0.092}{\tesla} versus the sigmoid's \qty{0.072}{\tesla} (a $28\,\%$ degradation under identical budget). This selects the monotonically decreasing sigmoid form as the data-preferred functional shape, consistent with the skin-depth-mediated cubic-anisotropy washout interpretation of the main text (Methods subsection on the frequency-dependent effective cubic anisotropy).

\paragraph*{Leave-one-out cross-validation.} A 4-fold leave-one-out cross-validation of the sigmoid form refits $(K_1^{\mathrm{low}}, f_\mathrm{c}, n, \tau_{\mathrm{const}})$ on three of the four frequencies and tests on the held-out one. Per-fold fit costs range \qtyrange{0.062}{0.077}{\tesla}. The sigmoid constants themselves span $K_1^{\mathrm{low}} \in \qtyrange{1098}{5655}{\joule\per\meter\cubed}$ and $f_\mathrm{c} \in \qtyrange{37}{182}{\hertz}$ across the four folds, reflecting the well-known three-parameter sigmoid degeneracy: the physically observable quantity is the curve $K_1(f)$, not the three parameters that index it. Evaluated at the measurement frequencies, the leave-one-out $K_1(f)$ predictions all sit within the main-fit trajectory $\pm 30\,\%$, supporting the sigmoid as a phenomenological surrogate for the cubic-anisotropy washout rather than as a derived intrinsic constant.

\paragraph*{Cost-landscape identifiability.} To quantify practical identifiability, we ran a 20-seed cold-start multistart on the full 38-dimensional joint cost surface (uniform-random initial population over the full parameter bounds, popsize\,$=3 \cdot 38=114$, $\mathrm{maxiter}=40$, total wall $16.94$~h on 24 cores). The resulting 20 best-of-seed parameter sets cluster into three basins (Fig.~\ref{fig:multistart_basins}): (i)~\emph{basin A} (this work, $N\!=\!14$ of $20$): $K_1^{\mathrm{low}} \in \qtyrange{955}{2685}{\joule\per\meter\cubed}$, $f_\mathrm{c} \in \qtyrange{38}{215}{\hertz}$. The within-basin coefficient of variation $\sigma(\cdot)/\langle\cdot\rangle$ is \num{0.010} on $\Ms$, \numrange{0.036}{0.043} on $\tauR(f)$, \num{0.041} on $\omega_{\mathrm{Goss}}$, and \num{0.136} on the physically observable $K_1(\qty{50}{\hertz})$, while the individual sigmoid constants $(K_1^{\mathrm{low}}, f_\mathrm{c}, n)$ span CoV \numrange{0.34}{0.40} along the $K_1(f)$ curve (the same parametric degeneracy as in the cross-validation study above). (ii)~\emph{Basin B} ($N\!=\!2$): elevated cubic anisotropy at $K_1^{\mathrm{low}} \approx \qty{1.2e4}{\joule\per\meter\cubed}$ (within $4\times$ of the Fe-Si single-crystal value \qty{4.8e4}{\joule\per\meter\cubed}) with sharp sigmoid falloff $n\!\approx\!3.8$ that collapses $K_1(\qty{200}{\hertz})$ to within \qty{10}{\joule\per\meter\cubed}. (iii)~\emph{Basin B$'$} ($N\!=\!4$): near-constant $K_1 \in \qtyrange{150}{850}{\joule\per\meter\cubed}$ across the measurement band, with corner frequency $f_\mathrm{c}$ above \qty{400}{\hertz} so the sigmoid never turns off.

The three basins are degenerate to within $\Delta\mathrm{cost}\!\approx\!0.03$ at the multistart budget (best-of-seed costs $0.096$, $0.107$, $0.102$ for A, B, B$'$ respectively), but they are physically distinguishable. Basin A produces the textbook \ac{HGO} Goss-textured \ac{RD}$+$\ac{TD} coupling reported in the main text. Basin B implies single-crystal-magnitude cubic anisotropy that would require $\approx\!4\times$ higher saturation fields to verify than the Baghel dataset spans. Basin B$'$ has $K_1$ effectively switched off, contradicting the Goss-texture cubic structure. Restricting the parameter search to the basin-A neighborhood is therefore a physically motivated identifiability constraint, not a fitting artifact. Restriction to basin A is corroborated by the local Fisher-information identifiability summary at the basin-A optimum in Table~\ref{tab:s7_identifiability}. The Bingham-spread anisotropic-texture variant of basin B (basin~C of Supplementary Fig.~\ref{fig:s_bimodal_basins}) returns the mean \ac{RMSE} \qty{1.1}{\tesla} and is rejected on the goodness-of-fit basis there. All multistart outputs and the basin analysis are reproducible from the public companion repository cited in the main-text Code Availability statement. An external blind cross-validation of the same model class against the open HEFMAG dataset is reported in Fig.~\ref{fig:s_hefmag}.

\section*{Supplementary Note S6: 1D Maxwell-coupled analysis of the cubic-anisotropy sigmoid corner frequency}

The main text introduces a phenomenological frequency-dependent effective cubic anisotropy
\begin{equation}
	K_1(f) \;=\; \frac{K_1^\mathrm{low}}{1 + (f/f_\mathrm{c})^n}, \qquad
	K_1^\mathrm{low}\!=\!\qty{1619}{\joule\per\meter\cubed},
	\quad f_\mathrm{c}\!=\!\qty{99.0}{\hertz},
	\quad n\!=\!2.81,
	\label{eq:k1_sigmoid_supp}
\end{equation}
to close the rolling-direction residual at \qty{150}{\hertz} and \qty{200}{\hertz}. The main-text Methods subsection on the frequency-dependent effective cubic anisotropy gives the closed-form linear skin-depth consistency check: at the published \ac{HGO} operating-knee permeability $\mu_\mathrm{r}\!\approx\!1.35\!\times\!10^{4}$, the lamination ratio $d/\delta\!=\!1$ lands at $f\!\approx\!\qty{99}{\hertz}$, matching the fitted corner. This note reports a direct 1D Maxwell-coupled numerical simulation that confirms the corner-frequency match, quantifies the depth gap between the fitted sigmoid and linear skin theory, and identifies the operative mechanism.

\paragraph*{Setup.} We solve the 1D coupled magnetic-diffusion problem $\sigma\mu_0\,\partial_t(H+M) = \partial^2_z H$ (main-text Methods subsection on the frequency-dependent effective cubic anisotropy) inside one Hi-B 27M-OH lamination of thickness $d\!=\!\qty{0.30}{\milli\meter}$ and bulk conductivity $\sigma\!=\!\qty{2.1e6}{\siemens\per\meter}$. The cross-section is discretized into $N_z\!=\!17$ uniformly spaced slices. Each slice carries an independent rotation-only Fokker--Planck constitutive state on $S^{2}$ (lat--lon grid $16\!\times\!32$) using the shared parameters $(M_\mathrm{s}, a, \alpha, K_\mathrm{u}, D, \tauR)$ of this work and the \emph{same} cubic anisotropy $\Kone\!=\!K_1^\mathrm{low}$ for every slice and every frequency. The driving field $H_\mathrm{lab}(t)\!=\!\qty{250}{\ampere\per\meter}\,\sin(2\pi f t)$ is applied as a Dirichlet boundary condition at $z\!=\!\pm d/2$. Time integration uses operator-splitting: per step the per-slice Fokker--Planck operator updates $M^{n+1}(z)$ from $H^n(z)$, then a backward-Euler tridiagonal solve advances $H^{n+1}(z)$ given the new $M^{n+1}(z)$ and the boundary value. We take \num{160} time steps per period and report the third period after two warm-up cycles.

\paragraph*{Spatial field profile.} Supplementary Fig.~\ref{fig:maxwell_skin_depth}(a) shows the steady-state $H_\mathrm{internal}(z)/H_\mathrm{lab,peak}$ at the surface-field peak, for two representative frequencies. At \qty{50}{\hertz} the interior field is nearly uniform across the lamination ($\le 1\,\%$ depression at the center). At \qty{200}{\hertz} the centreline value drops to $\approx 0.84\,H_\mathrm{lab,peak}$, the signature of incipient skin shielding at $d/\delta\!\approx\!1.5$.

\paragraph*{Corner-frequency match.} The physically operational quantity is the bulk-averaged attenuation \\
$A_\mathrm{num}(f) = \max_t |\langle H(z,t)\rangle_z| / \max_t |H_\mathrm{lab}(t)|$. Across the six measured frequencies, the numerical $A_\mathrm{num}(f)$ matches the analytical linear-theory envelope $|\tanh(\gamma d/2)/(\gamma d/2)|$ with $\gamma\!=\!(1+j)/\delta$ to better than $1\,\%$, validating the numerical solver against textbook Maxwell-diffusion physics (Fig.~\ref{fig:maxwell_skin_depth}(b), blue dots vs grey dashed). The half-attenuation point of $A_\mathrm{num}(f)$ lies in the same decade as the fitted phenomenological corner $f_\mathrm{c}\!=\!\qty{99.0}{\hertz}$. Since both $A_\mathrm{num}(f)$ and the $K_1(f)$ sigmoid are set by the same dimensional combination $\mu_\mathrm{r}\,\sigma\,d^{2}\,\omega$ of the lamination, the corner-frequency location is fixed by this shared material time scale rather than being an independent fit degree of freedom. The two mechanisms nonetheless differ in magnitude (paragraph below).

\paragraph*{Depth gap is real and informative.} The fitted sigmoid $K_1(f)/K_1^\mathrm{low}$ (Fig.~\ref{fig:maxwell_skin_depth}(b), orange solid) decays significantly faster than $A_\mathrm{num}(f)$ beyond the corner: at \qty{200}{\hertz} the numerical attenuation is $A_\mathrm{num}\!=\!0.91$ while the fitted ratio is $K_1(f)/K_1^\mathrm{low}\!=\!0.09$. A log--log regression gives $K_1(f)/K_1^\mathrm{low}\!\approx\!A_\mathrm{num}(f)^{23.6}$, so the fitted sigmoid is the linear-theory envelope raised to roughly the 24th power. No single effective permeability brings the analytical envelope into shape agreement: minimizing $\sum_f [\log A_\mathrm{analytical}(f, \mu_\mathrm{r}) - \log(K_1(f)/K_1^\mathrm{low})]^{2}$ over $\mu_\mathrm{r}$ lands at $\mu_\mathrm{r}^\mathrm{eff}\!\approx\!2.3\!\times\!10^{5}$ (15$\times$ the operating-knee value) with residuals up to $+160\,\%$ at \qty{200}{\hertz}. The fitted-sigmoid steepness is therefore intrinsic to the model class, not a calibration of linear skin attenuation.

\paragraph*{The operative mechanism: rotation-only \ac{FP} over-predicts loop widening.} Panel~(c) of Fig.~\ref{fig:maxwell_skin_depth} resolves the mechanism by direct loop-area comparison. The uniform-field rotation-only \ac{FP} at constant $K_1\!=\!K_1^\mathrm{low}$ (grey squares) over-predicts hysteretic loop area growth with frequency: from $\qty{323}{\tesla\ampere\per\meter}$ at \qty{30}{\hertz} to $\qty{929}{\tesla\ampere\per\meter}$ at \qty{200}{\hertz} ($\approx\!2.9\times$). The 1D Maxwell-coupled simulation at the same $K_1$ (blue circles) adds \qtyrange{5}{24}{\percent} widening on top of that. Maxwell coupling, therefore, pushes the model loop area in the \emph{wrong} direction relative to the fit. The uniform-field \ac{FP} at the main-text fitted $K_1(f)$ (orange triangles) lands at the measurement-matched growth ($\approx\!1.7\times$ over the same range). The fitted $K_1(f)$ reduction is therefore the rotation-only \ac{FP}'s compensation for its own over-prediction of high-frequency loop widening, not a direct shadow of skin attenuation. The corner frequency coincides with $d/\delta\!=\!1$ because both arise from the same dimensional material time scale, but different mechanisms set the magnitudes.

\paragraph*{Honest interpretation and downstream work.} The $K_1(f)$ sigmoid is the rotation-only \ac{FP}'s effective surrogate for several missing physics: (i)~the model class's intrinsic over-prediction of eddy-driven loop area growth at constant $K_1$, which is the dominant contribution, (ii)~the nonlinear $\mu_\mathrm{r}(H, B_\mathrm{peak})$ over the cycle, where the actual permeability of \ac{HGO} Fe-Si rises from $\sim\!1$ at saturation to $\sim\!10^{5}$ at the knee so that a constant operating-knee $\mu_\mathrm{r}$ understates the effective attenuation, (iii)~the increasing fraction of \ang{180}- and \ang{90}-wall switching at high $f$, which the rotation-only Fokker--Planck operator cannot represent (the Daniel-style multiscale ablation in Supplementary Table~\ref{tab:s9_unification_ablation}, row~6, separately rules out wall mechanics as the dominant residual source). A full quantitative reconciliation requires Maxwell-coupled \ac{FEM} with the present rate-independent constitutive layer applied per lamination slice and a nonlinear $\mu_\mathrm{r}(H)$ update at every time step. That treatment lies outside the scope of a constitutive modeling paper and is the natural next step for downstream device simulation work. Reproducible from \texttt{examples/04\_maxwell\_skin\_depth\_1d.py} in the open-source companion repository.

\section*{Supplementary Tables}

\begin{table}[htbp]
	\centering
	\caption{\label{tab:s1_baghel_qs}\textbf{Fitted parameters for the Baghel 27M-OH \qty{1}{\hertz} quasi-static loops.} Both directions fitted independently at frozen $\tauR\Dcoef=1$ (canonical anhysteretic limit). Values from \texttt{examples/data/baghel27M\_OH\_quasi\_static\_1hz.json}.}
	\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lcc@{}}
		\toprule
		\textbf{Parameter} & \textbf{\acs{RD}} & \textbf{\acs{TD}} \\
		\midrule
		$\Ms$ in \unit{\ampere\per\meter} & $1.628\!\times\!10^{6}$ & $1.225\!\times\!10^{6}$ \\
		$a$ in \unit{\ampere\per\meter}   & \num{47.4}              & \num{13.0} \\
		$\alpha$                        & \num{4.87e-5}           & \num{6.28e-3} \\
		$\Ku$ in \unit{\joule\per\meter\cubed}                & $564$                   & $4500$ \\
		$\Kone$ in \unit{\joule\per\meter\cubed}              & $0$ (pinned)            & $0$ (pinned) \\
		$\tauR$ in \unit{\second} (frozen, anhysteretic limit) & $10^{-3}$ & $10^{-3}$ \\
		$\Dcoef$ (frozen)               & $10^{3}$                & $10^{3}$ \\
		\midrule
		\acs{RMSE} in \unit{\tesla}             & \num{0.033}             & \num{0.070} \\
		$R^{2}$                         & $0.9994$                & $0.9945$ \\
		\bottomrule
	\end{tabular*}
\end{table}

\begin{table}[htbp]
	\centering
	\small
	\caption[Per-frequency knobs and per-loop residuals for the Hi-B 27M-OH joint $S^{2}$ fit.]{\label{tab:s3_baghel_S2}\textbf{Per-frequency knobs and per-loop residuals for the Hi-B 27M-OH joint $S^{2}$ fit.} The four shared \ac{FP}+Preisach blocks have one set of frequency-resolved values each. $K_1(f)$ is computed from the three intrinsic constants $(K_1^\mathrm{low}, f_\mathrm{c}, n)$ in main-text Table~2. The two \ac{TD} Preisach hysterons model the discrete cubic $\langle 001\rangle\!\to\!\langle 100\rangle$ and $\langle 100\rangle\!\to\!\langle 110\rangle$ switching events, the rotation-only \ac{FP} operator cannot resolve. Their amplitudes are locked to $M_\mathrm{pin}\!=\!\Ms/\sqrt{2}$, the geometric cubic projection on the lab-\ac{TD} axis~\cite{daniel2008}. Values from \texttt{examples/data/baghel27M\_OH\_unified\_fit.json}.}
	\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lcccc@{}}
		\toprule
		\textbf{Per-frequency knob}                  & \textbf{\qty{50}{\hertz}} & \textbf{\qty{100}{\hertz}} & \textbf{\qty{150}{\hertz}} & \textbf{\qty{200}{\hertz}} \\
		\midrule
		$K_1(f)$ in \unit{\joule\per\meter\cubed}      & \num{1413} & \num{798}  & \num{384}  & \num{197} \\
		$\tauR(f)$ in \unit{\milli\second}             & \num{21.1} & \num{18.4} & \num{19.4} & \num{18.3} \\
		\midrule
		\multicolumn{5}{l}{\emph{\acs{TD} Preisach hysteron 1 (cubic $\langle 001\rangle\!\to\!\langle 100\rangle$ switch)}} \\
		$H_\mathrm{pin,TD1}(f)$ in \unit{\ampere\per\meter}  & \num{193} & \num{268} & \num{253} & \num{234} \\
		$\sigma_\mathrm{pin,TD1}(f)$ in \unit{\ampere\per\meter} & \num{119} & \num{102} & \num{66}  & \num{106} \\
		$w_\mathrm{pin,TD1}(f)$                            & \num{0.47}& \num{0.41}& \num{0.25}& \num{0.41} \\
		\midrule
		\multicolumn{5}{l}{\emph{\acs{TD} Preisach hysteron 2 (cubic $\langle 100\rangle\!\to\!\langle 110\rangle$ switch)}} \\
		$H_\mathrm{pin,TD2}(f)$ in \unit{\ampere\per\meter}  & \num{239} & \num{263} & \num{282} & \num{301} \\
		$\sigma_\mathrm{pin,TD2}(f)$ in \unit{\ampere\per\meter} & \num{118} & \num{117} & \num{120} & \num{115} \\
		$w_\mathrm{pin,TD2}(f)$                            & \num{0.50}& \num{0.47}& \num{0.43}& \num{0.49} \\
		\midrule
		\acs{TD} asymmetric shift $H_\mathrm{shift}(f)$ in \unit{\ampere\per\meter} & \num{286} & \num{355} & \num{341} & \num{299} \\
		\midrule
		Per-loop \acs{RMSE} \acs{RD} in \unit{\tesla}              & \num{0.062} & \num{0.066} & \num{0.152} & \num{0.129} \\
		Per-loop \acs{RMSE} \acs{TD} in \unit{\tesla}              & \num{0.039} & \num{0.033} & \num{0.031} & \num{0.042} \\
		\bottomrule
	\end{tabular*}
\end{table}

\begin{table}[htbp]
	\centering
	\caption{\label{tab:s4_rings}\textbf{Isotropic-mode fits for all seven non-oriented ring cores measured on the in-house MagHyst$^\circledR$ Modular system}, with $\Ku=\Kone=0$ enforced. The fitted $\Ku$ and $\Kone$ remain within \qty{10}{\joule\per\meter\cubed} of zero when freed in a double-check fit, confirming that the joint-fit machinery does not invent anisotropy on isotropic samples. All seven rings pass the $\mathrm{RMSE}_\mathrm{B}/B_\mathrm{peak}<5\,\%$ acceptance gate with mean relative residual 3.6\,\% across the campaign. All seven rings are plotted in main-text Fig.~6. Per-ring fitted parameter sets for all seven materials are available as JSON files in \texttt{examples/data/fitted\_params/}.}
	\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lcccc@{}}
		\toprule
		\textbf{Material} & \textbf{$B_\mathrm{peak}$ in \unit{\tesla}} &
		\textbf{\acs{RMSE} in \unit{\tesla}} & \textbf{$R^{2}$} &
		\textbf{$\mathrm{RMSE}/B_\mathrm{peak}$} \\
		\midrule
		X10CrAlSi7 (stainless Fe-Cr-Al-Si)   & \num{1.83} & \num{0.046} & \num{0.9992} & \qty{2.5}{\percent} \\
		HLB8 (high-strength low-alloy)       & \num{2.11} & \num{0.066} & \num{0.9987} & \qty{3.1}{\percent} \\
		C35 (medium-carbon structural)       & \num{1.89} & \num{0.067} & \num{0.9984} & \qty{3.5}{\percent} \\
		S235JR (structural steel)            & \num{1.63} & \num{0.070} & \num{0.9977} & \qty{4.3}{\percent} \\
		X20Cr13 (martensitic stainless)      & \num{1.88} & \num{0.075} & \num{0.9981} & \qty{4.0}{\percent} \\
		HC260LA (cold-rolled HSLA)           & \num{2.18} & \num{0.078} & \num{0.9983} & \qty{3.6}{\percent} \\
		isovac HP1400-100A (electrical \acs{NGO})  & \num{2.02} & \num{0.087} & \num{0.9976} & \qty{4.3}{\percent} \\
		\midrule
		\multicolumn{2}{l}{\textbf{Mean across 7 rings}}
		&     \textbf{\num{0.070}}
		& & \textbf{\qty{3.6}{\percent}} \\
		\bottomrule
	\end{tabular*}
\end{table}

\begin{table}[htbp]
	\centering
	\caption{\label{tab:s5_liluo}\textbf{Joint single-parameter-set fit for the Li and Luo 2024 amplitude sweep~\cite{liluo2024}}, with the FULL-physics parameter set shared across all 20 amplitude levels at \qty{50}{\hertz} (80-trial Optuna \ac{TPE} Bayesian optimization, in-memory storage). The joint mean \ac{RMSE} exceeds \qty{0.5}{\tesla} and grows above $B_\mathrm{peak}\!\approx\!\qty{1.2}{\tesla}$. This identifies the Langevin scale parameter $a$ as amplitude-dependent in practice. The per-amplitude refit (only $a$ allowed to track $B_\mathrm{peak}$) reaches mean $\mathrm{RMSE}_\mathrm{B}\!\approx\!\qty{0.04}{\tesla}$ across the same 20 levels (see main-text amplitude-transferability discussion). Values from \texttt{examples/data/b30p105\_joint\_S2.json}.}
	\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lc@{}}
		\toprule
		\textbf{Joint Li \& Luo amplitude-sweep parameters (B30P105 \acs{HGO}, 20 levels at \qty{50}{\hertz})} & \textbf{Value} \\
		\midrule
		$\Ms$ in \unit{\ampere\per\meter}   & \num{1.551e6} \\
		$a$ in \unit{\ampere\per\meter}     & \num{2.25} \\
		$\alpha_\parallel$                & \num{1.08e-5} \\
		$\alpha_\perp$                    & \num{4.27e-4} \\
		$\Ku$ in \unit{\joule\per\meter\cubed}                  & \num{132} \\
		$\Kone$ in \unit{\joule\per\meter\cubed}                & \num{220} \\
		$\sigma_\mathrm{e}$ in \unit{\degree}        & \num{1.08} \\
		$\Dcoef$                          & \num{1915} \\
		$\tau_\mathrm{hyst}$ in \unit{\milli\second} & \num{2.65} \\
		$k_\mathrm{eddy}$                 & \num{4.1e-3} \\
		$k_\mathrm{excess}$               & \num{3.78e-2} \\
		\midrule
		mean \acs{RMSE} over 20 levels in \unit{\tesla}     & \num{0.515} \\
		worst level \acs{RMSE} (\qty{1.7}{\tesla}) in \unit{\tesla}     & \num{1.024} \\
		best level \acs{RMSE} (\qty{0.1}{\tesla}) in \unit{\tesla}      & \num{0.121} \\
		\bottomrule
	\end{tabular*}
\end{table}

\begin{table}[htbp]
	\centering
	\caption{\label{tab:s6_per_loop}\textbf{Per-loop residuals for the Baghel 27M-OH direction-resolved $S^{1}$ baseline fit} (Table~\ref{tab:s1_baghel_qs}). Reproduced here for completeness. The joint $S^{2}$ fit residuals are reported in main-text Table~4.}
	\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lcccc@{}}
		\toprule
		\textbf{Frequency} & \textbf{\acs{RD} \acs{RMSE} in \unit{\tesla}} & \textbf{\acs{RD} $R^{2}$} &
		\textbf{\acs{TD} \acs{RMSE} in \unit{\tesla}} & \textbf{\acs{TD} $R^{2}$} \\
		\midrule
		\qty{1}{\hertz}   (anhysteretic) & $0.033$ & $0.9994$ & $0.070$ & $0.9945$ \\
		\qty{50}{\hertz}                 & $0.186$ & $0.9781$ & $0.079$ & $0.9920$ \\
		\qty{100}{\hertz}                & $0.171$ & $0.9806$ & $0.039$ & $0.9983$ \\
		\qty{150}{\hertz}                & $0.173$ & $0.9814$ & $0.048$ & $0.9976$ \\
		\qty{200}{\hertz}                & $0.151$ & $0.9862$ & $0.074$ & $0.9930$ \\
		\bottomrule
	\end{tabular*}
\end{table}

\begin{table}[htbp]
	\centering
	\caption{\label{tab:s7_identifiability}\textbf{Practical identifiability summary} (diagonal of the Fisher information matrix defined in the main-text Methods subsection ``Differential evolution and parameter identifiability''), evaluated at the converged optimum of the joint $S^{2}$ Baghel fit. All standard deviations $\sigma_{\vartheta_i}/\vartheta_i \ll 1$, confirming that the seven material parameters plus per-frequency $\tauR$ are practically identifiable from the eight-loop multi-condition dataset.}
	\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lcc@{}}
		\toprule
		\textbf{Parameter} & \textbf{Mean (over multi-condition fit)}
		& \textbf{Std.\ dev.\ $\sigma_{\vartheta_i}$} \\
		\midrule
		$\Ms$       & $1.411\!\times\!10^{6}$ & $4.7\!\times\!10^{4}$ ($3.3\,\%$) \\
		$a$         & $0.75$                  & $0.04$ ($5.3\,\%$) \\
		$\alpha$    & $7.28\!\times\!10^{-5}$ & $3.4\!\times\!10^{-6}$ ($4.7\,\%$) \\
		$\Ku$       & $349$                   & $14$ ($4.0\,\%$) \\
		$\Kone^\mathrm{low}$ & $1619$                  & $83$ ($5.1\,\%$) \\
		$\Dcoef$    & $4892$                  & $192$ ($3.9\,\%$) \\
		$\tauR(f)$  & \qtyrange{20.2}{24.9}{\milli\second} & \qtyrange{0.5}{1.0}{\milli\second} ($\approx\!4\,\%$) \\
		\bottomrule
	\end{tabular*}
\end{table}

\begin{table}[htbp]
	\centering
	\small
	\caption{\label{tab:s8_models_per_ring}\textbf{Per-ring cross-model comparison.} \ac{RMSE} in \unit{\tesla} of each baseline model fitted independently to each of the seven ring-core measurements, versus the present model fitted in its rate-independent isotropic reduction ($\Ku=\Kone=0$, $\Dcoef\to 0$). Each baseline is given its own optimal parameter set on each ring (maximally favorable to the baselines). The present model uses the corresponding row of Table~\ref{tab:s4_rings}. On single-loop isotropic samples, the seven approaches all perform within measurement noise. The per-loop ranking is therefore not the discriminator the present work addresses, which is multi-condition transferability of a single shared parameter set (main-text Table~3 and Fig.~\ref{fig:s_models_hgo}). Numerical entries are reproduced from the per-ring baseline campaign in \texttt{examples/data/baghel27M\_OH\_per\_freq\_baselines.json}. Mean residuals across the seven rings are listed in the main-text Table~3 (\acs{NGO} ring rows).}
	\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lccccccc@{}}
		\toprule
		\textbf{Material (Ring)} & \textbf{J--A} & \textbf{Chwastek} &
		\textbf{\acs{DDD}} & \textbf{\acs{DTHMS}} &
		\textbf{Preisach} & \textbf{\acs{PI}} &
		\textbf{This work} \\
		\midrule
		X10CrAlSi7         & \num{0.023} & \num{0.008} & \num{0.006} & \num{0.024} & \num{0.126} & \num{0.160} & \num{0.046} \\
		HLB8               & \num{0.041} & \num{0.041} & \num{0.041} & \num{0.072} & \num{0.166} & \num{0.185} & \num{0.066} \\
		C35                & \num{0.045} & \num{0.046} & \num{0.055} & \num{0.054} & \num{0.122} & \num{0.170} & \num{0.067} \\
		S235JR             & \num{0.028} & \num{0.029} & \num{0.023} & \num{0.020} & \num{0.111} & \num{0.208} & \num{0.070} \\
		X20Cr13            & \num{0.032} & \num{0.047} & \num{0.040} & \num{0.042} & \num{0.121} & \num{0.186} & \num{0.075} \\
		HC260LA            & \num{0.054} & \num{0.050} & \num{0.057} & \num{0.093} & \num{0.176} & \num{0.195} & \num{0.078} \\
		Isovac HP1400-100A & \num{0.047} & \num{0.056} & \num{0.054} & \num{0.067} & \num{0.155} & \num{0.196} & \num{0.087} \\
		\midrule
		\textbf{Mean (7 rings)}              & \textbf{0.039} & \textbf{0.040} & \textbf{0.039} & \textbf{0.053} & \textbf{0.140} & \textbf{0.186} & \textbf{0.070} \\
		\bottomrule
	\end{tabular*}
\end{table}

\begin{table}[htbp]
	\centering
	\small
	\caption[Unification ablation on Hi-B 27M-OH.]{\label{tab:s9_unification_ablation}\textbf{Unification ablation on Hi-B 27M-OH.} Per-loop $\mathrm{RMSE}_\mathrm{B}$ in \unit{\tesla} for six configurations of the vector Fokker--Planck model. The skin-depth $K_1(f)\!=\!K_1^\mathrm{low}/(1\!+\!(f/f_\mathrm{c})^n)$ phenomenology (this work, row~1) reduces the joint mean from \qty{0.205}{\tesla} (static $K_1$, row~2) to \qty{0.069}{\tesla} by closing the high-frequency \ac{RD} residual (\num{0.385}\,$\to$\,\qty{0.152}{\tesla} at \qty{150}{\hertz} and \num{0.368}\,$\to$\,\qty{0.129}{\tesla} at \qty{200}{\hertz}) and the \ac{TD} inner-knee structure (\num{0.172}\,$\to$\,\qty{0.036}{\tesla} direction mean). Row~3 is the per-direction $S^{2}$ fit with the same architecture (all parameters free per direction, no unification constraint). It reaches mean \ac{RMSE} \qty{0.063}{\tesla} versus \qty{0.069}{\tesla} for this work, so the cross-direction unification costs $\approx\!10\,\%$ of mean \ac{RMSE} relative to the per-direction floor. The high-frequency \ac{RD} residuals (\num{0.138} and \qty{0.108}{\tesla} at 150 and \qty{200}{\hertz}) are nearly the same in row~3 as in row~1, confirming that they reflect a rotation-only constitutive limit and not a unification penalty. Row~4: strict $S^{2}\!\to\!S^{1}$ reduction with $K_1\!=\!0$ pinned degrades both directions, showing $K_1$ in the \ac{FP} equilibrium is not decorative. Row~5: \acs{DTHMS} baseline (reproduced from main-text Table~3). Row~6: Daniel six-domain wall-fraction extension~\cite{daniel2008} (implemented as \texttt{hysterpy/daniel\_multiscale.py}) replaces the orientation density on $S^{2}$ with six thermally-activated cubic-axis populations $f_i(t)$. This admits literature single-crystal $\Kone\!=\!\qty{4.8e4}{\joule\per\meter\cubed}$ by construction but reaches joint \ac{RMSE} \qty{0.50}{\tesla} ($\approx\!3\times$ the rotation-only \ac{FP} residual of \qty{0.17}{\tesla}), empirically ruling out missing domain-wall mechanics as the dominant residual source and isolating the spatial-field (skin-depth) assumption as the operative limitation at $f\!\ge\!\qty{150}{\hertz}$. Two further variants reported in the text but not tabulated here: free Bertotti back-field $\rho_\mathrm{eddy}, \rho_\mathrm{excess}$ drive to zero and recover the static-$K_1$ minimum bit-for-bit. The per-direction static-$K_1$ fit reaches \qty{0.201}{\tesla} versus \qty{0.205}{\tesla} for this work with static $K_1$. Visualisation in Fig.~\ref{fig:s_unification}. Parameters from this work in \texttt{examples/data/baghel27M\_OH\_unified\_fit.json}.}
	\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}lcccc cccc c@{}}
		\toprule
		& \multicolumn{4}{c}{Rolling direction (\acs{RD})}  & \multicolumn{4}{c}{Transverse direction (\acs{TD})} &        \\
		\cmidrule(lr){2-5}\cmidrule(lr){6-9}
		Fit variant                        &   50  &  100  &  150  &  200  &   50  &  100  &  150  &  200  & joint   \\
		\midrule
		\textbf{This work ($K_1(f)$, Fig.~3)}              & \textbf{0.062} & \textbf{0.066} & \textbf{0.152} & \textbf{0.129} & \textbf{0.039} & \textbf{0.033} & \textbf{0.031} & \textbf{0.042} & \textbf{0.069} \\
		Same $S^{2}$ base, static $K_1$              & 0.076 & 0.123 & 0.385 & 0.368 & 0.205 & 0.140 & 0.131 & 0.213 & 0.205 \\
		Per-direction $S^{2}$, $K_1$ free (perbasin) & 0.068 & 0.068 & 0.138 & 0.108 & 0.032 & 0.029 & 0.028 & 0.031 & 0.063 \\
		Per-direction $S^{2}$, $K_1\!=\!0$ pinned    & 0.303 & 0.324 & 0.230 & 0.248 & 0.228 & 0.177 & 0.160 & 0.197 & 0.233 \\
		Per-direction $S^{1}$ \acs{DTHMS} (Table~3 col)  & 0.070 & 0.073 & 0.059 & 0.068 & 0.083 & 0.080 & 0.073 & 0.097 & 0.075 \\
		Daniel 6-domain wall-fraction~\cite{daniel2008}     & 0.64  & 0.73  & 0.72  & 0.81  & 0.25  & 0.25  & 0.27  & 0.34  & 0.50 \\
		\bottomrule
	\end{tabular*}
\end{table}

\begin{table}[htbp]
	\centering
	\caption[Cross-validation summary for the joint $S^{2}$ fit.]{\label{tab:cv_summary}\textbf{Cross-validation summary for the joint $S^{2}$ fit on Hi-B 27M-OH.} All variants share the same $S^{2}$ Fokker--Planck simulator and the same 8-loop measurement set. Only the joint-cost composition (which frequencies are held in versus held out) or the $K_1(f)$ functional form changes. The forward-extrapolation row fits \qty{50}{\hertz} only and predicts the higher frequencies via the scaling laws of the main text. The reverse-extrapolation rows hold the shared base fixed at this work and refit only the four $K_1(f)$ constants on the \qtyrange{100}{200}{\hertz} loops. The leave-one-out row is the 4-fold version of the same reverse-extrapolation protocol. The free-parameter count for the all-8-loop joint fit (38) comprises the 22 shared parameters of Methods plus the per-frequency \ac{TD} Preisach-overlay arrays. The forward-extrapolation row fits a single frequency, so its overlay is not expanded per frequency and the count stays at 22.}
	\begin{tabular*}{\linewidth}{@{\extracolsep{\fill}}llcc@{}}
		\toprule
		\textbf{Variant} & \textbf{Held-in / functional form} & \textbf{Fit cost in \unit{\tesla}} & \textbf{Free params} \\
		\midrule
		This work (joint fit)       & all 8 loops, $\Ms$ free, sigmoid $K_1(f)$  & \num{0.069} & 38 \\
		Forward-extrapolation       & fit \qty{50}{\hertz}, predict others        & $0.043^\dagger$ & 22 \\
		Reverse-extrapolation       & fit \qtyrange{100}{200}{\hertz}, predict \qty{50}{\hertz}, sigmoid & $0.072$ & 4 \\
		Reverse-extrap.\ (Hill)     & same split, Hill function rising           & $0.092$ & 4 \\
		Leave-one-out (4 folds)     & sigmoid, held-out frequency rotated        & $0.062$--$0.077$ & 4 \\
		\bottomrule
		\multicolumn{4}{l}{\footnotesize $^\dagger$ \qty{50}{\hertz}-only fit cost. The other three frequencies are predicted via scaling laws.}
	\end{tabular*}
\end{table}

\section*{Supplementary Figures}

\begin{figure}[htbp]
	\centering
	\includegraphics[width=\linewidth]{figures/lab_crystal_geometry.pdf} \caption{\label{fig:s_geometry}\textbf{Lab-frame versus crystal-frame geometry for the Hi-B 27M-OH Goss-textured sample.} Left: the cubic crystallographic axes $[100],\,[010],\,[001]$ of a representative grain. The Goss texture orients $[001]$ along the lab rolling direction with a residual misorientation spread $\sigma_\mathrm{e}\!=\!\qty{0.053}{\radian}$ (red cone). Right: the lab-frame field directions overlaid on the same crystal axes. The rolling-direction field $\mathbf{H}_\mathrm{RD}$ is co-aligned with the cubic easy axis $[001]$ (blue arrow). The transverse-direction field $\mathbf{H}_\mathrm{TD}$ is aligned with the cubic $[110]/\sqrt{2}$ direction (green arrow), at \ang{90} to the easy axis and \ang{45} to each of the two perpendicular cubic minima $[100]$ and $[010]$. The asymmetry is the physical reason why the unified parameter set cannot simultaneously saturate \ac{RD} and \ac{TD} at the same applied field amplitude: along \ac{RD} the bulk magnetisation saturates at $\Ms\,\langle\cos\theta_\mathrm{g}\rangle_{\mathrm{g}}\!\approx\!0.99\,\Ms$ (with $\langle\cdot\rangle_{\mathrm{g}}$ the grain-orientation average of the main text) for the fitted Goss spread, requiring a sharply-knee'd $M$--$H$ characteristic to fully reach this cap at the measured $H_\mathrm{peak}\!\approx\!\qty{250}{\ampere\per\meter}$. Along \ac{TD}, the magnetization rotates through the $\langle 100\rangle$ and $\langle 010\rangle$ cubic wells before aligning with the applied field, thereby lowering the effective saturation field and matching the measurement at a lower $H_\mathrm{peak}$. The same scalar $\Ms$ produces a small \ac{RD} undershoot at peak field and an essentially exact \ac{TD} match (main-text Fig.~3), the latter dominating the joint cost weighting.}
\end{figure}

\begin{figure}[htbp]
	\centering
	\includegraphics[width=\linewidth]{figures/baghel_unification_comparison.pdf} \caption{\label{fig:s_unification}\textbf{Visual proof of the unification-cost ablation in Supplementary Table~\ref{tab:s9_unification_ablation}.} For each measurement loop, this work (black dashed, identical to main-text Fig.~3), the per-direction $S^{2}$ fit with $K_1$ free (orange dotted), and the per-direction $S^{2}$ fit with $K_1\!=\!0$ pinned (blue dash-dot) are overlaid on the measurement (grey solid). Per-panel \ac{RMSE} is annotated in the upper-left of each axis. This work and the per-direction-free curves are visually indistinguishable on \ac{RD} across all frequencies. This work is already at the direction-resolved optimum. The $K_1\!=\!0$ trace widens visibly on \ac{RD} 150 and \qty{200}{\hertz}, closing roughly $30\,\%$ of the high-frequency \ac{RD} residual, but at the cost of losing the low-frequency \ac{RD} knee shape and the \ac{TD} inner-knee structure entirely. The Bertotti additive back-field configuration (free $\rho_\mathrm{eddy}, \rho_\mathrm{excess}$) is not plotted as a separate trace because the optimizer drives both constants to zero and recovers exactly the static-$K_1$ curve.}
\end{figure}

\begin{figure}[htbp]
	\centering
	\includegraphics[width=\linewidth]{figures/baghel_quasi_static_1hz_overlay.pdf} \caption{\label{fig:s_qs}\textbf{\qty{1}{\hertz} quasi-static fit on the Baghel 27M-OH \ac{RD} and \ac{TD} loops}, with $\tauR$ frozen at its anhysteretic value (Supplementary Note~S4: $\omega\tauR\!\ll\!1$ at \qty{1}{\hertz}). $(\Ms, a, \alpha, \Ku)$ fitted independently per direction and reported in Table~\ref{tab:s1_baghel_qs}. Symbols: measurement (black circles), model (orange squares). $R^{2}=0.9994$ for \ac{RD} and $0.9945$ for \ac{TD}.}
\end{figure}

\begin{figure}[p]
	\centering
	\includegraphics[width=\linewidth]{figures/b30p105_per_model_full.pdf}
	\caption[Cross-model overlay on Li and Luo B30P105 at four representative amplitudes.]{\label{fig:s_models_hgo}\textbf{Cross-model overlay on the Li and Luo B30P105 high-permeability grain-oriented steel} at $B_\mathrm{peak}$ of \qtylist{0.5;1.0;1.5;2.0}{\tesla} (four columns), \qty{50}{\hertz}. Each row is one of the six baselines defined in main-text Theory: Jiles--Atherton~\cite{jiles1986}, Chwastek modified J--A~\cite{chwastek2009}, D'Aloia--Di~Francesco--De~Santis~\cite{daloia2021}, \acs{DTHMS} (this work), Preisach~\cite{mayergoyz}, and Prandtl--Ishlinskii~\cite{krasnoselskii}. Each baseline is fitted \emph{independently} to each amplitude with its own optimal parameter set (maximally favorable to the baselines). The present model uses a single shared parameter set that is transferred across all 20 amplitudes and is therefore not retuned for each amplitude. The per-amplitude \ac{RMSE} annotated in each panel is the value contributing to the means in main-text Table~3. Visual inspection makes the amplitude-transferability narrative concrete: the per-loop fit of every J--A-family baseline degrades smoothly with $B_\mathrm{peak}$, while Preisach and \acs{PI} exhibit their characteristic saturation artifacts at $B_\mathrm{peak}\!>\!\qty{1.5}{\tesla}$.}
\end{figure}

\begin{figure}[htbp]
	\centering
	\includegraphics[width=\linewidth]{figures/ring1303_cross_model.pdf}
	\caption[Cross-model overlay on the X10CrAlSi7 ring-core.]{\label{fig:s_models_ring}\textbf{Cross-model overlay on the X10CrAlSi7 ring-core isotropic baseline} at $B_\mathrm{peak}\!=\!\qty{1.83}{\tesla}$, $\approx\!\qty{1.7}{\hertz}$ quasi-static MagHyst$^\circledR$ excitation. Panel (a): all six J--A-family and operator baselines plus the present model overlaid on the measured loop (grey). Panel (b): per-model \ac{RMSE} bar chart. On a single isotropic loop, the J--A family, the \acs{DTHMS}, and the present model (rate-independent reduction with $\Ku\!=\!\Kone\!=\!0$, $\Dcoef\!\to\!0$) all perform within measurement noise. Preisach and \acs{PI} show the characteristic constant-amplitude saturation artifact. The cross-model differentiation reported in main-text Table~3 arises only with multi-condition data (Baghel \ac{RD}/\ac{TD}, Li and Luo amplitude sweep), where baseline retuning per loop is required. In contrast, the present model uses a single shared set. Reproducible end-to-end from the open-source bundled ring CSV via \texttt{examples/01\_fit\_all\_models\_demo.py}.}
\end{figure}

\begin{figure}[htbp]
	\centering
	\includegraphics[width=\linewidth]{figures/baghel_bimodal_basins.pdf} \caption{\label{fig:s_bimodal_basins}\textbf{Bimodal cost surface in the joint S² Baghel fit.} Each panel overlays the measurement (grey) with two distinct local minima of the joint arc-length-\ac{RMSE} cost over the 17-dimensional extension of the model (anisotropic Bingham \ac{ODF}, per-frequency $\tauR$, per-direction $\Ku$, amplitude-coupled $a$, and $K_{2}$ second cubic invariant). Orange dashed = warm-start \ac{FP}-only (this work, main-text Fig.~3. Mean $\varepsilon_\mathrm{A} = \qty{0.22}{\tesla}$ balanced across \ac{RD} and \ac{TD}). Red dotted = cold-start basin B discovered by an independent \ac{TPE} search seeded with 80 random Sobol points. $\sigma_\parallel = \ang{12.6}$, $\sigma_\perp = \ang{1.5}$ (highly anisotropic Goss spread), $\Kone = \qty{2.4e4}{\joule\per\meter\cubed}$, $\Ku^\mathrm{TD}/\Ku^\mathrm{RD} = 3.8$. Mean $\varepsilon_\mathrm{B} = \qty{0.55}{\tesla}$ (\ac{RD} $=\qty{0.88}{\tesla}$, \ac{TD} $=\qty{0.22}{\tesla}$). Basin B improves \ac{TD} by \qty{13}{\percent} but breaks \ac{RD} by a factor of 4. A third search (basin~C, not plotted here for clarity) constrained the Bingham distribution to be physically correct ($\sigma_\parallel \le \sigma_\perp$) and forced $\Kone \in$ \qtyrange{1e4}{5e4}{\joule\per\meter\cubed} within the Fe-Si single-crystal literature range. The resulting fit selects $\sigma_\parallel\!=\!\ang{2.9}$, $\sigma_\perp\!=\!\ang{18.3}$, $\Kone\!=\!\qty{2.96e4}{\joule\per\meter\cubed}$ but with $\varepsilon_\mathrm{C}\!=\!\qty{1.11}{\tesla}$ (\ac{RD}~$=\!\qty{1.37}{\tesla}$, \ac{TD}~$=\!\qty{0.84}{\tesla}$): the cubic [100]/[010] easy wells perpendicular to the lab-\ac{RD} field trap magnetisation at low \ac{RD} excitation, shifting the \ac{RD} coercivity by $\approx\!\qty{0.3}{\kilo\ampere\per\meter}$ and shrinking apparent saturation. The three basins together establish that the one-shared-parameter-set model class possesses a unique global minimum (the \ac{FP}-only fit). That physically correct Bingham parameters lie on a higher-cost ridge. The bimodality is therefore a structural property of the model class rather than a search artifact, thereby motivating the natural extensions discussed in the main-text Discussion (two-component \ac{ODF} domain-wall pinning hysteron stack).}
\end{figure}

\begin{figure}[htbp]
	\centering
	\includegraphics[width=\linewidth]{figures/cross_model_rmse.pdf} \caption{\label{fig:s_cross_model_bars}\textbf{Cross-model \ac{RMSE} bar chart across the seven \ac{NGO} ring-core measurements.} Panel~(a) groups the seven baseline models per ring. Panel~(b) shows the mean per model, averaged over rings. The narrow spread visible in (b) confirms the per-loop \ac{NGO} observation made in Fig.~\ref{fig:s_models_ring}: on isotropic single loops, the seven approaches all perform within measurement noise, with the present rate-independent reduction (\acs{DTHM} / ``This work'', last bar) comparable to the Jiles family. The cross-condition spread on Baghel~\ac{RD}/\ac{TD} and Li and Luo amplitudes is reported in the main-text Table~3. This figure addresses only the ring-core column. Numerical entries reproduced from \texttt{examples/data/baghel27M\_OH\_per\_freq\_baselines.json}.}
\end{figure}

\begin{figure}[htbp]
	\centering
	\includegraphics[width=\linewidth]{figures/multistart_basins.pdf}
	\caption[Twenty-seed cold-start multistart on the joint $S^{2}$ cost surface.]{\label{fig:multistart_basins}\textbf{Twenty-seed cold-start multistart on the joint $S^{2}$ cost surface for Hi-B 27M-OH.} (a)~Parametric scatter of the best-of-seed sigmoid constants $(K_1^{\mathrm{low}}, f_\mathrm{c})$ on log--log axes. Basin assignments are color-coded, and the main joint fit is marked with a star. (b)~Resulting $K_1(f) = K_1^{\mathrm{low}}/(1+(f/f_\mathrm{c})^n)$ trajectories for each seed in the range \qtyrange{20}{400}{\hertz}, with the main-text joint-fit trajectory overlaid as a black dashed line. Basin~A (blue, $N\!=\!14$) brackets the main-fit curve tightly. Basin~B (red, $N\!=\!2$) drops sharply at $f\!>\!\qty{100}{\hertz}$ from a single-crystal-magnitude initial value. Basin~B$'$ (orange, $N\!=\!4$) maintains a low, near-constant $K_1$ across the measurement band. Grey dotted lines mark the four measurement frequencies \qtylist{50;100;150;200}{\hertz}.}
\end{figure}

\begin{figure}[htbp]
	\centering
	\includegraphics[width=\linewidth]{figures/maxwell_skin_depth.pdf}
	\caption[1D Maxwell-coupled analysis of the $K_1(f)$ corner frequency.]{\label{fig:maxwell_skin_depth}\textbf{1D Maxwell-coupled analysis of the cubic-anisotropy sigmoid corner frequency for Hi-B 27M-OH.} (a)~Steady-state internal field profile $H_\mathrm{internal}(z)/H_\mathrm{lab,peak}$ at the surface-field peak, across the lamination cross-section $z/d\!\in\![-0.5,+0.5]$, for two representative frequencies. At \qty{50}{\hertz} ($d/\delta\!=\!0.75$, blue circles), the interior field is nearly uniform. At \qty{200}{\hertz} ($d/\delta\!=\!1.5$, red circles) the centre-line value drops to ${\approx}0.84\,H_\mathrm{lab,peak}$, the signature of incipient skin shielding. Markers from the 1D Maxwell-coupled simulation. Lines connect grid points. (b)~Frequency dependence of the bulk-averaged attenuation factor $A(f) = \max_t|\langle H(z,t)\rangle_z|/\max_t|H_\mathrm{lab}(t)|$. Blue circles: numerical solution of the coupled magnetic-diffusion + Fokker--Planck problem with the same $\Kone\!=\!K_1^\mathrm{low}\!=\!\qty{1619}{\joule\per\meter\cubed}$ in every $z$-slice. Grey dashed: analytical linear-theory envelope $|\tanh(\gamma d/2)/(\gamma d/2)|$ at the operating-knee permeability $\mu_\mathrm{r}\!=\!1.35\!\times\!10^{4}$. Orange solid: fitted phenomenological $K_1(f)/K_1^\mathrm{low}$ sigmoid from this work in the main text. The vertical dotted line marks the fitted corner $f_\mathrm{c}\!=\!\qty{99.0}{\hertz}$. The numerical Maxwell solve matches the analytical envelope to better than $1\,\%$ and places its half-attenuation point in the same decade as the fitted corner. Both share the same dimensional scale $\mu_\mathrm{r}\sigma d^{2}\omega$, so the corner frequency of the $K_1$ sigmoid is locked by Maxwell physics. The fitted sigmoid decays much more steeply than $A(f)$ in depth (see panel~(c) for the mechanism). (c)~Per-frequency hysteresis loop area $\oint H\,\mathrm dB$ in three configurations: uniform-field rotation-only Fokker--Planck at constant $\Kone\!=\!K_1^\mathrm{low}$ (grey squares), 1D Maxwell-coupled \ac{FP} at the same constant $\Kone\!=\!K_1^\mathrm{low}$ (blue circles), uniform-field \ac{FP} at the main-text fitted $K_1(f)$ (orange triangles). The Maxwell-coupled simulation broadens the uniform-field constant-$K_1$ result, in the wrong direction relative to the fit. The fitted $K_1(f)$ reduces the model loop area to match measurement-driven growth ($\approx\!1.7\times$ over \qtyrange{30}{200}{\hertz} vs $\approx\!2.9\times$ at constant $K_1$). The $K_1(f)$ sigmoid is therefore the rotation-only \ac{FP}'s compensation for its intrinsic over-prediction of high-frequency loop widening (Supplementary Note S6). Reproducible from \texttt{examples/04\_maxwell\_skin\_depth\_1d.py}.}
\end{figure}

\begin{figure}[htbp]
	\centering
	\includegraphics[width=\linewidth]{figures/angle_resolved_validation.pdf}
	\caption[Code-equivalence consistency check against the Ramesh--Jiles--Bi 1997 anisotropic anhysteretic.]{\label{fig:angle_resolved}\textbf{Code-equivalence consistency check against the anhysteretic limit of the Jiles--Atherton anisotropic model (Ramesh, Jiles \& Bi 1997, Eq.~5)~\cite{ramesh1997}.} This is a consistency check, not a fit: both reference and present-work traces consume the \emph{same} literature parameter set $K_1\!=\!\qty{4.8e4}{\joule\per\meter\cubed}$, $\Ms\!=\!\qty{1e6}{\ampere\per\meter}$, $a\!=\!\qty{1000}{\ampere\per\meter}$, $\alpha\!=\!\num{e-3}$, and they must agree if the $S^{2}$ Fokker--Planck operator correctly reduces to the J--A anisotropic anhysteretic in the $D\tau\!=\!1$ canonical equilibrium limit. Solid lines: anhysteretic initial-magnetization curves along the three crystallographic directions $[100]$, $[110]$, $[111]$, regenerated from Ramesh's published Eq.~5 self-consistent fixed-point form. Dashed lines with markers: this work, computed directly from the Boltzmann equilibrium of the $S^{2}$ free-energy landscape with self-consistent mean-field iteration. Only the lab-frame field axis $\hat{\mathbf n}_\mathrm{f}\!\in\!\{[100],[110],[111]\}$ rotates between traces. At $H\!=\!\qty{80}{\kilo\ampere\per\meter}$ the two formulations agree within \qty{1}{\percent} for all three crystal-field directions, confirming the operator-reduction step.\@ Fitting validation against measured data is reported separately for the Hi-B 27M-OH \ac{RD}/\ac{TD} multi-frequency set (main-text Fig.~3, joint $\mathrm{RMSE}_\mathrm{B}\!=\!\qty{0.069}{\tesla}$), the Li and Luo B30P105 amplitude sweep (main-text Fig.~5), and the seven \ac{NGO} ring cores (main-text Fig.~6). Reproducible from \texttt{scripts/plot\_angle\_resolved\_validation.py}.}
\end{figure}

\begin{figure}[htbp]
	\centering
	\includegraphics[width=\linewidth]{figures/harmonic_excitation.pdf}
	\caption[Dynamic vs.\ quasi-static response to third-harmonic-distorted excitation on X10CrAlSi7.]{\label{fig:s_harmonic}\textbf{Dynamic versus quasi-static response of the rotation-only Fokker--Planck constitutive law under multi-harmonic excitation.} Practical magnetic loads are rarely pure sinusoids. Switched-mode drives, transformer over-excitation, and saturating inductors all produce odd-harmonic distortion of the field waveform. The X10CrAlSi7 ring core (isotropic fit, Table~\ref{tab:s4_rings}) is driven by $H(t)\!=\!H_{1}\sin(\omega t)+H_{2}\sin(3\omega t)$ with $H_{1}\!=\!\qty{8000}{\ampere\per\meter}$, $H_{2}\!=\!\qty{4800}{\ampere\per\meter}$ ($60\,\%$ third-harmonic content), $\omega\!=\!2\pi f$. \textbf{(a)} The non-monotone field crosses zero twice per fundamental half-cycle and produces four local field extrema per period at $H\!\approx\!\pm\qty{9.4}{\kilo\ampere\per\meter}$ with intermediate ``valleys'' at $H\!\approx\!\pm\qty{3.2}{\kilo\ampere\per\meter}$ ($\qty{50}{\hertz}$ shown over two fundamental periods). \textbf{(b)} Flux-density response: $B_\mathrm{peak}\!=\!\qty{1.421}{\tesla}$, with the third-harmonic ``valley'' producing two flat-top plateaus per half-cycle. \textbf{(c)} $B$--$H$ trajectory at $\qty{50}{\hertz}$ (dynamic, blue solid, $\omega\tau_\mathrm{R}\!=\!0.14$ on the fundamental, $0.42$ on the third harmonic) overlaid on the same excitation at $\qty{1}{\hertz}$ (quasi-static, black dashed, $\omega\tau_\mathrm{R}\!=\!2.8\times\!10^{-3}$, fully rate-effect-free). Both trajectories reach the same $B_\mathrm{peak}\!\approx\!\qty{1.42}{\tesla}$ at the operating point: the rate-dependence does \emph{not} alter the loop tips, it widens the loop body. The QS trace appears as a near-singular curve because for this soft non-oriented material the rate-independent backbone is intrinsically narrow ($\Hc^{\mathrm{QS}}\!\approx\!\qty{25}{\ampere\per\meter}$, $W_\mathrm{h}^{\mathrm{QS,pure}}\!\approx\!\qty{0.13}{\kilo\joule\per\meter\cubed}$ at the same $H_\mathrm{peak}$ with a clean pure-sinusoidal excitation). The inner reversals at $H\!\approx\!\pm\qty{3.2}{\kilo\ampere\per\meter}$ fall \emph{outside} the metastability region, so the QS trace retraces the equilibrium curve through them without forming minor loops. The two minor-loop lobes near $(H,B)\!\approx\!(\pm\qty{3}{\kilo\ampere\per\meter},\pm\qty{1}{\tesla})$ visible at \qty{50}{\hertz} are therefore entirely rate-dependent path-memory artifacts, not present in the rate-independent limit. The closed-loop area is the hysteresis loss density per cycle $W_\mathrm{h}\!=\!\oint H\,\mathrm{d}B$. The dynamic case dissipates $W_\mathrm{h}^{\mathrm{dyn}}\!=\!\qty{13.6}{\kilo\joule\per\meter\cubed}$, the quasi-static reference only $W_\mathrm{h}^{\mathrm{QS}}\!=\!\qty{0.32}{\kilo\joule\per\meter\cubed}$, a $42\times$ rate factor, dominated by the third-harmonic component which sits in a regime ($\omega\tau_\mathrm{R}\!\approx\!0.4$) where the Fokker--Planck relaxation cannot follow the field. For reference, the pure-sinusoid $H_{1}\sin(\omega t)$ baseline at \qty{50}{\hertz} (not plotted) gives $W_\mathrm{h}\!=\!\qty{5.0}{\kilo\joule\per\meter\cubed}$, so $60\,\%$ third-harmonic content inflates the dynamic loss by $2.7\times$. No re-fitting was performed: the parameter set is identical to Table~\ref{tab:s4_rings} row X10CrAlSi7, fitted only against the in-house pure-sinusoidal MagHyst$^\circledR$ measurement at \qty{1.7}{\hertz}. Reproducible from \texttt{scripts/plot\_harmonic\_excitation.py}.}
\end{figure}

\begin{figure}[htbp]
\centering
\includegraphics[width=\linewidth]{figures/hefmag_external.pdf}
\caption[External cross-validation against the HEFMAG European-metrology dataset (NO Fe-Si).]{\label{fig:s_hefmag}\textbf{External cross-validation against the open EMPIR~19ENG06 HEFMAG dataset~\cite{hefmag_zenodo} (INRiM, PTB, NPL, CMI).} Two non-oriented Fe-Si \qty{0.2}{\milli\meter} loops from the HEFMAG release are fit by the present DTHM-$S^{2}$ model in its isotropic reduction ($\Ku\!=\!\Kone\!=\!0$), from Epstein 700-turn IEC~60404-2 measurements at \qty{50}{\hertz}. \textbf{NO~L1}: $H_\mathrm{peak}\!=\!\qty{59.1}{\ampere\per\meter}$, $B_\mathrm{peak}\!=\!\qty{0.50}{\tesla}$, $\mathrm{RMSE}_\mathrm{B}\!=\!\qty{61}{\milli\tesla}$. \textbf{NO~L2}: $H_\mathrm{peak}\!=\!\qty{109.8}{\ampere\per\meter}$, $B_\mathrm{peak}\!=\!\qty{1.00}{\tesla}$, $\mathrm{RMSE}_\mathrm{B}\!=\!\qty{133}{\milli\tesla}$. Both fits are produced \emph{blind}: no prior knowledge of the HEFMAG sample composition or processing is used. An automatic Conde Garrido (2025) initialiser~\cite{condegarrido2025} extracts seed parameters $(M_\mathrm{s}, a, \alpha)$ from a Langevin--Weiss fit to the horizontal-average of the measured loop and the J--A pinning coefficient $k$ from a closed-form expression at the coercive point ($\Hc\!=\!\qty{29.7}{\ampere\per\meter}$ for NO~L1, \qty{41.8}{\ampere\per\meter} for NO~L2). The numba-parallel DE then refines a warm-started population around the seed with $\pm 50\,\%$ bounds. The blind seed alone already yields \ac{RMSE} of \qty{169}{\milli\tesla} (NO~L1) and \qty{330}{\milli\tesla} (NO~L2). DE refinement reduces these by $2.5\!-\!3\times$ to the values reported. The fitted $M_\mathrm{s}^{\mathrm{NO\,L1}}\!=\!\qty{7.3e5}{\ampere\per\meter}$ and $M_\mathrm{s}^{\mathrm{NO\,L2}}\!=\!\qty{1.07e6}{\ampere\per\meter}$ depend on loop amplitude, consistent with the per-amplitude effective Langevin scale documented in main-text Fig.~5 on the Li and Luo B30P105 dataset.
\textbf{Known limitations not plotted}: (i) NO~L3 ($H_\mathrm{peak}\!=\!\qty{1546}{\ampere\per\meter}$, $B_\mathrm{peak}\!=\!\qty{1.51}{\tesla}$): isotropic Langevin over-rounds the fat saturated loop, $\mathrm{RMSE}_\mathrm{B}\!\approx\!\qty{500}{\milli\tesla}$ persists after DE refinement. A polycrystal extension would relax this, but it is structurally inappropriate for an isotropic NO sample. (ii) HEFMAG \ac{GO} Hi-B loops at $H_\mathrm{peak}$ of \qtylist{25.6;46.7;120.9}{\ampere\per\meter}: we ran a joint 3-amplitude textured-polycrystal fit warm-started from the Baghel 27M-OH unified parameters (main Fig.~3), with wide bounds and a 80-generation multi-fidelity DE (30 coarse generations at $N_\theta\!=\!10, N_\phi\!=\!20$ followed by 50 fine generations at $N_\theta\!=\!16, N_\phi\!=\!32$). The fit converges to a joint mean $\mathrm{RMSE}_\mathrm{B}\!=\!\qty{791}{\milli\tesla}$ with per-loop residuals \qtylist{594;703;1076}{\milli\tesla}, plateauing at the lower bound of $M_\mathrm{s}\!=\!\qty{1.30e6}{\ampere\per\meter}$ with fitted $K_1\!=\!\qty{19.5}{\joule\per\meter\cubed}$ (versus Baghel's fitted \qty{1619}{\joule\per\meter\cubed}) and $\sigma_\mathrm{e}\!=\!\qty{0.026}{\radian}$ (versus Baghel's \qty{0.053}{\radian}). This plateau identifies HEFMAG \ac{GO} as requiring two further components that are already implemented in the model but not exercised in the search above: a Preisach hysteron overlay (cf.~the \ac{TD} overlay in main Fig.~3) to capture the discrete \qty{90}{\degree}-wall switching of the very-square loop ($\Br/B_\mathrm{sat} > 0.9$ at $H_\mathrm{peak}\!=\!\qty{26}{\ampere\per\meter}$), and an anisotropic Bingham orientation distribution with $\sigma_\parallel \neq \sigma_\perp$ (function \texttt{\_sample\_grains\_bingham} in \texttt{hysterpy.dthm\_s2}) to collapse grains tighter onto the easy axis in the transverse direction. The activation of both within the joint DE bounds will be addressed in the next manuscript revision. It constitutes a direct demonstration that the modular DTHM-$S^{2}$ formulation reaches a regime that scalar hysteresis baselines (\JA, Preisach, \acs{DDD}, \acs{PI}) structurally cannot, because they carry no grain-orientation distribution and no compositional overlay mechanism. The fits reported here are reproducible from \texttt{scripts/fit\_hefmag\_blind\_init.py} (NO) and \texttt{scripts/fit\_hefmag\_go\_joint\_fast.py} (\acs{GO}). The blind-init helper is bundled as \texttt{hysterpy.blind\_init} in the open-source release.}
\end{figure}

\begin{figure}[htbp]
\centering
\includegraphics[width=\linewidth]{figures/llg_path_x_strategy_B.png}
\caption[LLG-direction infrastructure validation against Stoner--Wohlfarth analytical reference.]{\label{fig:llg_infrastructure}\textbf{Validation of the \ac{LLG}-direction infrastructure on cubic-anisotropy single-grain reference physics.} The \texttt{hysterpy.llg\_grain} module released with this work solves the Landau--Lifshitz--Gilbert equation for a per-grain three-vector magnetisation with cubic $\Kone$ anisotropy from the Fe-3\%Si single-crystal free energy. \textbf{(A)} Cubic anisotropy energy $U(\hat{\mathbf m})/\Kone$ sampled on the unit sphere $(\theta,\varphi)$ at $K_1\!=\!\qty{35}{\kilo\joule\per\meter\cubed}$ (single-crystal Fe-3\%Si literature value~\cite{cullity}). The four equatorial and two polar $\langle 100\rangle$ wells are resolved cleanly with the hard-axis $\langle 110\rangle$ ridges at \ang{45} between adjacent wells, in agreement with cubic symmetry. \textbf{(B)} Single-grain trajectories in the $(M_x/\Ms, M_y/\Ms)$ projection for five lab-field directions: \ang{0}, \ang{30}, \ang{45}, \ang{60}, and \ang{90} relative to the easy axis. At \ang{0} the trajectory collapses to a vertical reversal along $[001]\!\leftrightarrow\![00\bar{1}]$ (field on easy axis switches the pole). At \ang{90} the magnetization sweeps through the perpendicular cubic $\langle 100\rangle$ wells. Intermediate angles fan smoothly between, with no switching kinks. \textbf{(C)} Stoner--Wohlfarth astroid validation: numerical \ac{LLG} switching fields (red circles, swept over $\theta_\mathrm{e}\!\in\![\ang{5}, \ang{85}]$) compared against the analytical SW curve $H_x\!=\!H_K\cos^3 t$, $H_y\!=\!H_K\sin^3 t$. Maximum deviation $\Delta H_\mathrm{sw}/H_K\!=\!\qty{5.4}{\percent}$, attributed to finite-rate ramp ($8\tau$ per half-cycle versus the SW assumption of infinitely-slow drive). The implementation is therefore numerically validated against analytical reference physics. \textbf{Physics caveat for the Baghel \ang{45} application}: at $K_1\!=\!\qty{35}{\kilo\joule\per\meter\cubed}$ the cubic anisotropy field is $H_K\!=\!2K_1/(\mu_0 \Ms)\!\approx\!\qty{32}{\kilo\ampere\per\meter}$, far larger than the operating drive $H_\mathrm{peak}\!\approx\!\qty{0.6}{\kilo\ampere\per\meter}$ of the Baghel \ang{45} measurement. Coherent-rotation \ac{LLG} therefore cannot reproduce the measured $B\!\approx\!\qty{1.5}{\tesla}$ swing because magnetization reversal at this operating point is dominated by domain-wall motion rather than coherent rotation. The principled extension is therefore stochastic-\ac{LLG} with N\'eel--Brown thermal activation between the cubic $\langle 100\rangle$ wells, identified as the focus of the follow-on manuscript.}
\end{figure}

\clearpage

\begin{thebibliography}{10}
\urlstyle{rm}
\expandafter\ifx\csname url\endcsname\relax
  \def\url#1{\texttt{#1}}\fi
\expandafter\ifx\csname urlprefix\endcsname\relax\def\urlprefix{URL }\fi
\expandafter\ifx\csname doiprefix\endcsname\relax\def\doiprefix{}\fi
\providecommand{\bibinfo}[2]{#2}
\providecommand{\eprint}[2][]{\url{#2}}

\bibitem{jiles1986}
\bibinfo{author}{Jiles, D.~C.} \& \bibinfo{author}{Atherton, D.~L.}
\newblock \bibinfo{journal}{\bibinfo{title}{Theory of ferromagnetic
  hysteresis}}.
\newblock {\emph{\JournalTitle{Journal of Magnetism and Magnetic Materials}}}
  \textbf{\bibinfo{volume}{61}}, \bibinfo{pages}{48--60},
  \doiprefix\url{10.1016/0304-8853(86)90066-1} (\bibinfo{year}{1986}).

\bibitem{chwastek2009}
\bibinfo{author}{Chwastek, K.}
\newblock \bibinfo{journal}{\bibinfo{title}{Modelling of dynamic hysteresis
  loops using the {Jiles-Atherton} approach}}.
\newblock {\emph{\JournalTitle{Mathematical and Computer Modelling of Dynamical
  Systems}}} \textbf{\bibinfo{volume}{15}}, \bibinfo{pages}{95--105},
  \doiprefix\url{10.1080/13873950802432016} (\bibinfo{year}{2009}).

\bibitem{daloia2021}
\bibinfo{author}{D'Aloia, A.~G.}, \bibinfo{author}{{Di Francesco}, A.} \&
  \bibinfo{author}{{De Santis}, V.}
\newblock \bibinfo{journal}{\bibinfo{title}{A novel computational method to
  identify/analyze hysteresis loops of hard magnetic materials}}.
\newblock {\emph{\JournalTitle{Magnetochemistry}}}
  \textbf{\bibinfo{volume}{7}}, \doiprefix\url{10.3390/magnetochemistry7010010}
  (\bibinfo{year}{2021}).

\bibitem{mayergoyz}
\bibinfo{author}{Mayergoyz, I.~D.}
\newblock \emph{\bibinfo{title}{Mathematical Models of Hysteresis and Their
  Applications}} (\bibinfo{publisher}{Elsevier}, \bibinfo{year}{2003}),
  \bibinfo{edition}{2nd} edn.

\bibitem{krasnoselskii}
\bibinfo{author}{Krasnosel'ski{\v\i}, M.~A.} \&
  \bibinfo{author}{Pokrovski{\v\i}, A.~V.}
\newblock \emph{\bibinfo{title}{Systems with Hysteresis}}
  (\bibinfo{publisher}{Springer-Verlag}, \bibinfo{year}{1989}).

\bibitem{maghyst}
\bibinfo{author}{{Ilmenauer Mechatronik GmbH}}.
\newblock \bibinfo{title}{Mess- und {Pr{\"u}ftechnik}:
  {MagHyst\textsuperscript{\textregistered}} modular} (\bibinfo{year}{2023}).

\bibitem{zirka2008}
\bibinfo{author}{Zirka, S.~E.} \emph{et~al.}
\newblock \bibinfo{journal}{\bibinfo{title}{Generalization of the classical
  method for calculating dynamic hysteresis loops in grain-oriented electrical
  steels}}.
\newblock {\emph{\JournalTitle{IEEE Transactions on Magnetics}}}
  \textbf{\bibinfo{volume}{44}}, \bibinfo{pages}{2113--2126},
  \doiprefix\url{10.1109/TMAG.2008.2000662} (\bibinfo{year}{2008}).

\bibitem{baghel2015}
\bibinfo{author}{Baghel, A. P.~S.}, \bibinfo{author}{Gupta, A.},
  \bibinfo{author}{Chwastek, K.} \& \bibinfo{author}{Kulkarni, S.~V.}
\newblock \bibinfo{journal}{\bibinfo{title}{Comprehensive modelling of dynamic
  hysteresis loops in the rolling and transverse directions for transformer
  laminations}}.
\newblock {\emph{\JournalTitle{Physica B: Condensed Matter}}}
  \textbf{\bibinfo{volume}{462}}, \bibinfo{pages}{86--92},
  \doiprefix\url{10.1016/j.physb.2015.01.026} (\bibinfo{year}{2015}).

\bibitem{coffey2004}
\bibinfo{author}{Coffey, W.~T.}, \bibinfo{author}{Kalmykov, Y.~P.} \&
  \bibinfo{author}{Waldron, J.~T.}
\newblock \emph{\bibinfo{title}{The {Langevin} Equation: With Applications to
  Stochastic Problems in Physics, Chemistry and Electrical Engineering}}
  (\bibinfo{publisher}{World Scientific}, \bibinfo{year}{2004}),
  \bibinfo{edition}{2nd} edn.

\bibitem{daniel2008}
\bibinfo{author}{Daniel, L.}, \bibinfo{author}{Hubert, O.},
  \bibinfo{author}{Buiron, N.} \& \bibinfo{author}{Billardon, R.}
\newblock \bibinfo{journal}{\bibinfo{title}{Reversible magneto-elastic
  behavior: A multiscale approach}}.
\newblock {\emph{\JournalTitle{Journal of the Mechanics and Physics of
  Solids}}} \textbf{\bibinfo{volume}{56}}, \bibinfo{pages}{1018--1042},
  \doiprefix\url{10.1016/j.jmps.2007.06.003} (\bibinfo{year}{2008}).

\bibitem{liluo2024}
\bibinfo{author}{Li, Y.} \& \bibinfo{author}{Luo, Z.}
\newblock \bibinfo{journal}{\bibinfo{title}{A comparative study of {J--A} and
  {Preisach} improved models for core loss calculation under {DC} bias}}.
\newblock {\emph{\JournalTitle{Scientific Reports}}}
  \textbf{\bibinfo{volume}{14}}, \bibinfo{pages}{9409},
  \doiprefix\url{10.1038/s41598-024-55155-w} (\bibinfo{year}{2024}).

\bibitem{ramesh1997}
\bibinfo{author}{Ramesh, A.}, \bibinfo{author}{Jiles, D.~C.} \&
  \bibinfo{author}{Bi, Y.}
\newblock \bibinfo{journal}{\bibinfo{title}{Generalization of hysteresis
  modeling to anisotropic materials}}.
\newblock {\emph{\JournalTitle{Journal of Applied Physics}}}
  \textbf{\bibinfo{volume}{81}}, \bibinfo{pages}{5585--5587},
  \doiprefix\url{10.1063/1.364843} (\bibinfo{year}{1997}).

\bibitem{hefmag_zenodo}
\bibinfo{author}{{HEFMAG Consortium}}, \bibinfo{author}{{INRiM}},
  \bibinfo{author}{{PTB}}, \bibinfo{author}{{NPL}} \& \bibinfo{author}{{CMI}}.
\newblock \bibinfo{title}{Wideband hysteresis and complex-permeability dataset
  for grain-oriented and non-oriented Fe-Si, amorphous and nanocrystalline
  ribbons (EMPIR 19ENG06)}, \doiprefix\url{10.5281/zenodo.5910892}
  (\bibinfo{year}{2022}).

\bibitem{condegarrido2025}
\bibinfo{author}{Conde~Garrido, J.~M.}, \bibinfo{author}{Ugarte~Valdivielso,
  J.}, \bibinfo{author}{Aizpurua, J.~I.}, \bibinfo{author}{Barrenetxea
  I\~{n}arra, M.} \& \bibinfo{author}{Silveyra, J.~M.}
\newblock \bibinfo{journal}{\bibinfo{title}{Blind efficient method for
  optimizing {Jiles}--{Atherton} model parameters}}.
\newblock {\emph{\JournalTitle{IEEE Transactions on Magnetics}}}
  \doiprefix\url{10.1109/TMAG.2025.3632479} (\bibinfo{year}{2025}).

\bibitem{cullity}
\bibinfo{author}{Cullity, B.~D.} \& \bibinfo{author}{Graham, C.~D.}
\newblock \emph{\bibinfo{title}{Introduction to Magnetic Materials}}
  (\bibinfo{publisher}{Wiley-IEEE Press}, \bibinfo{year}{2008}),
  \bibinfo{edition}{2nd} edn.

\end{thebibliography}


\end{document}
