# The Michaelis–Menten reaction at low substrate concentrations:
Pseudo-first-order kinetics and conditions for timescale separation

Justin Eilertsen
Mathematical Reviews
American Mathematical Society
416 4th Street
Ann Arbor, MI, 48103
e-mail<EMAIL_ADDRESS>

Santiago Schnell
Department of Biological Sciences and
Department of Applied and Computational Mathematics and Statistics
University of Notre Dame
Notre Dame, IN 46556
e-mail<EMAIL_ADDRESS>

Sebastian Walcher
Mathematik A, RWTH Aachen
D-52056 Aachen, Germany
e-mail<EMAIL_ADDRESS>
###### Abstract

We demonstrate that the Michaelis–Menten reaction mechanism can be accurately
approximated by a linear system when the initial substrate concentration is
low. This leads to pseudo-first-order kinetics, simplifying mathematical
calculations and experimental analysis. Our proof utilizes a monotonicity
property of the system and Kamke’s comparison theorem. This linear
approximation yields a closed-form solution, enabling accurate modeling and
estimation of reaction rate constants even without timescale separation.
Building on prior work, we establish that the sufficient condition for the
validity of this approximation is $s_{0}\ll K$, where $K=k_{2}/k_{1}$ is the
Van Slyke–Cullen constant. This condition is independent of the initial enzyme
concentration. Further, we investigate timescale separation within the linear
system, identifying necessary and sufficient conditions and deriving the
corresponding reduced one-dimensional equations.

MSC (2020): 92C45, 34C12, 34E15
Keywords: Comparison principle, pseudo-first-order kinetics, total quasi-
steady state approximation, monotone dynamical system

##  1 Introduction and overview of results

The irreversible Michaelis–Menten reaction mechanism describes the mechanism
of action of a fundamental reaction for biochemistry. Its governing equations
are:

$\begin{array}[]{rccrclcl}\dot{s}&=&f_{1}(s,c):=&-k_{1}e_{0}s&+&(k_{1}s+k_{-1})c&&\\\
\dot{c}&=&f_{2}(s,c):=&k_{1}e_{0}s&-&(k_{1}s+k_{-1}+k_{2})c&&\\\ \end{array}$
(1)

with positive parameters $e_{0}$, $s_{0}$, $k_{1}$, $k_{-1}$, $k_{2}$, and
typical initial conditions $s(0)=s_{0}$, $c(0)=0$. Understanding the behavior
of this system continues to be crucial for both theoretical and experimental
purposes, such as identifying the rate constants $k_{1}$, $k_{-1}$ and
$k_{2}$. However, the Michaelis–Menten system (2) cannot be solved via
elementary functions, hence approximate simplifications have been proposed
since the early 1900’s.

In 1913, Michaelis and Menten [12] introduced the partial equilibrium
assumption for enzyme-substrate complex formation, effectively requiring a
small value for the catalytic rate constant $k_{2}$. Later, for systems with
low initial enzyme concentration ($e_{0}$), Briggs and Haldane [3] derived the
now-standard quasi-steady-state reduction (see Segel and Slemrod [16] for a
systematic analysis). Both these approaches yield an asymptotic reduction to a
one-dimensional equation. A mathematical justification is obtained from
singular perturbation theory. Indeed, most literature on the Michaelis–Menten
system focuses on dimensionality reduction, with extensive discussions on
parameter combinations ensuring this outcome.

We consider a different scenario: low initial substrate concentration
($s_{0}$). This presents a distinct case from both the low-enzyme and partial-
equilibrium scenarios, which yield one-dimensional reductions in appropriate
parameter ranges. Instead, with low substrate, we obtain a different
simplification:

$\begin{pmatrix}\dot{s}\\\
\dot{c}\end{pmatrix}=\begin{pmatrix}-k_{1}e_{0}&k_{-1}\\\
k_{1}e_{0}&-(k_{-1}+k_{2})\end{pmatrix}\cdot\begin{pmatrix}s\\\
c\end{pmatrix}.$ (2)

This linear Michaelis–Menten system (2) was first studied by Kasserra and
Laidler [11], who proposed the condition of excess initial enzyme
concentration ($e_{0}\gg s_{0}$) for the validity of the linearization of (1)
to (2). This aligns with the concept of pseudo-first-order kinetics in
chemistry, where the concentration of one reactant is so abundant that it
remains essentially constant [17].

Pettersson [14] further investigated the linearization, adding the assumption
that any complex concentration accumulated during the transient period remains
too small to significantly impact enzyme or product concentrations. Later,
Schnell and Mendoza [15] refined the validity condition $s_{0}\ll K_{M}$
(where $K_{M}=(k_{-1}+k_{2})/k_{1}$ is the Michaelis constant) to be
sufficient for the linearization (alias dictus the application of pseudo-
first-order kinetics to the Michaelis–Menten reaction). However, Pedersen and
Bersani [13] found this condition overly conservative. They introduced the
total substrate concentration $\bar{s}=s+c$ (used in the total quasi-steady
state approximation [2]) and proposed the condition $s_{0}\ll K_{M}+e_{0}$,
reconciling the Schnell–Mendoza and Kasserra–Laidler conditions.

In this paper, we demonstrate that the solutions of the Michaelis–Menten
system (1) admit a global approximation by solutions of the linear
differential equation (2) whenever the intial substrate concentration is
small. Our proof utilizes Kamke’s comparison theorem for cooperative
differential equations.111This concept is distinct from cooperativity in
biochemistry. This theorem allows us to establish upper and lower estimates
for the solution components of (1) via suitable linear systems. Importantly,
these bounds are valid whenever $s_{0}<K$, where $K=k_{2}/k_{1}$ is the Van
Slyke–Cullen constant, and the estimates are tight whenever $s_{0}\ll K$. This
effectively serves as sufficient condition for the validity of the linear
approximation Michaelis–Menten system (2). We further prove that solutions of
the original Michaelis–Menten system (1) converge to solutions of the linear
system (2) as $s_{0}\to 0$, uniformly for all positive times. Moreover, the
asymptotic rate of convergence is of order $s_{0}^{2}$.

Additionally, we explore the problem of of timescale separation within the
linear Michaelis–Menten system (2). While a linear approximation simplifies
the system, it does not guarantee a clear separation of timescales. In fact,
there are cases where no universal accurate one-dimensional approximation
exists. Since both eigenvalues of the matrix in (2) are real and negative, the
slower timescale will dominate the solution’s behavior. True timescale
separation occurs if and only if these eigenvalues differ significantly in
magnitude. Importantly, we find that this separation is not always guaranteed.
There exist parameter combinations where a one-dimensional approximation would
lack sufficient accuracy. Building upon results from [5], we establish
necessary and sufficient conditions for timescale separation and derive the
corresponding reduced one-dimensional equations.

##  2 Low substrate: Approximation by a linear system

To set the stage, we recall some facts from the theory of cooperative
differential equation systems. For a comprehensive account of the theory, we
refer to the H.L. Smith’s monograph [18, Chapter 3] on monotone dynamical
systems. Consider the standard ordering, here denoted by $\prec$, on
$\mathbb{R}^{n}$, given by

$\begin{pmatrix}x_{1}\\\ \vdots\\\
x_{n}\end{pmatrix}\prec\begin{pmatrix}y_{1}\\\ \vdots\\\
y_{n}\end{pmatrix}\Longleftrightarrow x_{i}\leq y_{i}\text{ for }1\leq i\leq
n.$

Now let a differential equation

$\dot{x}=f(x)\quad\text{on
}U\subseteq\mathbb{R}^{n},\quad\emptyset\not=U\text{ open }$ (3)

be given, with $f$ continuously differentiable. We denote by $F(t,y)$ the
solution of the initial value problem (3) with $x(0)=y$, and call $F$ the
local flow of the system.

Given a convex subset $D\subseteq U$, we say that (3) is cooperative on $D$ if
all non-diagonal elements of the Jacobian are nonnegative; symbolically

$Df(x)=\begin{pmatrix}*&+&\cdots&&+\\\ +&*&+&\cdots&+\\\
\vdots&+&\ddots&&\vdots\\\ \vdots&&&\ddots&+\\\
+&\cdots&\cdots&+&*\end{pmatrix}\text{ for all }x\in D.$

As shown by M.W. Hirsch ([8, and subsequent papers]) and other researchers
(see [18, and the references therein]), cooperative systems have special
qualitative features. Qualitative properties of certain monotone chemical
reaction networks were investigated by De Leenheer et al. [4]. In the present
paper, we consider a different feature of such systems. Our interest lies in
the following simplified version of Kamke’s comparison theorem [9].222Kamke
[9] considered non-autonomous systems with less restrictive properties. In
particular the statement also holds on P-convex sets. A subset $D$ of
$\mathbb{R}^{n}$ is called P-convex if, for each pair $y,\,z\in D$ with
$y\prec z$, the line segment connecting $y$ and $z$ is contained in $D$. Thus,
P-convexity is a weaker property than convexity.

######  Proposition 1.

Let $D\subseteq U$ be a convex positively invariant set for (3), and assume
that this system is cooperative on $D$. Moreover let $\dot{x}=g(x)$ be defined
on $U$, with continuously differentiable right hand side, and local flow $G$.

  * •

If $f(x)\prec g(x)$ for all $x\in D$, and $y\prec z$, then $F(t,y)\prec
G(t,z)$ for all $t\geq 0$ such that both solutions exist.

  * •

If $g(x)\prec f(x)$ for all $x\in D$, and $z\prec y$, then $G(t,z)\prec
F(t,y)$ for all $t\geq 0$ such that both solutions exist.

###  2.1 Application to the Michaelis–Menten systems (1) and (2)

For the following discussions, we recall the relevant derived parameters

$K_{S}:=k_{-1}/k_{1},\quad K_{M}:=(k_{-1}+k_{2})/k_{1}\quad
K:=k_{2}/k_{1}=K_{M}-K_{S},$ (4)

where $K_{S}$ is the complex equilibrium constant, $K_{M}$ is the Michaelis
constant, and $K$ is the van Slyke-Cullen constant.

The proof of the following facts is obvious.

######  Lemma 1.

  1. (a)

The Jacobian of system (1) is equal to

$\begin{pmatrix}\ast&k_{1}s+k_{-1}\\\ k_{1}(e_{0}-c)&\ast\end{pmatrix}.$

The off-diagonal entries are $\geq 0$ in the physically relevant subset

$D=\left\\{(s,c);\,s\geq 0,\,e_{0}\geq c\geq 0,\,s+c\leq s_{0}\right\\}$

of the phase space, and $D$ is positively invariant and convex. Thus system
(1) is cooperative in $D$.

  2. (b)

In $D$, one has

$\begin{array}[]{rcccl}f_{1}(s,c)&\geq&g_{1}(s,c)&:=&-k_{1}e_{0}s+k_{-1}c\\\
f_{2}(s,c)&\geq&g_{2}(s,c)&:=&k_{1}e_{0}s-(k_{1}s_{0}+k_{-1}+k_{2})c,\end{array}$

hence by Kamke’s comparison theorem the solution of the linear system with
matrix

$G:=\begin{pmatrix}-k_{1}e_{0}&k_{-1}\\\
k_{1}e_{0}&-(k_{1}s_{0}+k_{-1}+k_{2})\end{pmatrix}=k_{1}\begin{pmatrix}-e_{0}&K_{S}\\\
e_{0}&-K_{M}(1+s_{0}/K_{M})\end{pmatrix}$

and initial value in $D$ provides component-wise lower estimates for the
solution of (1) with the same initial value.

  3. (c)

In $D$, one has

$\begin{array}[]{rcccl}f_{1}(s,c)&\leq&h_{1}(s,c)&:=&-k_{1}e_{0}s+(k_{1}s_{0}+k_{-1})c\\\
f_{2}(s,c)&\leq&h_{2}(s,c)&:=&k_{1}e_{0}s-(k_{-1}+k_{2})c,\end{array}$

hence by Kamke’s comparison theorem the solution of the linear system with
matrix

$H:=\begin{pmatrix}-k_{1}e_{0}&k_{1}s_{0}+k_{-1}\\\
k_{1}e_{0}&-(k_{-1}+k_{2})\end{pmatrix}=k_{1}\begin{pmatrix}-e_{0}&K_{S}(1+s_{0}/K_{S})\\\
e_{0}&-K_{M}\end{pmatrix}$

and initial value in $D$ provides component-wise upper estimates for the
solution of (1) with the same initial value.

######  Remark 1.

The upper estimates become useless in the case $k_{1}s_{0}>k_{2}$, because
then one eigenvalue of $H$ becomes positive. For viable upper and lower bounds
one requires that $s_{0}<K=k_{2}/k_{1}$.

The right hand sides of both linear comparison systems from parts (b) and (c)
of the Lemma converge to

$Df(0)=\begin{pmatrix}-k_{1}e_{0}&k_{-1}\\\
k_{1}e_{0}&-(k_{-1}+k_{2})\end{pmatrix}=k_{1}\begin{pmatrix}-e_{0}&K_{S}\\\
e_{0}&-K_{M}\end{pmatrix}$ (5)

as $s_{0}\to 0$. We will use this observation to show that a solution of (1)
with initial value in $D$ converges to the solution of (2) with the same
initial value on the time interval $[0,\,\infty)$, as $s_{0}\to 0$.

Figure 1: Illustration of Proposition 2. The solution to the Michaelis–Menten
system (1) converges to the solution of the linear Michaelis–Menten system (2)
as $s_{0}\to 0$. In all panels, the solid black curve is the numerical
solution to the Michaelis–Menten system (1). The thick yellow curve is the
numerical solution to linear Michaelis–Menten system (2). The dashed/dotted
curve is the numerical solution to the linear system defined by matrix $G$.
The dotted curve is the numerical solution to the linear system defined by
matrix $H$. All numerical simulations where carried out with the following
parameters (in arbitrary units): $k_{1}=k_{2}=k_{-1}=e_{0}=1$. In all panels,
the solutions have been numerically-integrated over the domain $t\in[0,T],$
where $T$ is selected to be long enough to ensure that the long-time dynamics
are sufficiently captured. For illustrative purposes, the horizontal axis (in
all four panels) has been scaled by $T$ so that the scaled time, $t/T$,
assumes values in the unit interval: $\frac{t}{T}\in[0,1]$. Top Left: The
numerically-obtained time course of $s$ with $s_{0}=0.5$ and $c(0)=0.0$. Top
Right: The numerically-obtained time course of $c$ with $s_{0}=0.5$ and
$c(0)=0.0$. Bottom Left: The numerically-obtained time course of $s$ with
$s_{0}=0.1$ and $c(0)=0.0$. Bottom Right: The numerically-obtained time course
of $c$ with $s_{0}=0.1$ and $c(0)=0.0$. Observe that the solution components
of (2) become increasingly accurate approximations to the solution components
of (1) as $s_{0}$ decreases.

###  2.2 Comparison Michaelis–Menten systems (1) and (2)

Both $G$ and $H$, as well as $Df(0)$, are of the type

$\begin{pmatrix}-\alpha&\beta\\\ \alpha&-\gamma\end{pmatrix}$ (6)

with all entries $>0$, and $\gamma>\beta$ (to ensure usable estimates). The
trace equals $-(\alpha+\gamma)$, the determinant equals
$\alpha(\gamma-\beta)$, and the discriminant is

$\Delta=(\alpha+\gamma)^{2}-4\alpha(\gamma-\beta)=(\alpha-\gamma)^{2}+4\alpha\beta.$

The eigenvalues are

$\lambda_{1,2}=\frac{1}{2}\left(-(\alpha+\gamma)\pm\sqrt{\Delta}\right)$

with eigenvectors

$v_{1,2}=\begin{pmatrix}\beta\\\
\frac{1}{2}(\alpha-\gamma\pm\sqrt{\Delta}).\end{pmatrix}$

The following is now obtained in a straightforward (if slightly tedious)
manner. We specialize the initial value to the typical case, with no complex
present at $t=0$.

######  Lemma 2.

The solution of the linear differential equation with matrix (6) and initial
value $\begin{pmatrix}s_{0}\\\ 0\end{pmatrix}$ is equal to

$\begin{array}[]{rccl}\begin{pmatrix}\widetilde{s}\\\
\widetilde{c}\end{pmatrix}&=&\dfrac{s_{0}}{2\beta\sqrt{\Delta}}&\left[(\alpha-\gamma-\sqrt{\Delta})\begin{pmatrix}\beta\\\
\frac{1}{2}(\alpha-\gamma+\sqrt{\Delta})\end{pmatrix}\,\exp(\left(\frac{1}{2}(-\alpha-\gamma+\sqrt{\Delta})t\right)\right.\\\
&&+&\left.(-\alpha+\gamma-\sqrt{\Delta})\begin{pmatrix}\beta\\\
\frac{1}{2}(\alpha-\gamma-\sqrt{\Delta})\end{pmatrix}\,\exp(\left(\frac{1}{2}(-\alpha-\gamma-\sqrt{\Delta})t\right)\right]\\\
&=&\dfrac{s_{0}}{2\sqrt{\Delta}}&\left[\begin{pmatrix}-\alpha+\gamma+\sqrt{\Delta}\\\
2\alpha\end{pmatrix}\,\exp(\left(\frac{1}{2}(-\alpha-\gamma+\sqrt{\Delta})t\right)\right.\\\
&&+&\left.\begin{pmatrix}\alpha-\gamma+\sqrt{\Delta}\\\
-2\alpha\end{pmatrix}\,\exp(\left(\frac{1}{2}(-\alpha-\gamma-\sqrt{\Delta})t\right)\right]\\\
\end{array}$

This enables us to write down explicitly the solution of (2), as well as upper
and lower estimates, in terms of familiar constants. We use Lemma 2 with
$\alpha=k_{1}e_{0}$, $\beta=k_{1}K_{S}$ and $\gamma=k_{1}K_{M}$, so
$\Delta=k_{1}^{2}\left((K_{M}-e_{0})^{2}+4K_{S}e_{0}\right)$, to obtain

$\begin{array}[]{rcl}\begin{pmatrix}s^{*}\\\
c^{*}\end{pmatrix}&=&\dfrac{s_{0}\exp(-\frac{1}{2}((K_{M}+e_{0})-\sqrt{(K_{M}-e_{0})^{2}+4K_{S}e_{0}})k_{1}t)}{2\sqrt{(K_{M}-e_{0})^{2}+4K_{S}e_{0}}}\times\\\
&&\qquad\qquad\qquad\begin{pmatrix}+(K_{M}-e_{0})+\sqrt{(K_{M}-e_{0})^{2}+4K_{S}e_{0}}\\\
2e_{0}\end{pmatrix}\\\
&+&\dfrac{s_{0}\exp(-\frac{1}{2}((K_{M}+e_{0})+\sqrt{(K_{M}-e_{0})^{2}+4K_{S}e_{0}})k_{1}t)}{2\sqrt{(K_{M}-e_{0})^{2}+4K_{S}e_{0}}}\times\\\
&&\qquad\qquad\qquad\begin{pmatrix}-(K_{M}-e_{0})+\sqrt{(K_{M}-e_{0})^{2}+4K_{S}e_{0}}\\\
-2e_{0}\end{pmatrix}\\\ \end{array}$ (7)

as the solution of (2), after some simplifications.

######  Proposition 2.

To summarize:

  1. (a)

The solution $\begin{pmatrix}s^{*}\\\ c^{*}\end{pmatrix}$ of the linear
differential equation (2) with initial value $\begin{pmatrix}s_{0}\\\
0\end{pmatrix}$ is given by (7).

  2. (b)

The solution $\begin{pmatrix}s_{\rm low}\\\ c_{\rm low}\end{pmatrix}$ of the
linear differential equation with matrix $G$ and initial value
$\begin{pmatrix}s_{0}\\\ 0\end{pmatrix}$ is obtained from replacing $K_{M}$ by
$K_{M}+s_{0}$ in (7).

  3. (c)

The solution $\begin{pmatrix}s_{\rm up}\\\ c_{\rm up}\end{pmatrix}$ of the
linear differential equation with matrix $H$ and initial value
$\begin{pmatrix}s_{0}\\\ 0\end{pmatrix}$ is obtained from replacing $K_{S}$ by
$K_{S}+s_{0}$ in (7).

  4. (d)

Given that $s_{0}<K$, for all $t\geq 0$ one has the inequalities

$\begin{array}[]{cc}s_{\rm up}\geq s\geq s_{\rm low},&\quad s_{\rm up}\geq
s^{*}\geq s_{\rm low},\\\ c_{\rm up}\geq c\geq c_{\rm low},&\quad c_{\rm
up}\geq c^{*}\geq c_{\rm low};\\\ \end{array}$ (8)

where $\begin{pmatrix}s\\\ c\end{pmatrix}$ denotes the solution of (1) with
initial value $\begin{pmatrix}s_{0}\\\ 0\end{pmatrix}$.

  5. (e)

Let $\mathcal{M}$ be a compact subset of the open positive orthant
$\mathbb{R}^{4}_{>0}$, abbreviate $\pi:=(e_{0},k_{1},k_{-1},k_{2})^{\rm tr}$
and $K^{*}:=\min\\{k_{2}/k_{1};\,\pi\in\mathcal{M}\\}$. Then, there exists a
dimensional constant $C$ (with dimension concentration-1), depending only on
$\mathcal{M}$, such that for all $\pi\in{\mathcal{M}}$ and for all $s_{0}$
with $0<s_{0}\leq K^{*}/2$ the estimates

$\begin{array}[]{rcl}\left|\dfrac{s-s^{*}}{s_{0}}\right|&\leq&C\cdot s_{0}\\\
\left|\dfrac{c-c^{*}}{s_{0}}\right|&\leq&C\cdot s_{0}\\\ \end{array}$

hold for all $t\in[0,\infty)$. Thus, informally speaking, the approximation
errors of $s$ by $s^{*}$, and of $c$ by $c^{*}$, are of order $s_{0}^{2}$.

###### Proof.

The first three items follow by straightforward calculations. As for part (d),
the first column of (8) is just a restatement of Lemma 1, while the second
follows from the observation that (mutatis mutandis) the right-hand side of
(2) also obeys the estimates in parts (b) and (c) of this Lemma.

There remains part (e). In view of part (d) it suffices to show that

$\begin{array}[]{rcl}\left|\dfrac{s_{\rm up}-s_{\rm
low}}{s_{0}}\right|&\leq&C\cdot s_{0}\\\ \left|\dfrac{c_{\rm up}-c_{\rm
low}}{s_{0}}\right|&\leq&C\cdot s_{0}\\\ \end{array}$

for some constant and, in turn, it suffices to show such estimates for both
$s_{\rm up}-s^{*}$, $s^{*}-s_{\rm low}$, and $c_{\rm up}-c^{*}$, $c^{*}-c_{\rm
low}$. We will outline the relevant steps, pars pro toto, for the upper
estimates.

By (7) and part (c), one may write

$\dfrac{1}{s_{0}}\begin{pmatrix}s_{\rm up}\\\ c_{\rm
up}\end{pmatrix}=B_{1}(s_{0},\pi)\exp(-A_{1}(s_{0},\pi)t)+B_{2}(s_{0},\pi)\exp(-A_{2}(s_{0},\pi)t),$

with the $B_{i}$ and $A_{i}$ continuously differentiable in a neighborhood of
$[0,K^{*}/2]\times\mathcal{M}\times[0,\infty)$. Moreover, one sees

$\dfrac{1}{s_{0}}\begin{pmatrix}s^{*}\\\
c^{*}\end{pmatrix}=B_{1}(0,\pi)\exp(-A_{1}(0,\pi)t)+B_{2}(0,\pi)\exp(-A_{2}(0,\pi)t).$

It suffices to show that (with the maximum norm $\|\cdot\|$)

$\|B_{i}(s_{0},\pi)\exp(-A_{i}(s_{0},\pi)t)-B_{i}(0,\pi)\exp(-A_{i}(0,\pi)t)\|\leq
s_{0}\cdot\,{\rm const.}$

in $\widetilde{\mathcal{M}}:=[0,K^{*}/2]\times\mathcal{M}\times[0,\infty)$,
for $i=1,\,2$.

From compactness of $[0,K^{*}/2]\times{\mathcal{M}}$, continuous
differentiability and the explicit form of $A_{i}$ and $B_{i}$ one obtains
estimates, with the parameters and variables in $\widetilde{\mathcal{M}}$:

  * •

There is $A_{*}>0$ so that $|A_{i}(s_{0},\pi)|\geq A_{*}$.

  * •

There is $B^{*}>0$ so that $\|B_{i}(s_{0},\pi)\|\leq B^{*}$.

  * •

There are continuous $\widetilde{B}_{i}$ so that
$B_{i}(s_{0},\pi)-B_{i}(0,\pi)=s_{0}\cdot\widetilde{B}_{i}(s_{0},\pi)$
(Taylor), hence there is $B^{**}>0$ so that
$\|B_{i}(s_{0},\pi)-B_{i}(0,\pi)\|\leq s_{0}\cdot B^{**}$.

  * •

By the mean value theorem there exists $\sigma$ between $0$ and $s_{0}$ so
that

$\exp\left(-A_{i}(s_{0},\pi)t\right)-\exp\left(-A_{i}(0,\pi)t\right)=-t\cdot\dfrac{\partial
A_{i}}{\partial
s_{0}}\,(\sigma,\pi)\cdot\exp\left(-A_{i}(\sigma,\pi)t\right)\,s_{0}.$

So, with the constant $A^{**}$ satisfying $A^{**}\geq|\dfrac{\partial
A_{i}}{\partial s_{0}}|$ for all arguments in $\widetilde{\mathcal{M}}$ one
gets

$|\exp\left(-A_{i}(s_{0},\pi)t\right)-\exp\left(-A_{i}(0,\pi)t\right)|\leq
A^{**}\cdot t\exp\left(-A_{i}(\sigma,\pi)t\right)\,s_{0}\leq A^{**}\cdot
t\exp\left(-A_{*}t\right)\leq\dfrac{A^{**}}{A_{*}}s_{0}.$

So

$\begin{array}[]{cl}&\|B_{i}(s_{0},\pi)\exp(-A_{i}(s_{0},\pi)t)-B_{i}(0,\pi)\exp(-A_{i}(0,\pi)t)\|\\\
\leq&\|B_{i}(s_{0},\pi)-B_{i}(0,\pi)\|\cdot\exp(-A_{i}(s_{0},\pi)t)\\\
&+\|B_{i}(0,\pi)\|\cdot|\exp(-A_{i}(s_{0},\pi)t)-\exp(-A_{i}(0,\pi)t)|\\\
\leq&s_{0}\cdot(B^{**}+B^{*}\cdot A^{**}/A_{*}),\\\ \end{array}$

and this proves the assertion. ∎

For illustrative purposes, a numerical example is given in Figure 1.

######  Remark 2.

Our result should not be mistaken for the familiar fact that in the vicinity
of the stationary point $0$ the Michaelis–Menten system (1) is (smoothly)
equivalent to its linearization (2); see, e.g. Sternberg [19, Theorem 2]. The
point of the Proposition is that the estimate, and the convergence statement,
hold globally.

######  Remark 3.

Some readers may prefer a dimensionless parameter (and symbols like $\ll$) to
gauge the accuracy of the approximation. In view of Lemma 1 and equation (5)
it seems natural to choose $\varepsilon:=s_{0}/K_{S}<s_{0}/K_{M}$, which
estimates the convergence of $G$ and $H$ to $Df(0)$, and set up the criterion
$s_{0}/K_{S}\ll 1$. But still the limit $s_{0}\to 0$ has to be considered.

######  Remark 4.

Since explicit expressions are obtainable from (7) for the $A_{i}$ and $B_{i}$
in the proof of part (e), one could refine it to obtain quantitative
estimates. This will not be pursued further here.

###  2.3 Timescales for the linear Michaelis–Menten system (2)

In contrast to familiar reduction scenarios for the Michaelis–Menten system
(1), low initial substrate linearization (2) does not automatically imply a
separation of timescales. Timescale separation depends on further conditions
that we discuss next.

Recall the constants introduced in (4) and note

$\Delta=k_{1}^{2}\left((K_{M}+e_{0})^{2}-4Ke_{0}\right).$

In this notation, the eigenvalues of the matrix in the linearized system (2)
are

$\lambda_{1,2}=\frac{k_{1}}{2}(K_{M}+e_{0})\left(-1\pm\sqrt{1-\dfrac{4Ke_{0}}{(K_{M}+e_{0})^{2}}}\right).$
(9)

This matrix is the Jacobian at the stationary point $0$, thus it is possible
to use the results from earlier work [5].

######  Remark 5.

Since timescales are inverse absolute eigenvalues, by [5, subsection 3.3], (9)
shows:

  1. (a)

The timescales are about equal whenever

$4Ke_{0}\approx(K_{M}+e_{0})^{2},$

which is the case, notably, when $K_{S}\approx 0$ (so $K_{M}\approx K$) and
$K_{M}\approx e_{0}$.

  2. (b)

On the other hand, one sees from (9) that a significant timescale separation
exists whenever (loosely speaking, employing a widely used symbol)

$\dfrac{4Ke_{0}}{(K_{M}+e_{0})^{2}}\ll 1.$ (10)

Notably this is the case when $e_{0}/K_{M}$ is small, or when $K/K_{M}$ is
small.

  3. (c)

Moreover, condition (10) is satisfied in a large part of parameter space (in a
well-defined sense) as shown in Eilertsen et al. [5, subsection 3.3, in
particular Figure 2].

######  Proposition 3.

We take a closer look at the scenario with significant timescale separation,
stating the results in a somewhat loose language.

  1. (a)

Given significant timescale separation as in (10), the slow eigenvalue is
approximated by

$\lambda_{1}\approx-
k_{1}\cdot\dfrac{Ke_{0}}{K_{M}+e_{0}}=-\dfrac{k_{2}e_{0}}{K_{M}+e_{0}},$ (11)

and the reduced equation is

$\dfrac{d}{dt}\begin{pmatrix}s\\\
c\end{pmatrix}=-\dfrac{k_{2}e_{0}}{K_{M}+e_{0}}\begin{pmatrix}s\\\
c\end{pmatrix}.$ (12)

  2. (b)

Given significant timescale separation as in (10), the eigenspace for the slow
eigenvalue is asymptotic to the subspace spanned by

$\widehat{v}_{1}:=\begin{pmatrix}K_{M}-\dfrac{Ke_{0}}{K_{M}+e_{0}}\\\
e_{0}\end{pmatrix}.$ (13)

###### Proof.

Given (10) one has

$\sqrt{(K_{M}+e_{0})^{2}-4Ke_{0}}=(K_{M}+e_{0})\sqrt{1-\dfrac{4Ke_{0}}{(K_{M}+e_{0})^{2}}}\approx(K_{M}+e_{0})\left(1-\dfrac{2Ke_{0}}{(K_{M}+e_{0})^{2}}\right)$

by Taylor approximation, from which part (a) follows.

Moreover, according to (7), the direction of the slow eigenspace is given by

$v_{1}^{*}=\begin{pmatrix}(K_{M}-e_{0})+\sqrt{(K_{M}-e_{0})^{2}+4K_{S}e_{0}}\\\
2e_{0}\end{pmatrix}=\begin{pmatrix}(K_{M}-e_{0})+\sqrt{(K_{M}+e_{0})^{2}-4Ke_{0}}\\\
2e_{0}\end{pmatrix}.$

By the same token, we have the approximation

$\begin{array}[]{rcc}v_{1}^{*}&\approx&\begin{pmatrix}(K_{M}-e_{0})+(K_{M}+e_{0})\left(1-\dfrac{2e_{0}K}{(K_{M}+e_{0})^{2}}\right)\\\
2e_{0}\end{pmatrix}\\\ &=&2\begin{pmatrix}K_{M}-\dfrac{Ke_{0}}{K_{M}+e_{0}}\\\
e_{0}\end{pmatrix}=2\widehat{v}_{1},\end{array}$ (14)

which shows part (b). ∎

######  Remark 6.

One should perhaps emphasize that Proposition 3 indeed describes a singular
perturbation reduction when $\lambda_{1}\to 0$. This can be verified from the
general reduction formula in Goeke and Walcher [7] (note that the
Michaelis–Menten system (2) is not in standard form with separated slow and
fast variables), with elementary computations.

Finally we note that Proposition 3 is applicable to the discussion of the
linearized total quasi-steady-state approximation. Details are discussed in
[6].

##  3 Conclusion

Our comprehensive analysis of the Michaelis–Menten system under low initial
substrate concentration addresses an important gap in the literature. We
obtain the sufficient333Necessity can easily be verified. See also the
discussion in [6] for examples. condition, $s_{0}\ll K$ for the global
approximation of the Michaelis–Menten system (1) by the linear
Michaelis–Menten system (2). This simplification, known as the pseudo-first-
order approximation, is widely used in transient kinetics experiments to
measure enzymatic reaction rates. Moreover, given the biochemical significance
of this reaction mechanism and the renewed interest in the application of the
total quasi-steady state approximation in pharmacokinetics experiments under
low initial substrate concentration [1, 20], the linearization analysis holds
clear practical relevance.

From a mathematical perspective, it is intriguing to observe a simplification
that arises independently of timescale separation or invariant manifolds. The
natural next step is to explore this approach with other familiar enzyme
reaction mechanism (e.g., those involving competitive and uncompetitive
inhibition, or cooperativity [see, for example, Keener and Sneyd [10]]).
However, a direct application of Kamke’s comparison theorem is not always
feasible. Enzyme catalyzed reaction mechanism with competitive or
uncompetitive inhibition do not yield cooperative differential equations, and
standard cooperative networks do so only within specific parameter ranges.
Investigating these systems under low initial substrate conditions will
require the development of further mathematical tools.

## Data availability

We do not analyse or generate any datasets, because our work proceeds within a
theoretical and mathematical approach.

## References

  * [1] H.-M. Back, H.-Y. Yun, S.K. Kim, J.K. Kim: Beyond the Michaelis-Menten: Accurate Prediction of In Vivo Hepatic Clearance for Drugs With Low KM. Clin. Transl. Sci. 13, 1199–1207 (2020).
  * [2] J.A.M. Borghans, R.J. de Boer, L.A. Segel: Extending the quasi-steady state approximation by changing variables. Bull. Math. Biol. 58, 43–63 (1996).
  * [3] G.E. Briggs, J.B.S. Haldane: A note on the kinetics of enzyme action, Biochem. J. 19, 338–339 (1925).
  * [4] P. De Leenheer, D. Angeli, E.D. Sontag: Monotone chemical reaction networks. J. Math. Chem. 41, 295–314 (2007).
  * [5] J. Eilertsen, S. Schnell, S. Walcher: On the anti-quasi-steady-state conditions of enzyme kinetics. Math. Biosci. 350, 108870 (2022).
  * [6] J. Eilertsen, S. Schnell, S. Walcher: The unreasonable effectiveness of the total quasi-steady state approximation, and its limitations. J. Theor. Biol., in press (preprint available in bioRxiv, 2023.09.28.559798) (2023).
  * [7] A. Goeke, S. Walcher: A constructive approach to quasi-steady state reduction. J. Math. Chem. 52, 2596–2626 (2014).
  * [8] M.W. Hirsch: Systems of differential equations which are competitive or cooperative. I. Limit sets. SIAM J. Math. Anal. 13, 167–179 (1982).
  * [9] E. Kamke: Zur Theorie der Systeme gewöhnlicher Differentialgleichungen II. Acta Math. 58, 57–85 (1932).
  * [10] J. Keener, J. Sneyd: Mathematical physiology I: Cellular physiology, Second Ed. Springer-Verlag, New York (2009).
  * [11] H.P. Kasserra, K.J. Laidler: Transient-phase studies of a trypsin-catalyzed reaction. Can. J. Chem. 48: 1793–-1802 (1970).
  * [12] L. Michaelis, M. L. Menten: Die Kinetik der Invertinwirkung. Biochem. Z., 49, 333–369 (1913).
  * [13] M.G. Pedersen, A.M. Bersani: Introducing total substrates simplifies theoretical analysis at non-negligible enzyme concentrations: Pseudo first-order kinetics and the loss of zero-order ultrasensitivity. J. Math. Biol., 60, 267–283 (2010).
  * [14] G. Pettersson A generalized theoretical treatment of the transient-state kinetics of enzymic reaction systems far from equilibrium. Acta Chem. Scand. B 32, 437–-446 (1978).
  * [15] S. Schnell, C. Mendoza: The condition for pseudo-first-order kinetics in enzymatic reactions is independent of the initial enzyme concentration. Biophys. Chem. 107: 165-–174 (2004).
  * [16] L.A. Segel, M. Slemrod: The quasi-steady-state assumption: A case study in perturbation. SIAM Review 31, 446–477 (1989).
  * [17] F. Silicio, M.D. Peterson: Ratio errors in pseudo first order reactions. J. Chem. Educ. 38: 576-–577 (1961).
  * [18] H.L. Smith: Monotone Dynamical Systems. Am. Math. Soc., Providence (1995).
  * [19] S. Sternberg: On the structure of local homeomorphisms of euclidean n-space. II. Amer. J. Math. 80, 623–631 (1958).
  * [20] N.A.T. Vu, Y.M. Song, Q. Tran, H.-Y. Yun, S.K. Kim, J-w. Chae and J.K. Kim: Beyond the Michaelis-–Menten: Accurate prediction of drug interactions through Cytochrome P450 3A4 induction. Clin. Pharmacol. Ther. 113, 1048–1057 (2023).


# Classification of Hyperspectral Images Using SVM with Shape-adaptive
Reconstruction and Smoothed Total Variation

Ruoning Li  Department of Mathematics, City University of Hong Kong  Kangning
Cui  Department of Mathematics, City University of Hong Kong  Hong Kong Centre
for Cerebro-Cardiovascular Health Engineering  Raymond H. Chan111Corresponding
Author<EMAIL_ADDRESS>Department of Mathematics, City University
of Hong Kong  Hong Kong Centre for Cerebro-Cardiovascular Health Engineering
Robert J. Plemmons  Departments of Mathematics and Computer Science, Wake
Forest University

###### Abstract

In this work, a novel algorithm called SVM with Shape-adaptive Reconstruction
and Smoothed Total Variation (SaR-SVM-STV) is introduced to classify
hyperspectral images, which makes full use of spatial and spectral
information. The Shape-adaptive Reconstruction (SaR) is introduced to
preprocess each pixel based on the Pearson Correlation between pixels in its
shape-adaptive (SA) region. Support Vector Machines (SVMs) are trained to
estimate the pixel-wise probability maps of each class. Then the Smoothed
Total Variation (STV) model is applied to denoise and generate the final
classification map. Experiments show that SaR-SVM-STV outperforms the SVM-STV
method with a few training labels, demonstrating the significance of
reconstructing hyperspectral images before classification.

Index Terms: Hyperspectral Image Classification, Image Reconstruction, Support
Vector Machines,
Smoothed Total Variation, Semi-supervised Learning

##  1 Introduction

Hyperspectral Images (HSIs) often have hundreds of electromagnetic reflectance
bands collected by aircraft or satellites that contain discriminative
information. HSIs can be represented by $X\in\mathbb{R}^{M\times N\times B}$,
where $M$ and $N$ represent the spatial size and $B$ represents the number of
spectral bands [1]. With the enhancement of the spectral resolution of HSIs,
machine learning approaches can be applied to investigate HSIs in land cover
identification, change recognition, and mineral mapping, etc. [2, 3, 4, 5].
The ground truth annotations of HSIs are expensive, which motivates us to take
less randomly picked labels for training [2].

In this work, SVM with Shape-adaptive Reconstruction and Smoothed Total
Variation (SaR-SVM-STV) is introduced to classify highly-mixed HSIs in a semi-
supervised way. In our previous work, SVM with Smoothed Total Variation (SVM-
STV) incorporated a convex variant of the Mumford-Shah denoising algorithm
with $\nu$-SVM Classifier to utilize spatial and spectral information [2]. It
has been shown that SVM-STV was highly competitive compared with state-of-the-
art algorithms, illustrating the effectiveness of incorporating spatial
information into classical pixel-wise algorithms. In this work, spatial
information is further utilized by a preprocessing step named Shape-adaptive
Reconstruction (SaR). The number of labeled pixels of each class needed to
achieve satisfying performance is less than 30.

This paper is organized as follows. Section 2 introduces the related works and
the motivation of this work. In Section 3, the workflow of the SaR-SVM-STV
Algorithm is introduced. Section 4 shows the numerical results of SaR-SVM-STV
compared with $\nu$-SVM and SVM-STV on two benchmark HSI datasets. Section 5
concludes the experiments and discusses the planned future work.

##  2 Related Works

Semi-supervised methods are desirable and computational feasible since they
only require a few labeled data for training. The classical semi-supervised
pixel-wise algorithms, such as SVMs, performed well for HSIs with small
computational cost [6]. A publicly available toolbox called LIBSVM collects
the most common SVM models that can be instantly utilized, promoting the
application of SVMs [7].

In general, pixel-wise classification methods can be enhanced by analyzing the
spatial dependency of nearby pixels, due to the common theoretical-based
assumption that the pixels are locally homogeneous in HSIs [8]. In [9],
spatial-spectral features are merged to train the SVM Classifier with Mercer’s
kernels. $K$-means and Principal Component Analysis (PCA) are used in [10] to
extract features, which are then fed into SVM with RBF kernels. The
performance of the above two methods grows significantly compared to pixel-
wise classifiers. The Smoothed Total Variation (STV) model has been used for
denoising and segmentation with superb performance [2, 11, 12]. The SVM-STV in
[2] firstly trains a $\nu$-SVM Classifier to estimate the probability maps of
each class, and then applies the STV model with $L_{1}$ and $L_{2}$
regularizers on the generated 3D probability tensor.

It is customary to employ spatial windows of pixels to extract features in
HSIs. Adopting square-shaped fixed-size windows is a straightforward solution.
A joint sparsity representation (JSR) proposed in [13] takes full advantage of
the spatial information by constructing JSR on local square patches. However,
pixels in a square window, particularly those on the margins or corners, may
not have similar spectra as the center pixel. Nevertheless, fixed-size windows
will fail when both spatially homogeneous and heterogeneous regions exist.
Nested sliding windows are designed in [14] to find the best nested square
window that maximizes the mean of the Pearson Correlation between window
center and other pixels, and then replace the spectra of center pixels with
the product of weighted Pearson Correlation vectors and spectra of pixels in
the best nested windows. In [15], the shape-adaptive (SA) regions of pixels
are learned by constructing local smooth irregular polygons in the first
principal component. The local polynomial approximation filtering (LPA) and
the intersection of confidence intervals rule (ICI) are applied to determine
the size and shape of the polygons [16]. Then a shape-adaptive JSR is applied
to explore the spatial and spectral information properly within SA regions
[15].

In this work, instead of using JSR, we reconstruct the spectrum of each pixel
using its SA region. The main difference between the SaR-SVM-STV proposed in
this work and SVM-STV is that the shape-adaptive reconstruction (SaR) is
applied to further utilize the spatial features of HSIs before training the
$\nu$-SVM. That is, we reconstruct the HSIs by SaR and reduce the dimension of
the reconstructed data by PCA, and then perform SVM-STV on the PCA reduced
data. The workflow of the SaR-SVM-STV model is summarized in Algorithm 1 and
the following section.

##  3 Methodology

Input: HSI, $\Upsilon$ (training set),

$\nu$ (tolerance of training error),

$\gamma$ (scaling factor of the RBF kernel),

$\beta_{1},\beta_{2}$ (regularization parameters)

Output: $\hat{Y}$ (classification map)

1 Reconstructing the spectrum of each pixel with its Shape-adaptive Regions:
$S_{c}=S_{nbs}\hat{p}(x_{c})$, and get the reconstructed data $\tilde{X}$;

2 Reducing the dimension of data by PCA: $X_{r}=$PCA$(\tilde{X})$;

3 Estimating the probability tensor $V$ by using $\nu$-SVMs on the HSI data
with parameters $\nu$ and $\gamma$: $V=\nu$-SVMs$(X_{r},\nu,\gamma)$;

4 Applying STV model on the probability tensor $V$ to get the denoised tensor
$U$ for each class:
$\min\limits_{{\rm\mathbf{u}}_{k}}\frac{1}{2}\|{\rm\mathbf{u}}_{k}-{\rm\mathbf{v}}_{k}\|_{2}^{2}+\beta_{1}\|{\rm\mathbf{\nabla
u}}_{k}\|_{1}+\frac{\beta_{2}}{2}\|{\rm\mathbf{\nabla u}}_{k}\|_{2}^{2},{\rm
s.t.,}\ {\rm\mathbf{u}}_{k}|_{\Upsilon}={\rm\mathbf{v}}_{k}|_{\Upsilon}$;

5 The classification map $\hat{Y}$ is obtained by
$\hat{Y}_{i,j}=$MaxIndex$(U_{i,j,:})$.

Algorithm 1 SVM with Shape-adaptive Reconstruction and Smoothed Total
Variation (SaR-SVM-STV)

###  3.1 Shape-adaptive Reconstruction

The Shape-adaptive Reconstruction (SaR) is introduced as an efficient tool to
further utilize spatial features as a preprocessing step. The shape-adaptive
(SA) region is generated to pick surrounding pixels that are the closest to
the target pixel by the LPA denoising and ICI rule [15, 16]. The Pearson
Correlation coefficients ${\rm Corr}(x_{c},x_{i})=\frac{{\rm
Cov}(x_{c},x_{i})}{\sqrt{{\rm Var}(x_{c})\cdot{\rm Var}(x_{i})}}$ are computed
to estimate the correlation between the center pixel $x_{c}$ and the other
$n-1$ pixels $x_{i}$, where $n$ represents the total number of pixels in
$x_{c}$’s SA region. Then a vector $p(x_{c})\in\mathbb{R}^{n\times 1}$ is
formed with the $i$-th entry Corr$(x_{c},x_{i})$ when $i=1,2,\cdots,n-1$, and
the $n$-th entry 1. The correlation coefficients are then rescaled as a weight
vector $\hat{p}(x_{c})=\frac{p(x_{c})}{{\rm Sum}(p(x_{c}))}$, and the spectrum
of the center pixel is reconstructed by $S_{c}=S_{nbs}\hat{p}(x_{c})$, where
$S_{nbs}\in\mathbb{R}^{B\times n}$ represents the spectra collection of the
$n$ pixels in the SA region of pixel $x_{c}$ and the $n$-th column represents
the spectrum of $x_{c}$. Thus, pixels inside SA regions will not contribute to
the spectrum of center pixel uniformly, with highly-correlated pixels
dominating the spectrum of center pixel as they are associated with large
weights in $\hat{p}$. PCA dimension reduction is performed to further reduce
the noise in HSIs and speed up the training of our proposed model. Experiments
show that using the data before or after PCA produces similar classification
results.

###  3.2 SVM with Smoothed Total Variation

The $\nu$-SVMs are trained by the randomly chosen labels $\Upsilon$ with one-
against-one strategy and five-fold cross-validation, and predict the 3D
probability map $V\in\mathbb{R}^{M\times N\times K}$. The $k$-th channel
$\mathbf{v}_{k}=V_{:,:,k},k=1,...,K$ represents the probability of the pixels
belonging to the class $k$ [2]. The toolbox of LIBSVM library is used to
implement $\nu$-SVM [7]. Then the STV model is adopted to reduce the noise in
the probability tensor $V$ by applying the alternating direction method of
multipliers (ADMM) with the convergence guaranteed [11, 12, 17]. The STV model
is shown in Step 4 of Algorithm 1. The $\mathbf{u}_{k}$ represents the
smoothed probability map of class $k$. Meanwhile, the labels of pixels in the
training set $\Upsilon$ should remain unchanged. The smoothed 3D tensor $U$ is
obtained after the STV noise reduction, and the final classification map
$\hat{Y}$ is then produced by taking the index of the class that has the
maximum probability of each pixel [2].

##  4 Numerical Results

Figure 1: The performance of $\nu$-SVM, SVM-STV, and SaR-SVM-STV on Pavia U
and Indian Pines. The two columns present the performance of Pavia U and
Indian Pines, and the three rows show the overall accuracy, the average
accuracy, and the kappa coefficient of the three algorithms we compare
against, respectively.

In this section, we show the results of SaR-SVM-STV algorithm compared with
$\nu$-SVM and SVM-STV on two publicly available datasets: Indian Pines and
University of Pavia (Pavia U). Details of the datasets are summarized in Table
1.

Table 1: Summary of the two HSI datasets

| Type | Size | No. Classes | Sensor | Location
---|---|---|---|---|---
Indian Pines | Forest | 145$\times$145$\times$200 | 16 | AVIRIS | Indiana, US
Pavia U | Urban | 610$\times$340$\times$103 | 9 | ROSIS | Pavia, Italy

The PCA dimension reduction is used to keep $99.9\%$ of the variance of the
reconstructed data. The best $\nu$ and $\gamma$ are selected by the cross-
validation of $\nu$-SVM. The best smoothing parameters $\beta_{1}$ and
$\beta_{2}$ in the STV model are tuned on the probability tensor $V$ that
maximize the sum of the overall accuracy (OA), the average accuracy (AA), and
the Kappa coefficient in 10 trials. For Pavia U, $\beta_{1}=0.2$, and
$\beta_{2}=1$. For Indian Pines, $\beta_{1}=0.2$, and $\beta_{2}=4$. The
source code used in this work can be found at: github.com/ckn3/SA-Recon.

Fig. 1 shows the tendency of performance while increasing the labeled pixels
randomly picked from each class. The three metrics improve substantially by
incorporating more training labels, and the SaR-SVM-STV model outperforms the
comparison methods in all cases. Fig. 2 shows relationship between the number
of labeled pixels involved for training and the count of misclassifications on
Indian Pines using SaR-SVM-STV. As the training size increases, the portion of
misclassifications reduces monotonically.

Fig. 3 indicates the efficacy of incorporating spatial information. By
comparing Figs. 3 (c) and (d), it can be found that although the spatial noise
is reduced, the SVM-STV model is misled by misclassified pixels in the
predicted map of $\nu$-SVM when only a few labels are used. Applying the SaR
and PCA algorithms smooth the HSIs before classification, and thus reduce the
number of misclassified pixels in the $\nu$-SVM step. The STV model can
further tick out the misclassified pixels by enforcing spatial connectivity in
probability maps as a postprocessing step, see Fig. 3 (e).

##  5 Conclusions

In [2], SVM-STV has been shown to outperform many other methods. In this work,
we show that adding a preprocessing step, i.e., the Shape-adaptive
Reconstruction that reconstruct HSIs using spatial information, can further
enhance the performance. The SaR preprocessing and STV postprocessing steps
can be easily merged with other semi-supervised or unsupervised algorithms
that do not utilize the spatial information of HSIs. The work planned in the
near future would be combining the image processing techniques with
unsupervised clustering methods [18, 19], and developing the active learning
extension where a few carefully-selected annotations are used to train the
model [20, 21].

Figure 2: Misclassification heatmaps of Indian Pines using SaR-SVM-STV with
different numbers of labeled pixels. Figures (a-b) show the ground truth (each
color represents a class) and false-color visualization of Indian Pines.
Figures (c-f) show the heatmaps when 30, 20, 10, and 5 labeled pixels per
class are used. The colorbar represents the count of misclassifications in 10
trials. Figure 3: Misclassification heatmaps of Pavia U by comparing the three
methods when 10 labeled pixels per class are used for training. Figures (a-b)
show the ground truth (each color represents a class) and false-color
visualization of Pavia U. Figures (c-e) show the heatmaps of $\nu$-SVM, SVM-
STV, and SaR-SVM-STV. The colorbar represents the count of misclassifications
in 10 trials.

## References

  * [1] M.. Eismann  “Hyperspectral remote sensing”  SPIE, 2012
  * [2] R.. Chan, K.. Kan, M. Nikolova and R.. Plemmons  “A two-stage method for spectral–spatial classification of hyperspectral images”  In _J. Math. Imaging Vis._ 62.6 Springer, 2020, pp. 790–807
  * [3] K. Cui and R.. Plemmons  “Unsupervised Classification of AVIRIS-NG Hyperspectral Images”  In _Proc IEEE Workshop Hyperspectral Image Signal Process Evol Remote Sens_ , 2021, pp. 1–5  IEEE
  * [4] S. Camalan  “Change Detection of Amazonian Alluvial Gold Mining Using Deep Learning and Sentinel-2 Imagery”  In _Remote Sens._ 14.7, 2022  DOI: 10.3390/rs14071746
  * [5] M.. Khan et al.  “Modern Trends in Hyperspectral Image Analysis: A Review”  In _IEEE Access_ 6, 2018, pp. 14118–14129  DOI: 10.1109/ACCESS.2018.2812999
  * [6] F. Melgani and L. Bruzzone  “Classification of hyperspectral remote sensing images with support vector machines”  In _IEEE Trans. Geosci. Remote Sens._ 42.8 IEEE, 2004, pp. 1778–1790
  * [7] C. Chang and C. Lin  “LIBSVM: a library for support vector machines”  In _ACM Trans. Intell. Syst. Technol._ 2.3 Acm New York, NY, USA, 2011, pp. 1–27
  * [8] J. Liu et al.  “Spatial-spectral kernel sparse representation for hyperspectral image classification”  In _IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens._ 6.6 IEEE, 2013, pp. 2462–2471
  * [9] G. Camps-Valls et al.  “Composite kernels for hyperspectral image classification”  In _IEEE Geosci. Remote. Sens. Lett._ 3.1 IEEE, 2006, pp. 93–97
  * [10] L. Shu, K. McIsaac and G.. Osinski  “Learning spatial–spectral features for hyperspectral image classification”  In _IEEE Trans. Geosci. Remote Sens._ 56.9 IEEE, 2018, pp. 5138–5147
  * [11] X. Cai, R. Chan and T. Zeng  “A two-stage image segmentation method using a convex variant of the Mumford–Shah model and thresholding”  In _SIAM J. Imaging Sci._ 6.1 SIAM, 2013, pp. 368–390
  * [12] X. Cai, R. Chan, M. Nikolova and T. Zeng  “A three-stage approach for segmenting degraded color images: Smoothing, lifting and thresholding (SLaT)”  In _J. Sci. Comput._ 72.3 Springer, 2017, pp. 1313–1332
  * [13] H. Zhang, J. Li, Y. Huang and L. Zhang  “A Nonlocal Weighted Joint Sparse Representation Classification Method for Hyperspectral Imagery”  In _IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens._ 7.6, 2014, pp. 2056–2065  DOI: 10.1109/JSTARS.2013.2264720
  * [14] J. Ren et al.  “An SVM-Based Nested Sliding Window Approach for Spectral–Spatial Classification of Hyperspectral Images”  In _Remote Sens._ 13.1 Multidisciplinary Digital Publishing Institute, 2021, pp. 114
  * [15] W. Fu et al.  “Hyperspectral Image Classification Via Shape-Adaptive Joint Sparse Representation”  In _IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens._ 9.2, 2016, pp. 556–567  DOI: 10.1109/JSTARS.2015.2477364
  * [16] A. Foi, V. Katkovnik and K. Egiazarian  “Pointwise shape-adaptive DCT for high-quality denoising and deblocking of grayscale and color images”  In _IEEE Trans. Image Process._ 16.5 IEEE, 2007, pp. 1395–1411
  * [17] S. Boyd, N. Parikh and E. Chu  “Distributed optimization and statistical learning via the alternating direction method of multipliers”  Now Publishers Inc, 2011
  * [18] S.. Polk, K. Cui, R.. Plemmons and J.. Murphy  “Diffusion and Volume Maximization-Based Clustering of Highly Mixed Hyperspectral Images”  In _arXiv preprint arXiv:2203.09992_ , 2022
  * [19] J.. Murphy and S.. Polk  “A multiscale environment for learning by diffusion”  In _Appl. Comput. Harmon. Anal._ 57 Elsevier, 2022, pp. 58–100
  * [20] J.. Murphy and M. Maggioni  “Unsupervised clustering and active learning of hyperspectral images with nonlinear diffusion”  In _IEEE Trans. Geosci. Remote Sens._ 57.3 IEEE, 2018, pp. 1829–1845
  * [21] S.. Polk, K. Cui, R.. Plemmons and Murphy J. M.  “Active Diffusion and VCA-Assisted Image Segmentation of Hyperspectral Images”  In _arXiv preprint arXiv:2204.06298_ , 2022


# Chiral magnetic phases in Moire bilayers of magnetic dipoles

Ignacio Tapia  Xavier Cazor  Paula Mellado  School of Engineering and
Sciences, Universidad Adolfo Ibáñez, Santiago, Chile

###### Abstract

In magnetic insulators, the sense of rotation of the magnetization is
associated with novel phases of matter and exotic transport phenomena. Aimed
to find new sources of chiral magnetism rooted in intrinsic fields and
geometry, we study twisted square bilayers of magnetic dipoles with easy plane
anisotropy. For no twist, each lattice settles in the antiferromagnetic zig-
zag magnetic state and orders antiferromagnetically to the other layer. The
moire patterns that result from the mutual rotation of the two square lattices
influence such zig-zag order, giving rise to several magnetic phases that
depict non-collinear magnetic textures with chiral motifs that break both time
and inversion symmetry. For certain moire angles, helical and toroidal
magnetic orders arise. These novel phases can be further manipulated by
changing the vertical distance between layers. We show that the dipolar
interlayer interaction originates an emergent twist-dependent chiral magnetic
field orthogonal to the direction of the zig-zag chains, which is responsible
for the internal torques conjugated to the toroidal orders.

##  I Introduction

Chiral magnetic textures are promising building blocks for storing information
and offer new venues for the transport of energy in magnetic insulators Casher
and Susskind (1974); Tomita et al. (2018); Thiaville and Fert (1992). The
onset of such chiral textures is usually conditioned to the presence of
antisymmetric spin interactions associated with spin-orbit effects such as the
Dzyaloshinskii-Moriya (DM) interaction Dzyaloshinskii (1960, 1958); Heyderman
and Stamps (2013); Togawa et al. (2016) or dipolar interactions. Indeed, for
lattices of magnetic dipoles featuring out-of-plane anisotropy, the dipolar
coupling among all point dipoles gives rise to antisymmetric interactions
manifesting in chiral internal magnetic fields Mellado et al. (2023); Paula
and Tapia (2023); Yu and Bauer (2021); Ray and Das (2021); Shindou et al.
(2013); Lucassen et al. (2019). Among magnetic chiral orders, the study of
skyrmions Nagaosa and Tokura (2013); Mohylna et al. (2022); Camosi et al.
(2018) has led to a considerable amount of research motivating the inclusion
of an elegant mathematical formalism in the field of topological defects,
which is suitable to characterize a plethora of helical and chiral structures
Schulz et al. (2012); Shaban et al. (2023). While finding materials that
realize skyrmionic lattices is an active field of research, finding fields
that can efficiently manipulate skyrmions and other topological textures as
domain walls is crucial to moving the field to the next stage. More recently,
ferrotoroidal phases Dong et al. (2015); Spaldin et al. (2008); Zimmermann et
al. (2014); Kuprov et al. (2022) have attracted much interest because of their
putative role in the manifestation of the magnetoelectric effect Astrov
(1960); Date et al. (1961); Siratori et al. (1992); Ederer and Spaldin (2007);
Mellado et al. (2020). Ferrotoroidicity has been defined as a fourth form of
ferroic order and consists of a vortex-like magnetic state with zero
magnetization and a spontaneously occurring toroidal moment. Because the
toroidic order has zero magnetization, the main challenge is finding a
conjugate field that couples to the toroidic ferroic order. Recent experiments
with permalloy nanoislands interacting via dipolar coupling indeed manifested
the ferrotoroidic phase and have achieved local control of the
ferrotoroidicity by scanning a magnetic tip, thus widening the possibilities
of control on chiral phases Lehmann and Lehmann (2022). In the same direction,
the ideas of “twistronics” where the control over electric phases is achieved
by the relative rotation between 2D layers of van der Waals materials Klebl et
al. (2022); Xiao et al. (2020), which include graphene-based structures Andrei
et al. (2021); He et al. (2021), transition-metal dichalcogenides Tran et al.
(2019); Xiao et al. (2021), 2D magnets Li and Cheng (2020); Begum Popy et al.
(2022), and thin films have motivated a related paradigm of twisted bilayer
structures in magnetic systems. This approach does not rely on the engineering
of correlated flat bands but can produce interesting new phases by combining
known nontrivial properties of constituent monolayers Heyderman and Stamps
(2013). While in heterostructures, the spin-wave control is accomplished by
tuning the individual layer thicknesses and external fields, in magnetic
insulators, the relative twist of bilayers with antiferromagnetic and
ferromagnetic spin couplings was predicted to give rise to noncollinear
magnetic textures for small twist angles between layers Hejazi et al. (2020);
Yang et al. (2023); Li and Cheng (2020); Shimizu et al. (2021); Claassen et
al. (2022); Luo et al. (2022). Some of these predictions have been tested
recently in experiments Song et al. (2021).

Summary of results: motivated by the search of novel magnetic chiral phases
and based on the ideas of twistronics, here we study the dynamics of magnetic
dipoles placed at bilayers of square lattices and show that moire arrays of
magnetic dipoles give rise to novel chiral magnetic phases absent in the
untwisted counterparts. We find that depending on the twist angle $\theta$ and
the distance between the layers $h$, the system mutates from an
antiferromagnetic zig-zag pattern at zero twists to a set of chiral magnetic
configurations which go from periodic arrays of well-defined magnetic
antiferromagnetic toroidal patterns to disordered twisted phases passing
through helical periodic configurations. Toroidal correlations are easily
tuned by varying $\theta$ and $h$. We explain the chiral dipolar orders in
terms of a chiral internal magnetic field which arises from the relative twist
among the dipolar layers.

The rest of the paper is organized as follows. Section II describes square and
moire lattices, the system of dipoles, and the inter and intra-layer dipolar
interactions that arise. The different magnetic phases that emerge when the
twist angle and the distance between layers is varied are described in
sections III and IV. In section V we study the chiral magnetic fields that
arise due to the twist and which are responsible for the chiral phases
reported in section III. Section VI is devoted to concluding remarks. Appendix
A shows details of molecular dynamics simulations used to solve the dynamics
of each interacting dipole in the bilayer system, while Appendix B shows
details of calculations shown in V.

Figure 1: (a) The zig-zag order (ZZ) of the bilayer system for $\theta=0$.
(b-d) Moire patterns obtained for twist angles (b)
$\theta_{6,7}=8.79^{\circ}$, (c) $\theta_{2,3}=22.61^{\circ}$ and (d)
$\theta_{1,2}=36.86^{\circ}$. Black long arrows show the moire lattice
primitive vectors $v_{1}$ and $v_{2}$ in each case. (f-j) Real space dipolar
configurations obtained from the numerical solution of the equations of motion
for systems with (e) $\theta_{2,3}$ and $h=0.1$ featuring phase Helical I
(HI), (f) $\theta_{2,3}$ and $h=0.1$ featuring phase Toroidal I (TI), (g)
$\theta_{1,2}$ and $h=0.1$ featuring phase TI, (h) $\theta_{6,7}$ and $h=1$
featuring phase Helical II (HII), (i) $\theta_{2,3}$ and $h=1.0$ featuring
phase Toroidal II (TII) and (j) $\theta_{1,2}$ and $h=0.9$ featuring phase
TII. In all cases small blue arrows illustrate dipoles in layer 1 while red
small arrows dipoles in lattice 2. Yellow squares indicate the moire lattice
unit cell, while black squares indicate the magnetic unit cell. (k) Phase
Diagram in terms of the twist angle $\theta$ and the dimensionless interlayer
distance $h$. Figure 2: (a) Density plot of the magnetization of a dipolar
layer featuring the ZZ order.(b-g) Density plots of the $z$-component of
toroidal moment in twisted bilayers featuring phases (b) Helicoidal I:
$\theta_{6\,7}=8.78^{\circ}$, $h=0.1$. (c) Helicoidal II:
$\theta_{6\,7}=8.78^{\circ}$, $h=1.4$. (d) Toroidal I:
$\theta_{2\,3}=22.62^{\circ}$, $h=0.1$. (e) Toroidal II:
$\theta_{2\,3}=22.62^{\circ}$, $h=1.4$. (f) Toroidal I:
$\theta_{1\,2}=36.87^{\circ}$, $h=0.1$. (g) Toroidal II:
$\theta_{1\,2}=36.87^{\circ}$, $h=0.9$. Figure 3: Density plot of the
structure factor of the order parameters: (a) the magnetization for a single
layer with zig-zag order , and the toroidal moment for Toroidal I phase
obtained for (b) ($\theta_{23}=22.61^{\circ}$, $h=0.1$) and (c)
($\theta_{12}=36.87^{\circ}$, $h=0.1$). The dipolar energy obtained for
simulations with $h=0.1$ as a function of $\theta$. The solid line shows the
dipolar energy of the zig-zag configuration for different angles, while the
dots show the energies obtained from the relaxation dynamics. The color of the
dots indicates the phase obtained in the simulation: yellow for Helicoidal I
and red for Toroidal I.

##  II Model

The bilayer system consists of two identical square lattices with lattice
constant $a$ that extend in the $\hat{x}-\hat{y}$ plane and are stacked at a
distance $h$ along the orthogonal $\hat{z}$ direction. At the $n=n_{s}^{2}$
vertices of each square layer of side $L=an_{s}$, $n$ magnetic dipoles are
placed. They interact through the full long-range dipolar interaction with all
other dipoles of the bilayer system:

$\displaystyle
U_{d}=\frac{g}{2}\sum_{\mu,\nu}^{2}\sum_{i,j}^{n}\frac{\hat{m}_{i}^{\mu}\cdot\hat{m}_{j}^{\nu}-3(\hat{m}_{i}^{\mu}\cdot\hat{r}_{ij}^{\mu,\nu})(\hat{m}_{j}^{\nu}\cdot\hat{r}_{ij}^{\mu,\nu})}{|{\bm{r}}_{i}^{\mu}-{\bm{r}}_{j}^{\nu}|^{3}}$
(1)

where
$\hat{\bm{r}}_{ij}^{\mu\nu}=\frac{({\bm{r}}_{i}^{\mu}-{\bm{r}}_{j}^{\nu})}{|{\bm{r}}_{i}^{\mu}-{\bm{r}}_{j}^{\nu}|}$
denotes the unit vector joining dipole $\hat{\bm{m}}_{i}^{\mu}$ at site $i$
and layer $\mu$ and dipole $\hat{\bm{m}}_{i}^{\nu}$ at site $j$ and layer
$\nu$. $g=\frac{u_{0}m_{0}^{2}}{4\pi a^{3}}$ sets an energy scale Mellado et
al. (2023) and contains the physical parameters of the system, such as $a$,
the lattice constant along $x$ and $y$ directions, $u_{0}$, the magnetic
permeability, and $m_{0}$, the intensity of the magnetic moments. Each layer
has the full symmetry of the square Bravais lattice. They are set apart along
$\hat{z}$ by distance $z=ha$ where $h$ is a positive real number in the range
$h\in(0.1,1)$. Dipoles can rotate in the $\hat{x}-\hat{y}$ plane respect to a
local axis fixed at their vertices, and have unit vector
$\hat{m}_{i}^{\mu}=(m_{i,x}^{\mu},m_{i,y}^{\mu})=(\cos\alpha_{i}^{\mu},\sin\alpha_{i}^{\mu})$
where $\alpha_{i}^{\mu}$ defines the polar rotation angle of dipole i at layer
$\mu$. Hereafter the magnetic moments are normalized by $m_{0}$.

Global rotations between layer are performed with respect to axis $\hat{z}$
located at the center of the two layers. The twisted bilayer system forms a
crystal only for a discrete set of commensurate rotation angles. We focus on
commensurate twist geometries that are amenable to numerical study using the
standard Bloch representation. Hence we rotate the layers around specific
angles that respond to the equation Kariyado and Vishwanath (2019):
$\mathrm{\theta}_{p,q}=2\arctan(p/q)$, where p and q are natural numbers, and
the angle is in radians. These types of rotations generate systems called
moire geometries, where the patterns depend on $\mathrm{\theta}_{p,q}$
Kariyado and Vishwanath (2019). We note that any commensurate twist can be
described by an integer valued ‘twist’ vector $v=(p,q)$ as shown in
Figs.1(b-d). The corresponding moire unit cell comprises of
$\mathrm{2}(p^{2}+q^{2})$ sites. We rotate the two perfectly aligned square
lattices in opposite directions by $\mathrm{\theta}_{p,q}/2$, which leads to
the sites originally at the locations $v$ and $\bar{v}=(-p,q)$ to lie atop
each other. The unit cells for a selection of commensurate twist angles are
illustrated in Fig.1(b-d).

Writing the unit vectors characterizing each layer as ${\bf a}_{1}$ and ${\bf
a}_{2}$, the moire pattern induced by twist has a periodicity of ${\bf v}_{1}$
and ${\bf v}_{2}$, which is explicitly written as

${\bf v}_{i}=v_{i}\hat{z}\times{\bf a}_{i}/(\sin\frac{\theta}{2})$ (2)

Hermann (2012). The reciprocal lattice vectors of the first Brillouin zone of
the moire lattice in each case are given by ${\bf
G}_{i}=\frac{2\pi}{{\bf{v}_{i}}}$.

Next we consider the dipolar interaction between all dipoles of the bilayer at
zero twist. By separating $\mathrm{U}_{d}$ in intra and inter-layer
interactions $U_{d}^{\mu\mu}$ and $U_{d}^{\mu\nu}$, and identifying
contributions between pair of dipoles with director vector along the square
lattice symmetry axes $\hat{n}_{i}\in(\hat{x},\hat{y})$ and along any other
direction, it is easy to see that dipoles prefer to align ferromagnetically
with neighbours located along one of the $\hat{n}_{i}$ axes and
antiferromagnetically with dipoles along the other orthogonal axis. Given the
symmetry of the square lattice, the ground state magnetic configuration is
already determined by dipolar interactions at the nearest neighbors level.
Without loss of generality the square lattice can be seen as a collection of
chains of dipoles along the $\hat{n_{1}}$ direction that extend parallel along
the orthogonal direction $\hat{n_{2}}$. We can separate dipoles along the axis
$\hat{n_{2}}$ in two sublattices $s_{a}\in(s_{1},s_{2})$ which correspond to
the sets of odd and even positions along $\hat{n_{2}}$. Thus regardless the
layer, the dipolar energy of the system can be separated in the contribution
from dipoles interacting in the same sublattice and dipoles interacting in
different sublattices:
$U_{d}=\sum_{s_{a}}U_{d}^{(s_{a})}+U_{d}^{(s_{a},s_{b})}$. It is easy to
verify that parallel dipoles minimize their energy when they point in opposite
senses while collinear dipoles prefer to point in the same sense. If we choose
all dipoles to be parallel to $\hat{n_{1}}$, then in units of g the dipolar
interaction between dipoles in the same and different sublattices becomes
respectively $U_{d}^{(s_{a})}=\sum_{s_{a}}\sum_{i,j\in
s_{a}}\frac{\hat{m_{i}}^{s_{a}}\cdot\hat{m_{j}}^{s_{a}}}{(r_{ij})^{3}}(1-3(\frac{(r_{ij}^{n_{1}})^{2}}{(r_{ij})^{2}}))$
and $U_{d}^{(s_{a},s_{b})}=\sum_{s_{a},s_{b}}\sum_{i\in s_{a}j\in
s_{b}}\frac{\hat{m_{i}}^{s_{a}}\cdot\hat{m_{j}}^{s_{a}}}{(r_{ij})^{3}}(1-3(\frac{(r_{ij}^{n_{1}})^{2}}{(r_{ij})^{2}}))$.
A dipole in $s_{1}$, has two nearest neighbors in $s_{1}$ along $\hat{n}_{1}$
and two in $s_{2}$ along $\hat{n}_{2}$, and at the nearest neighbour level we
get
$U_{d}^{1nn}=\frac{\hat{m_{i}}^{s_{1}}\cdot\hat{m_{j}}^{s_{1}}}{(1+h^{2})^{3}}(1-\frac{3}{(1+h^{2})^{2}})+\frac{\hat{m_{i}}^{s_{1}}\cdot\hat{m_{j}}^{s_{2}}}{(1+h^{2})^{3}}$.
The energy is minimized when dipoles in the same sublattice are ferromagnetic
and dipoles in different sublattices are antiferromagnetically arranged. At
second nearest neighbours level, interactions happen only between dipoles in
different sublattices and though they favour the ferromagnetic interactions
between $s_{1}$ and $s_{2}$, these are proportional to
$\sim\frac{1}{(2+h^{2})^{3}}$ and therefore are overcomed by the
antiferromagnetic first neighbour interaction. Third nearest neighbours
interactions occur only between dipoles in the same sublattice and favour a
ferromagnetic coupling between dipoles consistent with the first neighbour
interaction.

The addition of another layer at zero twist angle further stabilizes the zig-
zag order in each layer. The relative orientation of dipoles in the two layers
at a distance $h$ is antiferromagnetic, which indeed corresponds to the most
favourable magnetic configuration between two parallel magnetic dipoles.
Magnetization and staggered magnetization in the bilayer system are defined
$\rm{\bf M}=\frac{1}{n}\sum_{i}^{n}{\bf m}_{i}$ and $\rm{\bf
N}=\frac{1}{n}\sum_{i}^{n}({\bf m}_{i}^{s_{1}}-{\bf m}_{i}^{s_{2}})$
respectively.

##  III Phase Diagram

To examine the possible equilibrium magnetic configurations of twisted dipolar
bilayers, we solved the equations of motion of each dipole when the lattices
are mutually rotated by $\theta/2$ in opposite senses and for a set of
interlayer distances $h$, using molecular dynamics simulations. Details are
presented in Appendix A. In the simulations, each dipole $i$ consists of a
magnetic bar with inertia moment $I$ and magnetic moment intensity $m_{0}$,
which rotates with angle $\alpha_{i}$ respect to its local axis of rotation.
The evolution of $\alpha_{i}$ is due to the torques $\mathcal{T}_{i}$ produced
by dipolar interaction of dipole $i$ with all other dipoles in the lattice.
The angular rotation has damping, product of the friction at the rotation axis
as shown by Eq.9.

The twist angle determines the moire pattern of the bilayer while $h$ changes
the strength of the interlayer interactions. Indeed, depending on $\theta$ and
$h$, dipoles settle into six magnetic phases. Real space dipolar
configurations obtained from the molecular dynamics simulations are shown in
Fig.1(a) and Fig.1(e-j). At $\theta=0$ and for all the range of $h$ the
bilayer settles in a zig-zag antiferromagnetic (ZZ) phase, in the
$\hat{x}-\hat{y}$ plane as shown in Fig.1(a). The antiferromagnetic zig-zag
phase is a colinear magnetic order along any of the axis $\hat{\bf
n}_{1},\hat{\bf n}_{2}\in(\hat{\bf x},\hat{\bf y})$, which doubles the unit
cell of the square lattice made of lattice vectors ${\bf a}_{1}=a(1,0)$ and
${\bf a}_{2}=a(0,1)$. As shown in Fig3(a) there are two possible order wave
vectors: the two X points at the centers of the edges of the square first
Brillouin zone (1BZ), which are half reciprocal vectors ${\bf
b}_{1}/2=(\frac{\pi}{a},0)$ and ${\bf b}_{2}/2=(0,\frac{\pi}{a})$. The ZZ
phase survives up to $\theta_{21,22}=\sim 2.7$.

Due to the $C_{4}$ symmetry of the square lattice, the maximum twist between
layers is reached at $\theta=45^{\circ}$. In this limit the ground state
magnetic state depends on the interlayer distance $h$. At $h\in(0.1,0.4)$
dipoles order in a chiral pattern denoted toroidal I (TI). TI is free of M and
N and is featured by dipoles arranged in periodic octagonal chiral loops and
stars as shown in Figs.1(f) and (g). The orientation of dipoles respect to the
center of each loop, determines its vorticity $\mathcal{V}$. The octagonal
loops in this phase have a vorticity which follows a zig-zag antiferromagnetic
pattern. In the range $h\in(0.5,0.9)$ the periodic dipolar loops disappear and
the system settles in a disordered state (D), until at $h=0.9$ it returns to a
canted zig-zag pattern with dipoles in both lattices a bit rotated respect to
the perfect collinear order of the untwisted bilayer. The previous scenario
remains for moire angles in the range $\theta_{p,q}\in(43,45)$.

At angles $\theta_{p,q}\in(39,31)$ a new phase arises for large values of $h$.
It is denoted Toroidal II (TII). Similar to TI it features chiral loops but
this time their vorticity follows an antiferromagnetic pattern as shown in
Fig.1(i) and Fig.1(j). The loops are intercalated by a square arrays of stars
which feature zero flux, and by other squared pattern of pairs of dipoles
arranged antiferromagnetically.

At moire angles $\theta_{p,q}\in(2,14)$ and for $h\in(0.1,0.4)$, dipoles
arrange in phase Helical I (HI) shown in Fig1(e). The magnetic configuration
is magnetization free and consists of a square pattern of a squared shaped
magnetic vortices that contain several dipoles and resemble in-plane
projections of spins in skyrmionic lattices. We describe this phase in more
detail in section IV. In this range of $\theta_{p,q}$ and for large h beyond
the D phase, the system orders once again and settles in phase Helical II
(HII) shown in Fig1(h). HII differs from HI in the size of the magnetic unit
cell. Larger interlayer distances are associated with smaller interlayer
interactions. In these cases the magnetic unit cell grows, decreasing the
order wavevector $q^{ord}$.

##  IV Antiferromagnetic toroidic orders.

In this section, we focus on helical and toroidal phases. Both break time and
inversion symmetries and depict chiral magnetic motifs made out of dipoles of
both layers. The main difference is the number of dipoles participating in the
pattern that repeats in the magnetic unit cell. Toroidal phases possess
periodic clusters that involve eight or fewer dipoles, while the periodic
motif in helical phases resembles skyrmion textures in the plane, with many
dipoles being part of a sizeable chiral cluster.

In order to characterize toroidal and helical chiral phases, we study their
ferrotoroidic order parameter T or toroidization Lehmann and Lehmann (2022),
which is derived from the toroidal moment $\bf t=\sum_{i}r_{i}\times m_{i}$,
where $r_{i}$ and $m_{i}$ are the position and magnetic moment of the i-th
dipole within the unit cell Spaldin et al. (2008). In the toroidal and helical
phases, toroidal moments $t_{j}$ belonging to magnetic unit cells $j$ can
align in particular ways to form a uniform toroidic state with a spontaneous
macroscopic toroidization, ${\bf T}=\sum{\bf t}_{j}$ directed along
$\pm\hat{z}$. The sense (sign) of the toroidal moment depends on the
handedness of the arrangement of magnetic moments in a magnetic unit cell.

The in-plane helical patterns realized by chiral canted dipoles in these
phases can be distinguished in terms of vorticity and helicity, which are used
to characterize skyrmion spin structuresNagaosa and Tokura (2013). Introducing
the polar coordinates $r=(r\cos\phi,r\sin\phi)$ at the center of each helical
motif, for a dipolar texture $m(r)=(\cos\alpha(\phi),\sin\alpha(\phi))$, the
vorticity is defined by the integer
$\mathcal{V}=[\alpha(\phi)]_{\phi=0}^{\phi=2\pi}/2\pi$ Nagaosa and Tokura
(2013). The helicity, $\gamma$ accounts for the phase appearing in
$\alpha(\phi)=m\phi+\gamma$ Nagaosa and Tokura (2013). In a toroidal phase
with all motifs having the same vorticity, $\gamma$ distinguishes octagonal
loops from octagonal stars like those shown in Fig.1(j). Vorticity becomes
helpful in distinguishing helical from toroidal phases. While toroidal phases
have $\mathcal{V}=1$, helical phases feature dipolar clusters with
$\mathcal{V}=-1$.

Figs.2(b-g) show real space density maps of the toroidal order parameter T in
moire dipolar bilayers featuring phases TI, TII, HI, and HII, where x and y
axes are in units of the square lattice constant. Blue and yellow stripes in
the magnetization density plot of Figs.2(a) clearly distinguish ZZ phases. In
toroidal phases, the periodic pattern of well-defined yellow and blue dots in
Figs.2(d-g) accounts for the antiferromagnetic toroidic order TI and for the
antiferromagnetic zig-zag pattern of toroidal loops realized at large
$\theta_{1,2}$ and $h\in(0.1,0.4)$. The helical order depicts concentric
squared stripes that alternate the sense of T as shown in Figs.2(b-c) .

The wave vector associated with zig-zag and AF orders can be identified in the
structure factor of T,

$S_{q}^{T}=\frac{1}{n}|\sum_{j}T_{j}e^{i{\bf q^{m}}\cdot{\bf r}_{j}}|^{2}$ (3)

where $q^{m}$ denotes wave vectors in the moire in 1BZ. In Fig.3(a) black dots
at the X points of 1BZ are due to the ZZ order in a single lattice.

For twisted bilayers with $\theta_{1,2}$, $\theta_{2,3}$ and for the case of
strong interlayer interactions $h=0.1$ the black dots at the M points of 1BZ
signal the AF toroidic order in Figs.3(b-c).

In order to show that the magnetic phases of Fig.1 correspond to the ground
states of the bilayer systems in the range of $\theta$ and $h$ studied here,
in Fig.3(d) the dipolar energy of ZZ phases is compared with that of toroidal
and helical phases in the limit of strong interlayer interactions $h=0.1$. For
$\theta>0.2$ it is apparent that the ZZ phase becomes too expensive, driving
dipoles of both layers to rotate to give way to toroidal and helical phases.

As a hallmark of a ferroic order as the toroidal order, there must be a
conjugate field to manipulate the order parameter. For ferrotoroidic
materials, this is a toroidal field, a magnetic vortex field violating both
space-inversion and time reversal symmetry, as we see next.

##  V Moire Torques

In the zig-zag equilibrium state at zero twists, the net torques on dipoles of
both layers cancel out when they point along any of the $\hat{n}_{i}$
directions. In the previous section, we demonstrated that a relative twist
between the layers affects such zig-zag pattern in each lattice and gives rise
to canted magnetic configurations whose patterns and motifs can be manipulated
by changing either $\theta$ or $h$. In this section, we show that the twist
between layers changes the distribution of internal magnetic fields in the
system and gives rise to an internal magnetic torque along the $z$ direction,
which is responsible for the twisted phases of Fig.1. Consider the origin of
the bilayer system $O^{\prime}=(0,0,0)$ at the center of both square lattices.
In units of $u_{0}/4\pi$ a dipole ${\bf m}_{i}^{\mu}$ in layer $\mu$ at
position ${\bf r}_{i}^{\mu}=(x_{i},y_{i},z_{i})^{\mu}$ respect to $O^{\prime}$
produces a magnetic field

${\bf B}(r_{i}^{\mu},r_{j}^{\nu})=\frac{3[{\bf m}_{i}^{\mu}\cdot({\bf
r}_{i}^{\mu}-{\bf r}_{j}^{\nu})]\cdot({\bf r}_{i}^{\mu}-{\bf
r}_{j}^{\nu})-{\bf m}_{i}^{\mu}|{\bf r^{\mu}}_{i}-{\bf
r}_{j}^{\nu}|^{2}}{|{\bf r}^{\mu}_{i}-{\bf r}_{j}^{\nu}|^{5}}$ (4)

at point ${\bf r}_{j}^{\nu}$ in layer $\nu$. An infinitesimal rotation of the
coordinate axes through angle $-\theta$ about $\hat{z}$ transforms the
components of a magnetic field like those of other vectors such as the
coordinates. In addition, $\mathbf{B}(\mathbf{r})$ depends on coordinates
$\mathbf{r}$, which must now be expressed in terms of the new ones (see
Appendix B for details). The overall infinitesimal change in $\mathbf{B}$,
combining the two changes becomes:

$\delta\mathbf{B}^{\mathrm{twist}}=\bm{\theta}\times\mathbf{B}-\bm{\theta}\cdot(\mathbf{r}\times\nabla)\mathbf{B},$
(5)

At zero twist, the bilayer system settles in the zig-zag phase. If we take
$\pm\hat{y}$ as the direction of the ferromagnetic collinear dipolar chains,
we have that the interlayer magnetic field created by a dipole i in layer
$\mu$, at points j of lattice $\nu$ is directed along the $\pm\hat{\bm{y}}$
direction (Appendix B). As a consequence, interlayer magnetic field and
dipoles are parallel and the total magnetic torque due to interlayer dipolar
interactions cancels out $\mathcal{T}(\theta=0)=\sum_{i,j}{\bf
m}_{j}^{\nu}\times{\bf B}_{i,(\theta=0)}^{\mu}=0$.

Consider next a finite twist $\theta$. Following Eq.5, the change in ${\bf
B}_{i}$ due to $\theta$ is given by:

$\delta\mathbf{B}_{i}^{\mathrm{twist}}=\theta
B_{i,(\theta=0)}(\hat{\bm{z}}\times\hat{\bm{y}})-\theta(\hat{\bm{z}}\cdot(\mathbf{r}\times\nabla))\mathbf{B}_{i},$
(6) $\delta\mathbf{B}_{i}^{\mathrm{twist}}=(\theta
B_{i,(\theta=0)},(-x\frac{\partial B_{i,(\theta=0)}}{\partial
y}+y\frac{\partial B_{i,(\theta=0)}}{\partial x}),0),$ (7)

Computing the total interlayer torque due to the twist, once again the
$\hat{\bm{y}}$ component of $\delta\mathbf{B}_{i}^{\mathrm{twist}}$ does not
produce torque. However the twist originates a magnetic field component along
$\hat{\bm{x}}$ which was absent at zero twist. The component orthogonal to the
dipolar chains gives rise to a torque along $\hat{\bm{z}}$,

$\delta\mathcal{T}^{\mathrm{twist}}(h,\theta)=\theta
B_{(\theta=0)}(\hat{\bm{z}})$ (8)

in units of g. Being directed along $\hat{\bm{z}}$,
$\delta\mathcal{T}^{\mathrm{twist}}$ rotates dipoles in the $\hat{x}-\hat{y}$
plane and gives rise to the toroidal phases of Fig.1.

##  VI Conclusion

Dipolar interactions in magnetic insulators with spins larger than one are
more frequent than thought a decade ago. Even when their long-range feature
contributes to stabilizing magnetic order, it is its anisotropy, the coupling
of lattice and magnetic degrees of freedom, which makes it of particular
interest for its role in non-collinear spin systems and topological magnetism.

Here, inspired by new research on twisted magnets, we have studied its effect
on the magnetic ground states of twisted square bilayers of dipoles. At zero
twist, the bilayers settle in a zig-zag antiferromagnetic order, and the
interlayer magnetic field has a single component along the zig-zag chains,
which gives zero net torque. We found that combining interlayer dipolar
interactions between dipoles with an easy plane rotation and an out-of-plane
relative rotation between dipolar layers gives rise to a net magnetic torque
orthogonal to the bilayer system. The twist torque is due to the emergence of
a new interlayer in-plane magnetic field produced by the relative twist
between layers. In particular, the component of the twist field directed
orthogonal to the zig-zag chains can rotate dipoles from their equilibrium
positions, and for certain moire angles, it gives rise to new magnetic orders
determined by the underlying moire pattern. The new magnetic phases are
magnetization-free and break time and inversion symmetry while depicting
chiral motifs that realize a nonzero toroidal moment. Thus, the emergent moire
field orthogonal to the zig-zag chains is conjugated to the toroidal moment
and, depending on the twist and the interlayer distance (that sets the
interlayer coupling), can tune them to order in antiferromagnetic or zig-zag
periodic fashions. We expect that our findings will motivate new research on
moire versions of other types of dipolar lattices. In particular we anticipate
that the study of spin waves in twisted bilayers of dipoles should manifest
nontrivial topological effects.

###### Acknowledgements.

P.M. X.C and I.T acknowledge support from Fondecyt under Grant No. 1210083.
This research was supported in part by the National Science Foundation under
Grant No. NSF PHY-1748958.

##  Appendix A Molecular dynamics simulations

We modeled the set of magnetic dipoles as a system of $n$ mechanical rods. The
position ${\bf r}_{i}$ ($i=1,\ldots,n$) of the center of mass of dipole $i$
remains fixed during the dynamics, and the orientation $\hat{\bf
m}_{i}=(\cos{\alpha_{i}},\sin{\alpha_{i}},0)$ of dipole $i$ is constrained to
the $\hat{x}-\hat{y}$ plane due to anisotropy. The magnetic state of the
system is specified by the dynamical angular variables $\alpha_{i}$. Magnetic
dipoles interact through the full long-range dipolar potential with all other
dipoles in the system. The equation of motion of each dynamical angle
$\alpha_{i}$ of each inertial dipole reads

$I\frac{d^{2}\alpha_{i}}{dt^{2}}=\mathcal{T}_{i}^{z}-\xi\frac{d\alpha_{i}}{dt}\,\text{,}$
(9)

where $I$ denotes the moment of inertia of the magnets and $\xi$ is the
damping for the rotation of dipoles respect to their local axis. The first
term at the right hand side of equation 9 is the $z$-component of the magnetic
torque due to the action of the internal magnetic field from dipolar
interaction between all dipoles in the system $\mathcal{T}_{i}^{z}=(\hat{\bf
m}_{i}\times\rm{\bf B}_{i})\cdot\hat{\bf z}$, where $\rm{\bf B}_{i}={\partial
U_{d}}/{\partial\hat{\bf m}_{i}}$ denotes the internal magnetic field produced
by all dipoles but the $i$-$th$ at the position of $\hat{\bf m}_{i}$ and the
potential energy defining the dipolar magnetic interactions is given by Eq.1.
Following previous experiments Mellado et al. (2023) we used in the
simulations dipolar rods with radius $r=0.79\times 10^{-3}$ m, length
$L=12.7\times 10^{-3}$ m, and mass $=0.19\times 10^{-4}$ kg, saturation
magnetization $M_{s}=1.05\times 10^{6}$ $\rm A/m$, inertial moment
$I=1.53\times 10^{-9}$ $\rm kg\,m^{2}$ and magnetic moment $m_{0}=0.026\times
10^{6}$ $\rm A\,m^{2}$. The second term on the RHS of Eq.9 corresponds to a
damping term. The viscosity $\xi(\alpha_{i})$ consists of a static and dynamic
contributions. Here we used experimental values reported in Mellado et al.
(2023). The static part the damping was computed directly by fitting the
relaxation of a single rod under a perturbation to the solution
$\alpha\sim\exp{-t/\tau_{D}}$. The estimated damping time $\tau_{D}\sim I/\xi$
of a single rod is $\tau_{D}=0.83$ s. The dynamics component is due to the
rotation of dipoles in the lattice with other magnets and therefore is
dependent on the orientation of the dipole, being maximum when it is
orientated perpendicular to $\hat{\bf x}$ axis.
We solved this set of coupled equations using a discretized scheme for the
integration of the differential equation, given by the recursion

$\theta_{i}(t+\Delta t)=f_{a}\theta_{i}(t)-f_{b}\theta_{i}(t-\Delta
t)+f_{c}\left(\hat{\bf z}\cdot\mathcal{T}_{i}(t)\right)\,\text{,}$ (10)

where the functions $f_{a}$, $f_{b}$ and $f_{c}$ corresponds to

$f_{a}=\frac{2}{1+\Delta t\,\eta(\theta_{i}(t))/2}\,\text{,}$ (11)
$f_{b}=\frac{1-\Delta t\,\eta(\theta_{i}(t))/2}{1+\Delta
t\,\eta(\theta_{i}(t))/2}\,\text{,}$ (12) $f_{c}=\frac{\Delta
t^{2}}{I(1+\Delta t\,\eta(\theta_{i}(t))/2)}\,\text{.}$ (13)

We initialized the simulations with two square lattices of $n_{s}\times n_{s}$
dipoles and lattice constant $a=2.2\times 10^{-2}\,\rm m$, separated by a
distance equal to $h$ (in units of a), and with an offset angle $\theta$
between the lattice vectors of each layer. The dynamical angles $\alpha_{i}$
where initialized in a random configuration. We solved the system of equations
Eq.9 numerically during a total time $t=3$ $\rm s$ using an integration step
of $\Delta t=4\times 10^{-6}$ $\rm s$. The dynamics are such that the damping
term drives the system to a magnetic configuration that locally minimizes the
dipolar interactions between magnetic dipoles.

We performed simulations with different interlayer distances $0.1<h<2$ and
different moire angles $0^{\circ}<\theta<45^{\circ}$ to characterize the
different stationary configurations that appear in the relaxation. We explored
system sizes up to $n=50$, and we used a replica method to avoid boundary
effects. The simulations were initialized with random magnetic configurations.

The dipolar interactions decay fast with the distance between dipoles $(i,j)$
as $r_{ij}^{-3}$. We implemented a cutoff in the computation of the
interactions, equal to $2a$, which shows results consistent with the dynamics
observed in simulations with full long-range interactions.

Fig.2 shows a coarse grain of the $z$-component of the toroidal moment T as a
function of real space. It can be written at each site as ${\bf T}_{i}={\bf
r}_{i}\times{\bf m}_{i}$ and it was calculated using the final magnetic
configurations obtained from the simulations. The value of the toroidal moment
at position $(x,y)$ was computed numerically by considering a cylinder with
center at position $(x,y)$ and a radius of $0.6a$, and then averaging the
toroidal moment of the dipoles that are located within the cylinder. We
considered cylinders spaced by $a/5$ to create the plots shown in Fig.2.

##  Appendix B Calculations details

###  B.1 Interlayer magnetic field in the ZZ phase

The magnetic field produced by a dipole i pointing along the $\hat{y}$ axis
and located at layer $\mu$ and sublattice $s_{1}$ in all points j of layer
$\nu$ can be separated into inter and intra-sublattice contributions:

${\bf
B}_{i}^{s_{1}}=\frac{(6mn,2m^{2}-4n^{2}-h^{2},3mh)}{(4n^{2}+m^{2}+h^{2})^{\frac{5}{2}}}$
${\bf
B}_{i}^{s_{2}}=\frac{(-3m(2n+1),-2m^{2}+(2n+1)^{2}+h^{2},-3mh)}{((2n+1)^{2}+m^{2}+h^{2})^{\frac{5}{2}}}$

in units of $g/m_{0}$ and where $n=(\Delta r_{ij}^{x})$, $m=(\Delta
r_{ij}^{y})$. After summation over all the j points of layer $\nu$, the x and
z components of $B_{i}$ cancel out to give only a non zero coponent along the
y direction,

$\displaystyle{\bf B}_{i}=\sum_{m}\big{[}\sum_{n\in
s_{1}}\frac{(0,2m^{2}-n^{2}-h^{2},0)}{(n^{2}+m^{2}+h^{2})^{\frac{5}{2}}}$ (14)
$\displaystyle+\sum_{n^{\prime}\in s_{2}}\frac{(0,-2m^{2}+n^{\prime
2}+h^{2},0)}{(n^{\prime 2}+m^{2}+h^{2})^{\frac{5}{2}}}\big{]}$

###  B.2 Effect of rotation of coordinates on a Magnetic field

An infinitesimal rotation of the coordinate axes through angle $-\theta$ about
$\hat{z}$ transforms the coordinates
$(x,y,z)\mapsto(x^{\prime},y^{\prime},z^{\prime})$ as follows Weyl (1929):

$x^{\prime}=x-y\theta,\quad y^{\prime}=y+x\theta,\quad z^{\prime}=z.$

For a rotation about a general axis $\hat{\mathbf{n}}$ through infinitesimal
angle $\theta$,

$\mathbf{r}^{\prime}=\mathbf{r}+\bm{\theta}\times\mathbf{r},$ (15)

where $\bm{\theta}=\theta\hat{\mathbf{n}}$. The inverse transform expresses
the old coordinates in terms of the new ones:

$\mathbf{r}=\mathbf{r}^{\prime}-\bm{\theta}\times\mathbf{r}^{\prime}.$

The magnetic field changes as a result of the rotation in two ways. First, it
is a vector quantity, so its components change like those of other vectors
such as the coordinate (15):

$\mathbf{B}^{\prime}=\mathbf{B}+\bm{\theta}\times\mathbf{B}.$ (16)

Second, $\mathbf{B}(\mathbf{r})$ depends on coordinates $\mathbf{r}$, which
must now be expressed in terms of the new ones $\mathbf{r}^{\prime}$:

$\displaystyle\mathbf{B}^{\prime}(\mathbf{r}^{\prime})$ $\displaystyle=$
$\displaystyle\mathbf{B}(\mathbf{r})=\mathbf{B}(\mathbf{r}^{\prime}-\bm{\theta}\times\mathbf{r}^{\prime})$
$\displaystyle=$
$\displaystyle\mathbf{B}(\mathbf{r}^{\prime})-[(\bm{\theta}\times\mathbf{r}^{\prime})\cdot\nabla^{\prime}]\mathbf{B}(\mathbf{r}^{\prime}).$

The overall infinitesimal change in $\mathbf{B}$, combining the two changes
(16) and (B.2), can be written in a compact form:

$\delta\mathbf{B}^{\mathrm{twist}}=\bm{\theta}\times\mathbf{B}-\bm{\theta}\cdot(\mathbf{r}\times\nabla)\mathbf{B}.$
(18)

where we have dropped the primes.

## References

  * Casher and Susskind (1974) A. Casher and L. Susskind, Physical Review D 9, 436 (1974).
  * Tomita et al. (2018) S. Tomita, H. Kurosawa, T. Ueda, and K. Sawada, Journal of Physics D: Applied Physics 51, 083001 (2018).
  * Thiaville and Fert (1992) A. Thiaville and A. Fert, Journal of magnetism and magnetic materials 113, 161 (1992).
  * Dzyaloshinskii (1960) I. E. Dzyaloshinskii, Soviet Physics JETP 10, 628 (1960).
  * Dzyaloshinskii (1958) I. Dzyaloshinskii, J. Phys. Chem. Solids 4, 241 (1958).
  * Heyderman and Stamps (2013) L. J. Heyderman and R. L. Stamps, Journal of Physics: Condensed Matter 25, 363201 (2013).
  * Togawa et al. (2016) Y. Togawa, Y. Kousaka, K. Inoue, and J.-i. Kishine, Journal of the Physical Society of Japan 85, 112001 (2016).
  * Mellado et al. (2023) P. Mellado, A. Concha, K. Hofhuis, and I. Tapia, Scientific reports 13, 1245 (2023).
  * Paula and Tapia (2023) M. Paula and I. Tapia, Journal of Physics: Condensed Matter 35, 164002 (2023).
  * Yu and Bauer (2021) T. Yu and G. E. Bauer, Chirality, Magnetism and Magnetoelectricity: Separate Phenomena and Joint Effects in Metamaterial Structures pp. 1–23 (2021).
  * Ray and Das (2021) S. Ray and T. Das, Physical Review B 104, 014410 (2021).
  * Shindou et al. (2013) R. Shindou, J.-i. Ohe, R. Matsumoto, S. Murakami, and E. Saitoh, Physical Review B 87, 174402 (2013).
  * Lucassen et al. (2019) J. Lucassen, M. J. Meijer, O. Kurnosikov, H. J. Swagten, B. Koopmans, R. Lavrijsen, F. Kloodt-Twesten, R. Frömter, and R. A. Duine, Physical review letters 123, 157201 (2019).
  * Nagaosa and Tokura (2013) N. Nagaosa and Y. Tokura, Nature nanotechnology 8, 899 (2013).
  * Mohylna et al. (2022) M. Mohylna, F. G. Albarracín, M. Žukovič, and H. Rosales, Physical Review B 106, 224406 (2022).
  * Camosi et al. (2018) L. Camosi, N. Rougemaille, O. Fruchart, J. Vogel, and S. Rohart, Physical Review B 97, 134404 (2018).
  * Schulz et al. (2012) T. Schulz, R. Ritz, A. Bauer, M. Halder, M. Wagner, C. Franz, C. Pfleiderer, K. Everschor, M. Garst, and A. Rosch, Nature Physics 8, 301 (2012).
  * Shaban et al. (2023) P. Shaban, I. Lobanov, V. Uzdin, and I. Iorsh, arXiv preprint arXiv:2304.08306 (2023).
  * Dong et al. (2015) S. Dong, J.-M. Liu, S.-W. Cheong, and Z. Ren, Advances in Physics 64, 519 (2015).
  * Spaldin et al. (2008) N. A. Spaldin, M. Fiebig, and M. Mostovoy, Journal of Physics: Condensed Matter 20, 434203 (2008).
  * Zimmermann et al. (2014) A. S. Zimmermann, D. Meier, and M. Fiebig, Nature communications 5, 4796 (2014).
  * Kuprov et al. (2022) I. Kuprov, D. Wilkowski, and N. Zheludev, Science Advances 8, eabq6751 (2022).
  * Astrov (1960) D. N. Astrov, J. Exptl. Theoret. Phys. (U.S.S.R.) 38, 984 (1960).
  * Date et al. (1961) M. Date, J. Kanamori, and M. Tachiki, Journal of the Physical Society of Japan 16, 2589 (1961).
  * Siratori et al. (1992) K. Siratori, K. Kohn, and E. Kita, Acta Phys. Pol. A 81, 431 (1992).
  * Ederer and Spaldin (2007) C. Ederer and N. A. Spaldin, Physical Review B 76, 214404 (2007).
  * Mellado et al. (2020) P. Mellado, A. Concha, and S. Rica, Physical Review Letters 125, 237602 (2020).
  * Lehmann and Lehmann (2022) J. Lehmann and J. Lehmann, Toroidal Order in Magnetic Metamaterials pp. 113–131 (2022).
  * Klebl et al. (2022) L. Klebl, Q. Xu, A. Fischer, L. Xian, M. Claassen, A. Rubio, and D. M. Kennes, Electronic Structure 4, 014004 (2022).
  * Xiao et al. (2020) Y. Xiao, J. Liu, and L. Fu, Matter 3, 1142 (2020).
  * Andrei et al. (2021) E. Y. Andrei, D. K. Efetov, P. Jarillo-Herrero, A. H. MacDonald, K. F. Mak, T. Senthil, E. Tutuc, A. Yazdani, and A. F. Young, Nature Reviews Materials 6, 201 (2021).
  * He et al. (2021) F. He, Y. Zhou, Z. Ye, S.-H. Cho, J. Jeong, X. Meng, and Y. Wang, ACS nano 15, 5944 (2021).
  * Tran et al. (2019) K. Tran, G. Moody, F. Wu, X. Lu, J. Choi, K. Kim, A. Rai, D. A. Sanchez, J. Quan, A. Singh, et al., Nature 567, 71 (2019).
  * Xiao et al. (2021) F. Xiao, K. Chen, and Q. Tong, Physical Review Research 3, 013027 (2021).
  * Li and Cheng (2020) Y.-H. Li and R. Cheng, Physical Review B 102, 094404 (2020).
  * Begum Popy et al. (2022) R. Begum Popy, J. Frank, and R. L. Stamps, Journal of Applied Physics 132 (2022).
  * Hejazi et al. (2020) K. Hejazi, Z.-X. Luo, and L. Balents, Proceedings of the National Academy of Sciences 117, 10721 (2020).
  * Yang et al. (2023) B. Yang, Y. Li, H. Xiang, H. Lin, and B. Huang, Nature Computational Science 3, 314 (2023).
  * Shimizu et al. (2021) K. Shimizu, S. Okumura, Y. Kato, and Y. Motome, Physical Review B 103, 184421 (2021).
  * Claassen et al. (2022) M. Claassen, L. Xian, D. M. Kennes, and A. Rubio, Nature communications 13, 4915 (2022).
  * Luo et al. (2022) Z.-X. Luo, U. F. Seifert, and L. Balents, Physical Review B 106, 144437 (2022).
  * Song et al. (2021) T. Song, Q.-C. Sun, E. Anderson, C. Wang, J. Qian, T. Taniguchi, K. Watanabe, M. A. McGuire, R. Stöhr, D. Xiao, et al., Science 374, 1140 (2021).
  * Kariyado and Vishwanath (2019) T. Kariyado and A. Vishwanath, Physical Review Research 1, 033076 (2019).
  * Hermann (2012) K. Hermann, Journal of Physics: Condensed Matter 24, 314210 (2012).
  * Weyl (1929) H. Weyl, Bulletin of the American Mathematical Society 35, 716 (1929).