# Data-driven scale identification in oscillatory dynamos

Anna Guseva
Department of Applied Mathematics, University of Leeds, Leeds LS2 9JT, UK
LERMA, Observatoire de Paris, Paris, France E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Parker’s mean-field model includes two processes generating large-scale
oscillatory dynamo waves: stretching of magnetic field lines by small-scale
helical flows, and by differential rotation. In this work, we investigate the
capacity of data-driven modal analysis, Dynamic Mode Decomposition, to
identify coherent magnetic field structures of this model. In its canonical
form, the only existing field scale corresponds to the dynamo instability. To
take into account multi-scale nature of the dynamo, the model was augmented
with coherent in time flow field, forcing small-scale magnetic field with a
faster temporal evolution. Two clusters of DMD modes were obtained: the “slow"
cluster, located near the dynamo wave frequency and associated with its
nonlinear self-interaction, and the “fast" cluster, centered around the
forcing frequency and resulting from the interaction between the wave and the
flow. Compared to other widely used methods of data analysis, such as Fourier
transform, DMD provides a natural spatiotemporal basis for the dynamo, related
to its nonlinear dynamics. We assess how the parameters of the DMD model, rank
and delay, influence its accuracy, and finally discuss the limitations of this
approach when applied to randomly forced, more complex dynamo flows.

###### keywords:

magnetic fields – dynamo – (magnetohydrodynamics) MHD

††pubyear: 2023††pagerange: Data-driven scale identification in oscillatory
dynamos–LABEL:lastpage

##  1 Introduction

Large-scale magnetic fields play a key role in formation and evolution of
planets, stars and accretion discs. In the Sun and some low-mass stars, these
fields evolve periodically in time, resulting in stellar magnetic cycles
(Muñoz-Jaramillo & Vaquero, 2019; Saikia et al., 2016). Mean-field
electrodynamics explains sustained large-scale astrophysical magnetic fields
through systematic stretching and twisting of magnetic field lines by
turbulent flows (Parker, 1955; Moffatt, 1978; Krause & Raedler, 1980). This
process is efficient when fluctuating flow $u^{\prime}$ and fluctuating field
$b^{\prime}$ correlate well, so that the mean electromotive force, or emf
$\varepsilon=\langle u^{\prime}\times b^{\prime}\rangle$, is non-zero. In
mean-field theory, the emf is assumed as a truncated expansion in the large-
scale magnetic field (Charbonneau, 2020). In isotropic turbulent flows non-
zero kinetic helicity, the first component of this expansion is
$\varepsilon_{i}=\alpha\delta_{ij}B_{j}$. In a nutshell, the dynamo can be
described with the one-dimensional $\alpha-\Omega$ model,

$\displaystyle A_{t}$ $\displaystyle=$ $\displaystyle\alpha
B+\eta(A_{xx}-l^{2}A),$ (1) $\displaystyle B_{t}$ $\displaystyle=$
$\displaystyle\Omega^{\prime}A_{x}-B^{3}+\eta(B_{xx}-l^{2}B),$ (2)

where $\eta$ is magnetic diffusivity (Proctor, 2007; Richardson & Proctor,
2010). The $x$-coordinate represents latitude or vertical direction in the
flow; $l$ is an inverse length scale of the field in the other spatial
directions. $\Omega$-effect corresponds to stretching of poloidal magnetic
field with potential $A$ into toroidal field $B$ by large-scale differential
rotation $\Omega^{\prime}$; $\alpha$-effect parametrizes generation of $A$
from $B$ by helical turbulent motions. Depending on the form of $\alpha$,
equations (1-2) admit steady or oscillatory linearly unstable solutions of
dipolar or quadrupolar parity (Richardson & Proctor, 2010). The oscillatory
solutions, travelling latitudinally and resembling solar magnetic activity,
are also known as dynamo waves. Mean-field model (1-2) is kinematic, i.e.
magnetic field grows linearly from a prescribed flow parametrized by
$\Omega^{\prime}$ and $\alpha$. It can be extended to include shear in the
solar tachocline (Charbonneau & MacGregor, 1997), flux transport by meridional
flows (Babcock, 1961; Leighton, 1969; Charbonneau, 2014; Cameron et al.,
2018), and dynamical feedback of the field on the flow through Lorentz force
(Tobias, 1997; Bushby, 2003). Large-scale dynamo cycles observed in direct
numerical simulations (DNS) of convective turbulence are often consistent with
the mean-field theory (Schrinner et al., 2011; Racine et al., 2011; Käpylä et
al., 2013).

Symmetry-preserving nonlinear term $B^{3}$ was introduced into (2) to model
saturation of the dynamo instability to a steady state. This can happen due to
expulsion of magnetic flux tubes from active dynamo regions by magnetic
buoyancy (Parker, 1979), the feedback from the Lorentz force slowing down
differential rotation (Gilman, 1983) or reducing stretching properties of the
flow ($\alpha$-effect) when magnetic energy becomes comparable to kinetic
energy of the flow (e.g., Jones et al., 2010). In the astrophysically relevant
limit of weak magnetic diffusion, generation of magnetic field at small scales
is favoured over large ones (e.g., Brandenburg & Subramanian, 2005), implying
very low saturation levels for large-scale magnetic fields (Vainshtein &
Cattaneo, 1992; Cattaneo & Hughes, 1996). However, numerical simulations of
Tobias & Cattaneo (2013) showed that strong shear, developing in many
astrophysical flows, promotes large-scale dynamo waves by suppressing
fluctuations at small scales. Shear may also facilitate temporally coherent
magnetic patterns in self-sustained magnetorotational turbulence without
helical or convective forcing (Nauman & Pessah, 2016). Thus, dynamos organize
at large scales through interactions between turbulence, shear, small-scale
and large-scale fields (see the reviews by Brandenburg (2018), Rincon (2019),
Tobias (2021) and references therein). To understand these nonlinear
interactions, it is necessary to unambiguously identify coherent dynamo
structures in multiscale dynamo models like state-of-art DNS (with about
$10^{6}$ degrees of freedom).

Spatiotemporal coherence of dynamo waves implies that separated in space fluid
parcels evolve synchronously while waves propagate through the flow, and there
is a particular spatial structure of magnetic and velocity fields associated
with this wave. Some insight into this structure can be gained from Fourier
transform or wavelet analysis; however, results are inevitably biased towards
a priori choice of spatial and temporal basis inherent in these methods.
Another way is to analyse eigenmodes of the dynamo, obtained by linearizing
its equations about some $n$-dimensional basic state. However, this
linearization is valid only in the vicinity of parameters where it was
performed and does not necessarily capture nonlinear behaviour of the system
away from this region. A more general approach is to obtain these structures
directly from dynamo data using the methods of data-driven modal
decomposition. Frequently applied for analysis of hydrodynamic turbulence
(Taira et al., 2020), they have not been widely employed in astrophysics; yet
they were found useful for torsional magnetic wave detection (Hori et al.,
2023). In essence, modal decomposition is a factorization of the matrix of
flow observables into matrices representing their spatial structure, temporal
dynamics, and amplitudes indicating relative importance of every spatial mode.
Proper Orthogonal Decomposition (POD) produces modes that give an optimal
representation of the flow in terms of energy compared to any other basis of
the same dimension (Lumley, 2007; Sirovich, 1987; Holmes et al., 2012); it can
mix different time and length scales in a single mode if their energy
contribution is comparable. Dynamic Mode Decomposition (DMD) seeks for the
closest approximation of the flow data in terms of a linear system, and is
efficient for analysing periodic and quasi-periodic dynamics (Schmid, 2010;
Bagheri, 2010; Schmid, 2022). POD and DMD are entirely data-driven and
therefore can be applied independently of boundary conditions and flow
complexity. Moreover, compared to empirical temporal analysis such as Hilbert-
Huang transform (Huang et al., 1998), these methods also provide physically
interpretable spatial structures. Unlike Fourier and wavelet transforms, no a
priori spatiotemporal basis is assumed; both spatial shape and temporal
behaviour of the modes result from the data analysis. Assumptions about
linearity and stationarity of the signals are relaxed for Dynamic Mode
Decomposition, as compared to Fourier transform. As a drawback, DMD requires
two empirical, user-defined parameters, rank and delay, and we will discuss
below the influence of these parameters on decomposition accuracy. POD and DMD
result in large modal bases for very chaotic signals; this is not the case for
dynamo mean-field models analysed here.

In this paper, we will test the applicability of these methods for detection
of scales and their interaction on two dynamo problems: Parker’s
$\alpha-\Omega$ model (1-2) and an extended dynamo model featuring both large
and small scales. This paper is structured as follows: in section 2 we
describe the data-driven approach, POD and DMD methods; in section 3 we apply
this methods to analysis of Parker’s dynamo waves; in section 4 we construct a
more complex dynamo model with scale separation, analyse its dynamics with
DMD. In section 5 we assess the accuracy of the DMD model. Section 6 concludes
the paper.

##  2 Data-driven approach

Figure 1: Time evolution of dynamo waves in system (1-2): magnetic potential
$A$ (top), toroidal component of magnetic field $B$ (bottom). The data are
collected until $t=100$; $D=30$.

Consider the original nonlinear system (1-2) or an analogous one which evolves
the magnetic potential $A$ and the magnetic field $B$ and depends on a set
$\bm{\mu}$ of dimensionless functions and parameters,

$\frac{d\bm{q}}{dt}=\bm{f}^{nonlinear}(\bm{q},t,\bm{\mu}),\quad\bm{q}(t)=(A,B),\quad\bm{\mu}=\alpha,\Omega^{\prime},\eta,l...$
(3)

The functional form of the $\alpha$-effect and the boundary conditions are

$\alpha=-\alpha_{0}\sin(2\pi x/L),\quad A=B=0\text{ at }x=0,L=10.$ (4)

Dynamo number $D=\Omega^{\prime}\alpha_{0}/\eta^{2}l^{3}$ is the dimensionless
parameter of the system defining frequency of waves; if
$\eta=l=\Omega^{\prime}=1$, then $D=\alpha_{0}$. We solve system (3)
numerically in Python spectral solver Dedalus (Burns et al., 2020), with
Chebyshev discretization of $N=256$ points in $x$ and Runge-Kutta time-
stepping method. The dynamo appears as a linear instability in the form of
traveling waves for $D>20$, as shown in figure 1. The initial conditions were
of dipolar parity, with $A$ symmetric and $B$ antisymmetric in $x$. This way,
the system is constrained to dipolar solutions; otherwise quadrupolar dynamos
with antisymmetric $A$ and symmetric $B$ could also appear.

During the simulation, we collect instantaneous spatial profiles (snapshots)
of the numerical solution of the system, $\bm{q}_{k}(x)\in R^{N}$, at time
steps $k=1,2,\dots,K$, with frequency sufficient to resolve the fastest flow
dynamics, and assemble them into the data matrix $Q$,

$Q=[\bm{q}_{1}\quad\bm{q}_{2}\cdots\bm{q}_{K}],$ (5)

for further analysis. Here, we performed POD and DMD analysis separately on
the data for $A$ and $B$, although treating $A$ and $B$ as a single system
state vector is also possible.

###  2.1 Proper Orthogonal Decomposition

Proper Orthogonal Decomposition (POD), also known as Principal Component
Analysis, or Karhunen–Loeve expansion, is based on the idea that spatial flow
correlations can be decomposed into orthogonal vectors ranked by their energy.
Given system state $\mathbf{q}$, defined on the domain $\mathbf{x}$, POD modes
are eigenfunctions of the integral equation

$\int\langle\mathbf{q(x)}\mathbf{q^{*}(x^{\prime})}\rangle\phi(\mathbf{x^{\prime}})\mathbf{dx^{\prime}}=\sigma^{2}\phi(\mathbf{x}),$
(6)

with autocorrelation function
$\langle\mathbf{q(x)}\mathbf{q^{*}(x^{\prime})}\rangle$ as a kernel, and
$\langle\dots\rangle$ denoting temporal average. The eigenvalues
$\sigma^{2}\geq 0$ correspond to the average energy content in each mode
$\phi$, and the first $r$ modes represent more of the system energy then any
other spatial basis of order $r$. A rigorous introduction to POD analysis for
fluid flows can be found in Holmes et al. (2012, pp. 86-128). For discretized
data (5), modes $\phi$ are eigenvectors of positive semi-definite correlation
matrix $Q^{*}Q$. In practice, it is convenient to find $\phi$ by factorizing
the data matrix $Q$ with Singular Value Decomposition (SVD),

$Q=\Phi\Sigma V^{*}\approx\Phi_{r}\Sigma_{r}V_{r}^{*},$ (7)

which generalizes the eigenvalue decomposition to non-square matrices. The
data set is thus represented as a product of the matrix $\Phi$ of the POD
modes $\phi_{i}$, the matrix $\Sigma$ of their singular values $\sigma_{i}$,
and their temporal evolution coefficients $V$. $Q$ is typically low-rank in
fluid flows, so most of the singular values are nearly zero; the corresponding
modes contain very little energy. A common practice is to retain a subset $r$
of modes with energy content larger than a threshold, and to discard the rest
as noise. To simplify our analysis, we interpolated all our data from non-
equispaced Chebyshev points to an equispaced grid; otherwise, an energy weight
matrix is required in (6,7) (see Schmidt & Colonius, 2020).

###  2.2 Dynamic Mode Decomposition

Unlike statistics-based POD, Dynamic Mode Decomposition (DMD) harnesses
temporal dynamical behaviour of the flow. It is closely related to Koopman
operator theory, proposing that the temporal dynamics of nonlinear system (3)
can be described with a linear operator acting on suitable flow observables
(Koopman, 1931). For an arbitrary nonlinear system, there is no generalized
method of constructing an exhaustive set of such observables and the operator
itself. DMD provides an approximation to eigenvectors and eigenvalues of
Koopman operator from the data sequence, without calculating the operator
explicitly (Tu et al., 2014). The reader is referred to the comprehensive
review by Schmid (2022) for a rigorous derivation of the DMD algorithm and its
links to the Koopman theory.

####  2.2.1 Standard DMD

In this work, we use the “exact" definition of the DMD algorithm proposed by
Tu et al. (2014). We seek to approximate nonlinear system (3), discretized in
time, with a linear operator $\mathcal{A}$ that maps the flow state at time
$t_{k}$ to the flow state at time $t_{k+1}$,

$\bm{q}_{k+1}=\mathcal{A}\bm{q}_{k},\quad\text{where}\quad\bm{q}_{k}=\sum_{i=1}^{n}\psi_{i}\lambda_{i}^{k}b_{i}(0),$
(8)

where $b_{i}(0)$ is the initial magnitude of the eigenvector $\psi_{i}$.
Discrete eigenvalues $\lambda_{i}$ define whether $\psi_{i}$ are decaying,
neutral, or growing in time. We augment the data matrix (5) with one more
snapshot of the system, and construct another matrix $Q^{\prime}$, such that

$Q^{\prime}=[\bm{q}_{2}\quad\bm{q}_{2}\cdots\bm{q}_{K+1}],\qquad
Q^{\prime}\approx\mathcal{A}Q.$ (9)

A least-squares approximation to $\mathcal{A}$ can be obtained by calculating
the matrix product of the pseudoinverse of $Q$ with $Q^{\prime}$. However, for
low-rank matrices like $Q$, pseudoinverse is ill-conditioned and leads to
spurious results (Press et al., 1992). A workaround proposed by Schmid (2010)
is to construct a lower, $r$-dimensional representation of matrix $Q$ by
truncating its singular value decomposition (7) at rank $r$. After this
operation, (9) is re-arranged as

$Q^{\prime}\approx\mathcal{A}\Phi_{r}\Sigma_{r}V_{r}^{*},\qquad\Phi_{r}^{*}\mathcal{A}\Phi_{r}=\Phi_{r}^{*}Q^{\prime}V_{r}\Sigma^{-1}_{r}=\mathcal{A}_{r},$
(10)

where $\mathcal{A}_{r}$ is a projection of $\mathcal{A}$ on the POD modes
$\Phi_{r}$. Its eigenvalues and eigenvectors can be calculated from
$\mathcal{A}_{r}\tilde{\psi}=\lambda\tilde{\psi}$. The full-state eigenvectors
in the physical space can be recovered from $r$-dimensional $\tilde{\psi}$
through

$\psi=\frac{1}{\lambda}Q^{\prime}V_{r}\Sigma^{-1}_{r}\tilde{\psi},\qquad\omega=\ln(\lambda)/\Delta
t.$ (11)

The second equation in (11) relates discrete-time eigenvalues $\lambda$ to
continuous-time eigenvalues $\omega$, both complex. In the following, we will
refer to $\psi$ and $\omega$ as DMD modes and DMD eigenvalues. If the real
part of the eigenvalue $\Re(\omega)<0$, the mode is dampened; if
$\Re(\omega)>0$, the mode will be growing. In a steady state system, its
significant modes are expected to be nearly neutral, $\Re(\omega)\approx 0$.
The imaginary part $\Im(\omega)$ is the temporal frequency of the mode.

As in section 2.1, threshold $r$ is user-defined. A common choice for systems
with little noise is to retain first $r$ modes so that
$\sum_{r}\sigma_{r}/\sum_{tot}\sigma=0.99$, retaining 99% information from the
original data set. An alternative cut-off criterion is
$\sum^{2}_{r}\sigma^{2}_{r}/\sum^{2}_{tot}\sigma^{2}=0.99$, bearing in mind
that $\sigma^{2}$, and not $\sigma$, represents the energy of POD modes in
(6). The latter criterion means that 99% of energy is retained and results in
a sparser spatial basis. The relative error of the DMD model can be estimated
as

$\epsilon=||Q-Q^{model}||_{2}/||Q||_{2},$ (12)

where $||\cdot||$ is $L_{2}$-norm and $Q^{model}$ is the approximation to $Q$
calculated according to (8).

###  2.3 High-order (Hankel) DMD

The robustness of DMD in ergodic dynamical systems can be improved by
embedding state vectors $q_{k}$ in a higher dimension. Using the method of
delays (Takens, 1981), we can construct Hankel matrix $Q_{H}$ from the
original data (5) as

$Q_{H}=\left[\begin{array}[]{cccc}\bm{q}_{1}&\bm{q}_{2}&\cdots&\bm{q}_{K-d+1}\\\
\bm{q}_{2}&\bm{q}_{3}&\cdots&\bm{q}_{K-d+2}\\\ \vdots&\vdots&&\vdots\\\
\bm{q}_{d}&\bm{q}_{d+1}&\cdots&\bm{q}_{K}\\\ \end{array}\right],$ (13)

where $d$ is the delay parameter, and augment (9) in the same way. This is
done because our initial choice of system observables, $\bm{q}_{k}$, is not
necessarily contained within the finite-dimensional invariant subspace of the
Koopman operator. In Hankel matrix (13), we increased the number of
observables: the matrix rank of $Q_{H}$ is higher then the rank of $Q$. Arbabi
& Mezic (2017) proved that in the limit $d\to\infty$ DMD eigenvalues and
eigenmodes approximating the system

$Q^{\prime}_{H}\approx\mathcal{A}Q_{H}.$ (14)

converge to Koopman modes and eigenvalues from the finite-dimensional
invariant subspace of $\mathcal{A}$. Hankel DMD is also viewed as a higher-
order approximation of (9) (Le Clainche & Vega, 2017).

Construction of (13) introduces another user-defined, delay parameter $d$, and
its choice is not straightforward. Takens (1981) stated that an attractor of
dimension $m$ can be embedded into a space with $d=2m+1$ dimensions. Broomhead
& King (1986), considering physically relevant time scales of the system,
suggested to choose $d$ so that the SVD spectrum of (13) converges as the
window length of delay $t_{k+d}-t_{k}$ is varied. In practice, $d$ is
increased until the desired accuracy of approximation is achieved in the sense
of $L_{2}$-norm (Fujii et al., 2019).

##  3 Parker’s dynamo waves

As a benchmark problem, we take the system (1-2) augmented with (4), as
described in the beginning of section 2. We analyse the steady-state data sets
of toroidal magnetic field $B$ and poloidal magnetic potential $A$, depicted
in figure 1.

(a)

(b)

(c)

(d)

Figure 2: (a) Singular value spectrum of the dynamo waves in figure 1, system
(1-2). Empty symbols, $\sigma_{i}/\sum_{i}\sigma$, filled symbols,
$\sigma^{2}_{i}/\sum_{i}\sigma^{2}$. Dashed line represent 99% cut-off for
both $\sigma$\- and $\sigma^{2}$-criteria. (b) The first four POD modes of
dynamo waves. Solid line, the first of the POD mode pairs in panel (a);
dashed: the second mode in the pair. (c) Reduced model of the dynamo waves in
figure 1 based on: top, $\phi^{B}_{0}$; middle, $\phi^{B}_{1}$, bottom:
$\phi^{B}_{0}$ and $\phi^{B}_{1}$ together. (d) DMD modes of the dynamo waves,
calculated using the first four POD modes from panel (c). Solid lines: the
first mode of the complex-conjugate modal pair. The real part of the mode is
represented by a darker color and the imaginary part by a lighter (gray)
color. Dashed, the same but for the second mode $\psi_{1}$ in the pair
$\psi^{*}_{1}=\psi_{0}$ (only for $\psi^{A}_{0,1}$).

###  3.1 POD decomposition

Figure 2(a) shows the singular value spectrum of $A$ and $B$. The spectrum
decays rapidly, indicating the low rank of the data and absence of noise. $B$
has more complex structure, with $4$ modes needed to describe 99% of the data
set according to $\sigma$-criterion, compared to $2$ modes for $A$. This
difference can be attributed to the nonlinear term $B^{3}$ in equation (2),
and therefore more complex dynamics. Potential $A$ is influenced by this
nonlinearity only indirectly through $\alpha$-effect. According to
$\sigma^{2}$-criterion, 99% of energy in both $A$ and $B$ is captured by the
first two modes. The singular values of $A$ and $B$ come in nearly equal
pairs, a common feature in fluid flows with symmetries (Deane et al., 1991).
The dipolar symmetry implies that

$x\to-x,\quad A\to A,\quad B\to-B.$ (15)

If the solution of (3) is invariant under the action of symmetry group $G$,
and $\phi(x)$ is an eigenfunction of (6), then for every element $g\in G$ such
as (15), $g\circ\phi$ is also an eigenfunction of (6) with the same eigenvalue
(Aubry et al., 1993, Proposition 3).

The POD modes, corresponding to the singular value pairs, are shown in Figure
2(b). Following the fields parity in figure 1, the $A$-modes $\phi^{A}_{0-3}$
are symmetric, and $B$-modes $\phi^{B}_{0-3}$ are antisymmetric with respect
to $x$. The most energetic modes $\phi^{A}_{0,1}$, $\phi^{B}_{0,1}$ are large-
scale and correspond to the length scale of the dynamo wave observed in figure
1. The spatial scale of the second, less energetic mode pair is approximately
two times smaller, and its temporal dynamics is faster. Each mode in the pair
$\phi^{B}_{0,1}$ (or $\phi^{A}_{0,1}$), weighted by its temporal coefficient
(see equation 24 in section 5), corresponds to a standing wave shifted in
phase with respect to the other (figure 2c, upper and middle panel). They
enhance each other at some spatial locations, and weaken at others, producing
travelling wave dynamics observed in figure 2 (c), bottom panel. The relative
error between this 2-mode model and toroidal field $B$ in figure 1, evaluated
as $L_{2}$ norm, is only 3.7%.

###  3.2 DMD decomposition

We calculate the DMD eigenmodes and eigenvalues of the dynamo waves in figure
1 using the algorithm in section 2.2. The previous section shows that
truncation of rank $r=4$ is sufficient to reproduce the data reliably. As they
are real, the matrix $\mathcal{A}_{r}$ in (10) is also real and its
eigenvalues come in complex conjugate pairs with nearly zero real parts, as
shown in table 1.

The first frequency pair, $\Im(\omega_{0,1})=\pm 3.558$, corresponds to the
frequency of dynamo linear instability at $D=30$, and is the main frequency
observed in figure 1. The second frequency is three times larger than the
first one, $\Im(\omega_{2,3})\approx 3\Im(\omega_{0,1})$, and can be
interpreted as an interaction of the first mode with itself through nonlinear
damping $B^{3}$ in (2), as the periodic temporal behaviour implies
$\sin^{3}(\omega)\sim\sin(3\omega)$. Note the implication of choosing $\sigma$
or $\sigma^{2}$ cut-off: if only 2 POD modes are retained according to
$\sigma^{2}$-criterion, the dynamics will be projected on the most energetic,
slowly oscillating modes, and the fast-oscillating modes will not be
recovered. The choice of cut-off $r$ is therefore a balance between desired
accuracy and basis sparsity.

The first two panels in figure 2(d) show the “slow" DMD modes of $A$ and $B$;
one can observe that the $x$-profile of the real parts of modes is similar to
one of the corresponding POD mode pairs in figure 2(b), while their imaginary
part is resembles the other. Like the eigenvalues, the DMD modes are also
complex conjugates, $\psi^{*}_{1}=\psi_{0}$, thus, their real parts are equal,
and their imaginary parts have equal magnitude along $x$ but opposite signs.
This allows to represent the travelling wave structures similar to the lower
panel in figure 2(c) with just one complex DMD mode $\psi_{0}$ and its
temporal coefficient. Indeed, consider a DMD mode
$\psi_{0}(x)=\psi^{r}+\psi^{i}\mathrm{i}$, and the corresponding projection of
the data set on this mode, $c_{0}(t)=c^{r}+c^{i}\mathrm{i}$. One-mode
approximation for the data matrix $Q$ would be

$Q\approx\psi^{r}c^{r}-\psi^{i}c^{i}+(\psi^{i}c^{r}+\psi^{r}c^{i})\mathrm{i}.$
(16)

When the corresponding contribution from the conjugate mode $\psi_{1}$ will be
added, the imaginary part of (16) will be cancelled, as the data matrix $Q$ is
real; the real part will be double of that of (16). Thus, it is enough to
track one mode from the complex-conjugate pair in $A$ and $B$, and DMD
provides a more compact vector basis for this system compared to POD. In the
following, we will focus on the DMD approach.

Mode | $\omega^{A}$ | $\omega^{B}$
---|---|---
$0$, $1$ | $-7.401\cdot 10^{-7}\pm 3.558\mathrm{i}$ | $-9.525\cdot 10^{-6}\pm 3.558\mathrm{i}$
$2$, $3$ | $-1.063\cdot 10^{-4}\pm 10.673\mathrm{i}$ | $-1.257\cdot 10^{-3}\pm 10.675\mathrm{i}$
Table 1: DMD eigenvalues of potential $A$ and toroidal magnetic field $B$,
obtained based on $r=4$ truncation of the data. The modes are ranked by their
growth rates (less decaying first). The eigenvalues were rounded up to third
digit.

##  4 Multi-scale dynamo model

Parker’s dynamo model (1-2) parametrises interaction of small-scale dynamo and
turbulence with $\alpha$-effect. Thus, it has only one independent time and
length scale, related to the dynamo linear instability at large scales. The
second, faster and smaller length scale, identified by POD and DMD, is the
result of magnetic field interacting with itself. It is also very weak, less
than 1% of the total magnetic energy, as can be seen in figure 2(a). However,
it is not the case for more realistic numerical models of dynamo where waves
arise naturally from shear turbulence. In those simulations, the flow is
driven by a combination of large-scale shear and a small-scale forcing,
coherent or randomized, and energy is distributed on a range of scales (Tobias
& Cattaneo, 2013; Pongkitiwanichakul et al., 2016; Nigro et al., 2017). More
advanced dynamo models allow to integrate large and small-scale dynamics
simultaneously, for example, with a shell model describing behaviour of
velocity and magnetic fields at small scales (Nigro & Veltri, 2011). In line
with such models, we augment system (1-2) with extra terms, mimicking the
multi-scale flow behaviour in the DNS, and study whether DMD can recognize
mixed small and large scales of the system. The augmented one-dimensional
kinematic dynamo equations are

$\displaystyle A_{t}$ $\displaystyle=$
$\displaystyle\varepsilon+\eta(A_{xx}-l^{2}A),$ (17) $\displaystyle B_{t}$
$\displaystyle=$
$\displaystyle\Omega^{\prime}A_{x}-B^{3}+\eta(B_{xx}-l^{2}B),$ (18)
$\displaystyle b_{t}$ $\displaystyle=$ $\displaystyle\mu u_{x}B-\gamma
b^{3}+\eta(b_{xx}-l^{2}b),$ (19) $\displaystyle\varepsilon$ $\displaystyle=$
$\displaystyle\varepsilon_{0}g(ub),\quad
g(f)=\frac{1}{\sqrt{2\pi}\sigma}\exp(-\frac{x^{2}}{2\sigma^{2}}),$ (20)
$\displaystyle u$ $\displaystyle=$ $\displaystyle u_{0}\sin(2\pi
k_{u}x/L\pm\omega_{u}t).$ (21)

Temporal and spatial evolution of the small-scale flow field $u$ is prescribed
by (21). Equation (19) describes the evolution of small-scale fluctuations $b$
of toroidal magnetic field and is inspired by the full induction equation for
magnetic fluctuations, see e.g. Tobias (2021, p. 29). Since the latter is
nonlinear in fluctuations due to the terms $\nabla\times(u^{\prime}\times
b^{\prime}-\overline{u^{\prime}\times b^{\prime}})$, nonlinearity affects all
scales of magnetic field, and dynamo should be able to saturate at both small
and large scales in statistically steady flows. Exact mechanism of this
saturation is still a subject of active research (Moffatt & Dormy, 2019). For
simplicity, we introduce the same form of nonlinearity $\gamma b^{3}$ into
(19), which allows to control the saturation amplitude of $b$. The dissipative
term in (19) is proportional to $b_{xx}$; $\mu u_{x}B$ models generation of
$b$ through interaction of the flow and the large-scale field. Parameters
$\gamma$ and $\mu$ are introduced to control the strength of nonlinearity and
the induction-like term, respectively, in this simplified mathematical model
of a multiscale dynamo. In a more realistic DNS simulation, signals from $B$
and $b$ would be mixed in multiscale magnetic field, $B_{s}=B+b$. The
interaction between $b$ and $u$ contributes to fluctuating electromotive
force, proportional to $ub$. Net emf $\varepsilon=\varepsilon_{0}g(ub)$ on the
scale of the dynamo waves in (17) is calculated by filtering product $ub$ with
a Gaussian filter (20) and ensuring $\varepsilon=0$ at the boundaries $x=0,L$.
Alternatively, time averaging or envelope smoothing can be employed to model
action of emf on large scales.

(a)

(b)

(c)

(d)

Figure 3: The data obtained by integrating in time multi-scale dynamo
equations (17-19). Without imposed $\alpha$-effect: multi-scale toroidal
magnetic field $B_{s}=B+b$ (a) and fluctuating toroidal magnetic field $b$
(b); flow and emf given by (20-21). System with controlled $\alpha$-effect:
$B_{s}$ (c), $b$ (d); flow and emf from (22-23). See text for description of
the system parameters.

Dynamo waves were consistently observed when $\varepsilon_{0}>800$ and the
rest of parameters set to $\Omega^{\prime}=\eta=l=1$, $\mu=5$, $u_{0}=0.5$. We
set $\gamma=0$, allowing the dynamo to saturate only through large scales.
Large value of $\sigma=30$ removes completely small-scale fluctuations from
$\varepsilon$ and thus inhibits non-oscillatory dynamo solutions (Richardson &
Proctor, 2010). Frequency $\omega_{u}$ and wave number $k_{u}$ of small-scale
$u$ were set to $\omega_{u}=13$, $k_{u}=13$ to ensure a faster temporal
dynamics with respect to the dynamo wave and at least a decade of scale
separation between the wave and the small-scale dynamo. Figures 3(a,c) show
results of numerical integration of (17-21) with these parameters and
$\varepsilon_{0}=1500$. The large-scale field $B$ is again dipolar, anti-
symmetric with respect to the equator; it is the dominant feature of the
dynamo waves observed in the total magnetic field $B_{s}$. However, it is
spatially and temporally modulated, allowing locally wider regions of positive
and negative field direction. These differences appear because spatial
distribution of net emf departs from the one in “classic" dynamo waves. On the
other hand, the fluctuating field $b$ in figure 3(c) is predominantly small-
scale, but yet with a systematic footprint from the large-scale $B$, symmetric
with respect to the $x/L=0.5$. Fast Fourier Transform (FFT) of multiscale
magnetic field $B_{s}$ over the temporal domain, with resulting spatial
Fourier modes integrated in $x$, shows that the temporal dynamics in this
system is complex and features numerous frequencies (figure 4a). The most
prominent peaks are clustered around the dominant dynamo ingredients: the
“slow" cluster corresponding to the dynamo wave and the “fast" cluster around
the flow frequency $\omega_{u}$. Other degrees of freedom are excited in
between these peaks, so the spectrum in the “slow" cluster is relatively
broad-band, which complicates identification of dynamo waves, small scales,
and their interaction. It is thus of interest to construct a multiscale dynamo
system featuring a well-defined $\alpha$-effect as in (4).

###  4.1 Imposing $\alpha$-effect

For this, we add one more component to the flow velocity,

$u=u_{0}(\sin(\pi x/L)+\sin(2\pi k_{u}x/L\pm\omega_{u}t)),$ (22)

so that magnetic fluctuations scale as $b\sim\cos(\pi x/L)B$. If electromotive
force is defined as

$\varepsilon=-\varepsilon_{0}g(ub)\sim\sin(2\pi x/L)B,$ (23)

then $\alpha$-effect is similar to (4), yet dynamo waves appear from
interaction of the flow and small-scale magnetic field. Formally, the extra
component in (22) can be interpreted as a large-scale one, although it is not
constant in $x$ compared to mean differential rotation $\Omega^{\prime}$ in
(18). Magnetic fluctuations will thus contain direct interaction of this
component in (22) and large-scale magnetic field in induction equation (19),
unlike in the mean-field theory. This inconsistency is resolved in the
previous multiscale dynamo system (17-21) without imposed $\alpha$-effect,
where large-scale shear interacts with $B$ indirectly through
$\Omega^{\prime}A_{x}$ and $u_{x}B$. Controlling $\alpha$-effect in this way
allows to obtain a multiscale dynamo with few degrees of freedom, where large
scales are similar to Parker’s dynamo waves (1-2).

We perform a numerical simulation of augmented model (17-19), together with
(22-23), with dipolar initial conditions. Model parameters were set to
$\varepsilon_{0}=30$, $\Omega^{\prime}=l=\mu=1$, $\eta=0.1$, $u_{0}=0.5$,
$\sigma=10$, $\omega_{u}=13$, $k_{u}=13$, resulting in “classic” large-scale
dynamo waves with a period $T\approx 20$, about ten times longer (figure 3b).
The small-scale nonlinearity was set non-zero, $\gamma=10$, in order to reduce
the magnitude of $b$ to values below those of large-scale $B$ (figure 3d), as
in dynamo without imposed $\alpha$-effect. The data were collected up to
$T=200$, covering about $10$ periods of wave evolution and thus resolving both
slow and fast dynamics. As previously, $b$ is mostly small-scale with a
symmetric footprint of the large-scale field, arising through forcing of $b$
with the product of anti-symmetric $u_{x}$ and anti-symmetric $B$. Note the
difference with multiscale dynamo system (17-21): the sign of the large-scale
footprint in $b$ is periodic in time, following the sign of
$u_{x}B\sim\cos(\pi x/L)B$. This is less apparent in figure 3(c) because the
product of $B$ and $u_{x}$ remains small-scale. The resulting multi-scale
magnetic field, $B_{s}=B+b$ is thus of mixed parity with respect to $x$;
however, since the magnitude of $b$ is about $2$ times smaller than the one of
$B$ and contains a significant energy contribution from small scales, this
asymmetry very weak. FFT transform of $B_{s}$ in time shows much simpler
temporal dynamics, with just a few dominant frequency peaks (figure 4b). In
the following, we will analyse this mixed signal $B_{s}$ with DMD, discarding
transients up to $t=50$.

###  4.2 DMD of the model with imposed $\alpha$-effect

(a)

(b)

(c)

(d)

Figure 4:  Mode identification of multiscale magnetic field $B_{s}$. (a)
Amplitudes and frequencies of the DMD modes, compared with FFT of the data
(dashed) of dynamo waves from (17-21), without imposing $\alpha$-effect
(figure 3a). Filled circles, standard DMD of rank $r=14$; filled squares,
Hankel DMD with $r=14$ and delay $d=20$; empty squares, Hankel DMD, $r=4$,
$d=20$. Dots corresponds to high-rank Hankel DMD, $r=40$, $d=300$. Only
positive DMD frequencies are shown. (b) FFT (dashed) and DMD spectrum of the
model with imposed $\alpha$-effect (figure 3b). Filled circles, standard DMD
of $r=14$; filled squares, Hankel DMD with rank $r=14$ and delay $d=9$; empty
squares, Hankel DMD, $r=6$, $d=9$; dots, high-rank Hankel DMD, $r=22$, $d=10$.
Crossed-out circles correspond to spurious eigenvalues that appear in high-
rank models ($r=22$) without increasing $d$. (c) Cumulative singular value
spectrum of the multi-scale dynamo model with (circles) and without (solid
lines) imposed $\alpha$-effect. Dark blue, $\sigma^{2}$; light blue, $\sigma$.
Vertical lines denote 99% cut-off in both spectra
($\sum_{i}\sigma_{i}(\sigma_{i}^{2})/\sum_{tot}\sigma(\sigma^{2})=0.99$ ,
dashed). (d) Profiles of DMD modes $\psi_{0}$, $\psi_{1}$, $\psi_{2}$ (slow
cluster) and $\psi_{4}$ (fast cluster). Solid, with imposed $\alpha$-effect;
dashed, without imposed $\alpha$-effect.

With imposed $\alpha$-effect Without imposed $\alpha$-effect Interpretation
DMD, HDMD HDMD FFT DMD HDMD HDMD $r=14$ $r=6$, $d=9$ $r=14$, $d=9$ $r=14$
$r=14$, $d=20$ $r=40$, $d=300$ Slow $\omega_{0}$ $0.282$ $0.282$ $0.282$
$0.209$ $0.257$ $0.244$ $0.229$ Dynamo wave, $\Im(\omega_{0})$ $\omega_{1}$
$0.845$ $0.847$ $0.670$ $0.687$ $\Im(\omega_{1})\approx 3\Im(\omega_{0})$
$\omega_{2}$ $1.435$ $1.426$ $1.131$ $1.267$ $1.145$ $\Im(\omega_{2})\approx
5\Im(\omega_{0})$ $1.591$ $1.499$ $0.687$ $\approx 7.6\Im(\omega_{0})$ $2.052$
$2.061$ $\approx 9.9\Im(\omega_{0})$ $2.512$ $2.519$ $\approx
12.0\Im(\omega_{0})$ Fast $11.850$ $11.237$ $11.781$ $11.854$
$\Im(\omega_{u})-\Im(\omega_{2})$ $\omega_{3}$ $11.871$ $12.152$ $12.311$
$12.313$ $\Im(\omega_{u})-\Im(\omega_{1})$ $\omega_{4}$ $12.715$ $12.715$
$12.718$ $12.772$ $12.790$ $12.758$ $12.771$ $\Im(\omega_{u})-\Im(\omega_{0})$
$\omega_{5}$ $13.279$ $13.282$ $13.282$ $13.232$ $13.263$ $13.238$ $13.229$
$\Im(\omega_{u})+\Im(\omega_{0})$ $\omega_{6}$ $13.694$ $13.847$ $13.693$
$13.687$ $\Im(\omega_{u})+\Im(\omega_{1})$ $14.153$ $14.108$ $14.145$
$\Im(\omega_{u})+\Im(\omega_{2})$ $26.031$ Higher-order interaction $5.404$,
$3.567$ Spurious

Table 2: Imaginary part of DMD eigenvalues (frequencies) of the total magnetic
field $B_{s}=B+b$ from the dynamo models with and without imposed
$\alpha$-effect, obtained with standard and Hankel DMD for different values of
rank and delay. The parameters are specified in the header; the frequencies
were rounded up to the third digit. Only frequencies of peaks with the
amplitude at least 0.001 of the dominant one are shown. FFT frequencies in the
dynamo with imposed $\alpha$-effect are the same as those recovered with
Hankel DMD with $r=14$, $d=9$.

As before, we perform dimension reduction of the mixed magnetic field $B_{s}$
with POD. The largest singular values still come in pairs, as $B_{s}$ visibly
retains some degree of symmetry. Circles in figure 4(c) denote the cumulative
contribution of the first $i$ POD modes to the full data set, both in terms of
$\sigma$ and $\sigma^{2}$. As expected, the cumulative energy of the modes
$\sum_{i}\sigma_{i}^{2}$ grows faster than $\sum_{i}\sigma_{i}$: only 5 modes
are necessary to reproduce 99% of energy in the flow, compared to 14 according
to $\sigma$-criterion. However, with such a small $r$ only the principal
dynamo wave with $\omega_{0}$ is captured correctly with DMD and spurious
eigenvalues appear. We thus perform DMD on the POD basis of rank $r=14$.

Two mode clusters are detected in the DMD spectrum of $B_{s}$: low-frequency
modes $\omega_{0-2}$ with $\Im(\omega)<3$ and high-frequency modes
$\omega_{3-6}$ with $\Im(\omega)>10$ (figure 4b). The mode $\psi_{0}$ with the
slowest oscillation frequency $\Im(\omega_{0})\approx 0.282$ corresponds to
the dynamo wave, with a shape similar to the Parker’s dynamo wave (figure 2c).
It is the large-scale oscillation visible in figure 3. The other frequencies
of the slow-evolution cluster, $\Im(\omega_{1})\approx\pm 3\Im(\omega_{0})$,
$\Im(\omega_{2})\approx\pm 5\Im(\omega_{0})$, correspond to the nonlinear
interactions of the mode $\psi_{0}$ with itself through the nonlinearities
$b^{3}$, $B^{3}$ in (18-19), as shown in section 3.2. The weaker frequency
component $\Im(\omega_{2})$ appears through secondary interaction of
$\omega_{0}$ and $\omega_{1}$ in this nonlinearity, consider, for example
$[\sin(\omega_{0})+\sin(3\omega_{0})]^{3}$. Despite being dynamically related,
the characteristic spatial length scale of the modes $\psi_{1}$, $\psi_{2}$ is
smaller then the length scale of the large-scale dynamo wave $\psi_{0}$
(figure 4d). The eigenvalues in the second, “fast", cluster are located around
the forcing frequency $\omega_{u}=13$, but none of them coincides with it
exactly. Instead, they are result of interaction of the modes in the slow-
evolution cluster, $\psi_{0,1}$, with the forcing $u$. Indeed, the nonlinear
terms in model (17-19) imply $b_{t}\propto
u_{x}B\propto\cos(\pm\omega_{u}t)\exp(\omega_{0,1}t)\propto\pm\omega_{u}\pm\Im(\omega_{0,1})$.
The time-dependence of the DMD modes is exponential according to (8) and (11),
so their interactions will contain all possible combinations of sine and
cosine functions. Length scales of these faster oscillating modes are much
smaller than length scales of the slow modes (see figure 4d).

Comparison of DMD frequencies and amplitudes with FFT results in figure 4(b)
allows to evaluate DMD performance. The optimal amplitudes of each DMD mode
were calculated using the best-fit of the data matrix $Q$ onto the DMD modes
Jovanović et al. (2014). Although standard DMD of rank $r=14$ captures both
fast and slow frequency clusters, the agreement with FFT is imperfect. The DMD
amplitudes of the modes $\omega_{2-6}$ are larger compared to the FFT
amplitudes; weaker DMD eigenvalues of the “fast" cluster, $\omega_{3}$ and
$\omega_{6}$, are shifted towards slower frequencies; these less energetic
modes are captured by standard DMD with less precision.

To improve DMD accuracy, we use the method of delays (Hankel DMD), as
described in section 2.3. We increase the delay parameter up to $d=9$, and
keep the rank $r=14$ to compare with the standard approximation. Hankel DMD
identifies the same dynamical components of $B_{s}$ as standard DMD; however,
the weak dynamical components with frequencies $\omega_{2}$, $\omega_{3}$ and
$\omega_{6}$ are now in a better agreement with FFT results (figure 4b). The
real parts of both weak and strong modes become considerably more neutral (see
figure 6a). Combined with reducing rank $r$, Hankel DMD leads to more compact
yet physical representation of the system dynamics. Consider a reduced-rank
DMD decomposition, $r=6$, of the augmented Hankel matrix $Q_{H}$ with delay
$d=9$. Although only one additional degree of freedom was allowed compared to
$r=5$ from $\sigma^{2}$ cut-off, the principal dynamical modes $\omega_{0}$,
$\omega_{4}$, $\omega_{5}$ with the largest amplitudes are identified
correctly (figure 4b). See table 2 for a list of detected DMD eigenvalues and
Appendix A for further discussion on DMD accuracy as a function of rank and
delay.

###  4.3 DMD of the dynamo without imposed $\alpha$-effect

After identifying the dominant interactions in the system with a well-defined
$\alpha$-effect, it becomes straightforward to perform the same analysis on
the original multiscale dynamo system (17-21). We analyse the modified dynamo
waves in figure 3(a), discarding transients. As mentioned before, FFT of these
data (figure 4a) also results in a “slow" frequency cluster around the dynamo
wave with $\Im(\omega^{\prime}_{0})=0.209$ and the “fast" one around
$\omega_{u}$. Here primes are used to distinguish these frequencies from those
with imposed $\alpha$-effect in section 4.2. Analysis of the dominant
frequency peaks of this spectrum in table 2 shows that the “fast" cluster
corresponds to quadratic interactions of $\omega_{u}$ with magnetic field.
Without a priori defined large-scale $\alpha$-effect, the spectrum of the
“slow" modes becomes more broad-band; nevertheless, principal interactions of
the dynamo wave with itself corresponding to $\Im(3\omega^{\prime}_{0})$,
$\Im(5\omega^{\prime}_{0})$, are still discernible. Higher frequencies of this
cluster, $\Im(9.9\omega^{\prime}_{0})$ , $\Im(12\omega^{\prime}_{0})$ are
potentially influenced by interactions through quadratic terms, but their
magnitudes are weak. After these interactions, the energy is further
transferred downscale, and so additional weaker frequency peaks appear around
$2\omega_{u}=26$ (table 2). As this frequency cluster has much lower energy,
it was not shown in figure 4a.

The singular value spectra of this system are similar to those with imposed
$\alpha$-effect (solid lines in figure 4c), with the difference that the first
pair of POD modes contains more energy, as the amplitude of the large-scale
field is stronger in this case (figure 3a,c). It is thus sufficient to use $4$
POD modes according to $\sigma^{2}$\- and $14$ modes according to
$\sigma$-criterion for subsequent DMD decomposition. In figure 4(a), we
compare resulting frequencies of DMD modes with the FFT spectrum of the field.
In this case, standard DMD decomposition with $r=14$ according to
$\sigma$-criterion reliably identifies only the large-scale dynamo wave with
$\Im(\omega_{0})$ and two strongest modal pairs of its quadratic interaction
with the flow, $\Im(\omega_{4})$, $\Im(\omega_{5})$. The rest of the DMD
eigenvalues either approximate several dynamical components of the spectrum
when located in between dominant peaks, or do not correspond to energetic flow
structures and thus are spurious (table 2). In this case, using Hankel DMD is
crucial to improve scale detection. Using delay of $d=20$ allows to resolve
the “slow" frequency $5\Im(\omega_{0})$ and its interactions with the flow
from the “slow" cluster, as well as to approximate second-order quadratic
interactions proportional to $2\omega_{u}$. As previously, Hankel DMD with
$d=20$ allows to use low, $\sigma^{2}$ model rank $r=4$ and to capture both
fast and slow dynamics using only two pairs of complex-conjugate modes.
Finally, increasing further the values of rank and delay, to $r=40$ and
$d=300$, allows to detect all energetic spectral peaks in figure 4(a).

Spatial shapes of the modes in the “fast" and the “slow" clusters highlight
some differences in dynamics of the two systems. In the system without imposed
$\alpha$-effect, DMD modes are slightly more symmetric with respect to
$x/L=0.5$, as $b$ has much less pronounced systematic time-periodic component
(figure 3). The first two complex-conjugate DMD modes with the dominant
frequency $\Im(\omega_{0})$ are able to reproduce the spatial modulation in
figure 3(a), indicating that it is related to the spatial shape of filtered
$\alpha$-effect. As the nonlinearity $\gamma b^{3}$ was set to zero in this
model, the “slow" modes $\omega_{0}$, $\omega_{1}$ and $\omega_{2}$ remain
large-scale (figure 4d), while the spatial distribution of small-scale modes
from the “fast" cluster is nearly identical both dynamo systems. Except for
these features, the systems with and without imposed $\alpha$-effect have
qualitatively similar dynamics.

##  5 Validity of linear approximation

(a)

(b)

(c)

(d)

(e)

Figure 5: Temporal evolution of DMD modes as a function of time. Solid,
instantaneous (actual) coefficients $c_{i}^{I}$; dashed, coefficients of the
linear model $c_{i}^{L}$. Only real part of $c_{i}$ is shown. (a) $c_{0}$,
dynamo wave, standard DMD with $r=14$. (b) $c_{2}$, corresponding to the mode
$\psi_{2}$, standard DMD with $r=14$. (c) $c_{6}$, corresponding to the mode
$\psi_{6}$, Hankel DMD, $r=14$, $d=9$. (d) $c_{2}$, Hankel DMD, $r=14$, $d=9$.
(e) $c_{2}$, Hankel DMD, $r=22$, $d=10$ (dashed, in red); $c_{2}^{I}$ obtained
with standard DMD is given in blue for comparison. See table 2 and figures
4(b-d) for mode description.

In the following, we analyse how Hankel DMD improves scale detection on the
data from multi-scale dynamo system with imposed $\alpha$-effect (17-20),
(22-23). As this system has fewer dynamical components (figure 4b), it allows
to track easily the influence of varying parameters of DMD model on detection
and properties of individual modes without increasing computational
complexity. We leave the corresponding analysis of the multi-scale dynamo
without imposed $\alpha$-effect for the future work.

DMD approximates the flow dynamics with a linear system; however, the actual
temporal evolution of DMD modes can be more complex than the one prescribed by
(8). According to (8) and (11), temporal coefficients $c_{i}^{L}(t)$ of the
DMD modes $\psi_{i}$ are defined by DMD eigenvalues through
$c_{i}^{L}(t)=\exp(\omega_{i}t)b_{i}(0)$. This linear model allows either
growing or decaying in time oscillating modes. However, we can also evaluate
actual instantaneous temporal coefficients of the modes, $c_{i}^{I}(t)$, from
the data. This is done using non-orthogonal projection of the data matrix $Q$
on the matrix of DMD modes $\Psi$,

$\bm{c}^{I}(t)=\Psi^{\dagger}Q.$ (24)

where $\Psi^{\dagger}$ denotes its pseudoinverse. By comparing $c_{i}^{I}(t)$
and $c_{i}^{L}(t)$, we can assess how linear is the dynamics of the $i$-th DMD
mode.

Figure 5(a) compares $c_{0}^{I}(t)$ and $c_{0}^{L}(t)$ of the dynamo wave mode
$\psi_{0}$. Excellent agreement between the two indicates that the dynamo wave
component evolves essentially linearly. The energetic modes from the “fast"
cluster, $\psi_{4}$ and $\psi_{5}$, also feature nearly linear dynamics but
faster frequencies (not shown here). These modes are represented relatively
accurately both by DMD and Hankel DMD. On the other hand, the actual temporal
dynamics of weaker modes $\psi_{2}$, $\psi_{3}$ and $\psi_{6}$ in both “slow"
and “fast" clusters is not entirely linear. In figure 5(b) we illustrate this
on the temporal coefficient $c_{2}^{I}(t)$ of the mode with $\omega_{2}$. This
mode, besides the principal oscillating frequency, exhibits a slower
modulation on the time scale comparable to the one of $\psi_{0}$. When
$|c_{2}^{I}(t)|$ is large, it is contaminated by a faster evolving signal. The
corresponding linear coefficient $c_{2}^{L}(t)$ is identified by the standard
DMD as decaying (see negative real parts of the corresponding DMD eigenvalues
in figure 6a), and rapidly becomes out-of-phase with the actual coefficient
$c_{2}^{I}(t)$. The $\psi_{2}$ component is thus quickly lost in the linear
DMD model.

Callaham et al. (2022) showed that POD and DMD modes of a flow with two
principal frequencies are nonlinearly correlated through triadic interactions,
creating mixed-frequency content. In our augmented system, the nonlinearity is
more complex, with both quadratic and cubic terms. Quantitative estimation of
these correlations is thus out of scope of this work. Nevertheless, we
identify the presence of mixed frequency dynamics from the temporal
coefficients qualitatively. From equations (18-19), the mode with eigenvalue
$\omega_{2}$ appears through cubic interaction between the dynamo wave
$\omega_{0}$ and its first nonlinearity with
$\Im(\omega_{1})=3\Im(\omega_{0})$ as
$[\exp(\omega_{0})+\exp(3\omega_{0})]^{3}$. When expanded, this gives the
following set of harmonics:
$[\exp(3\omega_{0})+3\exp(5\omega_{0})+3\exp(7\omega_{0})+\exp(9\omega_{0})]$.
In the FFT of $c_{2}^{I}(t)$ (not shown for brevity), all these harmonics are
present, and the frequency content of $7\Im(\omega_{0})$ is only 10 times
weaker than the principal signal frequency, $\Im(\omega_{2})\approx
5\Im(\omega_{0})$. The fast-oscillating contribution to the time series of
$c_{2}^{I}(t)$ corresponds to the frequencies of
$\Im(\omega_{u}\pm\omega_{2})$, i.e., interaction of $\psi_{2}$ and the small-
scale part of the flow $u$. Thus, DMD decomposition of rank $r=14$ mixes these
dynamical components in a single one. The weak modes $\psi_{3}$, $\psi_{6}$
from the “fast" modal cluster are also modulated, although less then
$\psi_{2}$ (figure 5c).

We introduced delay in the data with Hankel DMD to see whether this prediction
improves, while keeping the model rank $r=14$. It can influence the results in
two ways: first, by identifying modal frequency and amplitude more precisely;
and second, by obtaining a spatial modal basis that is closer to Koopman
eigenvectors of the system. In our case, the changes in spatial shape of mode
$\psi_{2}(x)$ were minor; but its eigenvalue $\omega_{2}$ was identified as
less dampened (figure 6a). The mode now decays much slower and its linear
temporal evolution represents the real system well at early times (figure 5d).
Eventually, it gets out of phase with the instantaneous modal coefficient,
indicating that the temporal signal of $\psi_{2}$-mode is still mixed-
frequency. Hankel DMD also recovers modes with $\omega_{3}$ and $\omega_{6}$
as nearly non-decaying and in-phase with the actual signal (figure 5c), when
standard DMD would result in rapidly decaying, out-of-phase behaviour similar
to figure 5(b).

Further improvement is achieved by increasing rank $r$ beyond the one given by
$\sigma$-criterion, as shown by high-rank DMD results in figure 4b (red dots).
The rank of $r=20$ and delay $d=6$ are large enough to separate the weakest
dynamical components with $7\Im(\omega_{0})$, $\Im(\omega_{u}\pm\omega_{2})$
from the mode $\psi_{2}$. Then, the temporal coefficient $c_{2}^{I}(t)$
becomes uncorrelated with them (figure 5e), the corresponding mode $\psi_{2}$
behaves linearly and according to (8). With smaller delays, these additional
dynamical components are not recovered and spurious eigenvalues appear
(crossed-out circles in figure 4b).

##  6 Conclusions

In this work, we used Proper Orthogonal Decomposition and Dynamic Mode
Decomposition for scale identification in one-dimensional dynamo models. Our
first benchmark model belongs to a family of $\alpha-\Omega$ dynamos suggested
by Parker (1955), Proctor (2007) and Richardson & Proctor (2010), among
others. We found that POD and DMD modes of this system have similar spatial
shapes; in fact, DMD can be viewed as a rotation of the POD basis (Callaham et
al., 2022). However, complex-conjugate DMD modes give a sparser dynamical
basis for the oscillating dynamo. We identified the modes corresponding to the
dynamo wave, and the modes appearing due to the nonlinear damping term $B^{3}$
in the model. Together, these two DMD modes form an accurate linear model of
the dynamo.

The shortcoming of the benchmark dynamo model is that it has only one
independent dynamo wave component. We constructed two qualitative augmented
dynamo models featuring both the dynamo wave and the small-scale magnetic
field forced by the flow, with and without imposed $\alpha$-effect. The aim of
these reduced models was not to reproduce the mean-field dynamo theory in a
predictive way but rather to create an clear benchmark for DMD analysis in
multiscale dynamos. By choosing the frequency of the flow, we separated
dynamically relevant scales in space and time in both models. DMD eigenvalues
of this dynamo were located in two regions of the complex plane: the “slow"
cluster near the dynamo wave with frequency $\Im(\omega_{0})$, and the “fast"
cluster, centered about the small-scale flow frequency $\omega_{u}$. Analysis
of the frequency distribution suggests that the “slow" cluster results from
the nonlinear interactions of the dynamo wave with itself through the cubic
terms $B^{3}$, $b^{3}$ in equations (18-19), while the “fast" cluster appears
through quadratic nonlinearity proportional to the product of $u$ and $B$ (or
$b$) in equations (17), (19). As nonlinearities in the models with and without
imposed $\alpha$-effect have similar functional form, the temporal dynamics of
both models is qualitatively similar. Note that DMD is not able to separate
$B$ and $b$ from the mixed signal $B_{s}$; instead, it identifies different
scales with the same temporal dynamics as a single mode. Thus, DMD modes from
the “slow" cluster in figure 4(d) are slightly asymmetric and take into
account the systematic large-scale footprint of the dynamo wave in $b$ through
term $u_{x}B$. In the dynamo without imposed $\alpha$-effect, the modes
corresponding to cubic nonlinear interactions remain uncontaminated with small
scales, because cubic nonlinearity for $b$ was suppressed in this model,
$\gamma=0$ compared to $\gamma=10$ in the case with imposed $\alpha$-effect.
The model without imposed $\alpha$-effect thus features a better separation of
spatial scales.

Furthermore, we investigated the influence of the DMD parameters, rank and
delay, on accuracy of scale detection. Fast-scale dynamics is not captured if
a more conservative $\sigma^{2}$-criterion for $r$ is applied. Thus, rank $r$
is important if dynamically relevant modes are not the most energetic ones. We
demonstrated that Hankel (high-order) DMD improves the accuracy of
approximation, recovering the amplitudes and frequencies of weak dynamical
components in the data with more precision. In the dynamo with imposed
$\alpha$-effect, we identified that the linearity assumption breaks for the
weakest DMD modes; Hankel DMD gives a better linear representation for these
modes, especially if the rank of the model is large enough to separate mixed-
frequency components in the these modes. This improvement relies on
spatiotemporal coherency of the corresponding dynamical scales. In dynamical
systems with a more broad-band spectrum, like the one without imposed
$\alpha$-effect, the errors of linear model (8) will be larger; nevertheless,
Hankel DMD still considerably improves accuracy of detection and modelling of
energetic length scales. In this case, further extensions of DMD method are of
interest - for example, the multiresolution DMD based on wavelet-like
techniques proposed by Dylewsky et al. (2019). We leave this analysis for the
future work.

Finally, we give some notes on extension of this approach to more realistic
DNS of shear dynamos. There, the small-scale flow is frequently forced by a
helical forcing with a defined length scale, or forcing wave number, but
randomized in time (Tobias & Cattaneo, 2013; Nigro et al., 2017). In this
case, DMD is expected to represent poorly the small-scale flow and field
components, because they are incoherent in time. Nevertheless, it should be
possible to extract the large-scale, coherent components of the magnetic and
velocity fields corresponding to the “slow" modal cluster from such
simulations. As a final remark, our augmented benchmark dynamo model has a
well-defined scale separation while many turbulent dynamos have continuous
energy spectra. These dynamos are expected to give continuous DMD spectra
rather than clustered.

## Acknowledgements

This project has received funding from the European Union’s Horizon 2020
research and innovation programme under the Marie Skłodowska-Curie grant
agreement No 890847. The author is grateful to Prof. Steven Tobias and Dr.
Calum Skene for discussions that helped to shape many aspects of this work.
Support from the ARC supercomputing center at the University of Leeds, Leeds
Institute for Fluid Dynamics, and the Isaac Newton Institute for Mathematical
Sciences, Cambridge, during the programme “Frontiers in dynamo theory: from
the Earth to the stars", is fully acknowledged. This work was supported by
EPSRC grant no EP/R014604/1.

## Data Availability

The data in this work were generated using publicly available code Dedalus
(Burns et al., 2020) and will be provided on request to the corresponding
author.

## References

  * Arbabi & Mezic (2017) Arbabi H., Mezic I., 2017, SIAM Journal on Applied Dynamical Systems, 16, 2096
  * Aubry et al. (1993) Aubry N., Lian W.-Y., Titi E. S., 1993, SIAM Journal on Scientific Computing, 14, 483
  * Babcock (1961) Babcock H., 1961, The Astrophysical Journal, 133, 572
  * Bagheri (2010) Bagheri S., 2010, PhD thesis, KTH
  * Brandenburg (2018) Brandenburg A., 2018, Journal of Plasma Physics, 84, 735840404
  * Brandenburg & Subramanian (2005) Brandenburg A., Subramanian K., 2005, Physics Reports, 417, 1
  * Broomhead & King (1986) Broomhead D. S., King G. P., 1986, Physica D: Nonlinear Phenomena, 20, 217
  * Burns et al. (2020) Burns K. J., Vasil G. M., Oishi J. S., Lecoanet D., Brown B. P., 2020, Physical Review Research, 2, 023068
  * Bushby (2003) Bushby P., 2003, Monthly Notices of the Royal Astronomical Society, 338, 655
  * Callaham et al. (2022) Callaham J. L., Brunton S. L., Loiseau J.-C., 2022, Journal of Fluid Mechanics, 938
  * Cameron et al. (2018) Cameron R. H., Duvall T., Schüssler M., Schunker H., 2018, Astronomy & Astrophysics, 609, A56
  * Cattaneo & Hughes (1996) Cattaneo F., Hughes D. W., 1996, Physical Review E, 54, R4532
  * Charbonneau (2014) Charbonneau P., 2014, Annual Review of Astronomy and Astrophysics, 52, 251
  * Charbonneau (2020) Charbonneau P., 2020, Living Reviews in Solar Physics, 17, 1
  * Charbonneau & MacGregor (1997) Charbonneau P., MacGregor K., 1997, The Astrophysical Journal, 486, 502
  * Deane et al. (1991) Deane A., Kevrekidis I., Karniadakis G. E., Orszag S., 1991, Physics of Fluids A: Fluid Dynamics, 3, 2337
  * Dylewsky et al. (2019) Dylewsky D., Tao M., Kutz J. N., 2019, Physical Review E, 99, 063311
  * Fujii et al. (2019) Fujii K., Takeishi N., Kibushi B., Kouzaki M., Kawahara Y., 2019, Scientific reports, 9, 16755
  * Gilman (1983) Gilman P. A., 1983, The Astrophysical Journal Supplement Series, 53, 243
  * Holmes et al. (2012) Holmes P., Lumley J. L., Berkooz G., Rowley C. W., 2012, Turbulence, coherent structures, dynamical systems and symmetry. Cambridge university press
  * Hori et al. (2023) Hori K., Jones C. A., Antuñano A., Fletcher L. N., Tobias S. M., 2023, Nature Astronomy, pp 1–11
  * Huang et al. (1998) Huang N. E., et al., 1998, Proceedings of the Royal Society of London. Series A: mathematical, physical and engineering sciences, 454, 903
  * Jones et al. (2010) Jones C. A., Thompson M. J., Tobias S. M., 2010, Space Science Reviews, 152, 591
  * Jovanović et al. (2014) Jovanović M. R., Schmid P. J., Nichols J. W., 2014, Physics of Fluids, 26, 024103
  * Käpylä et al. (2013) Käpylä P. J., Mantere M. J., Cole E., Warnecke J., Brandenburg A., 2013, The Astrophysical Journal, 778, 41
  * Koopman (1931) Koopman B. O., 1931, Proceedings of the National Academy of Sciences, 17, 315
  * Krause & Raedler (1980) Krause F., Raedler K.-H., 1980, Oxford: Pergamon Press
  * Le Clainche & Vega (2017) Le Clainche S., Vega J. M., 2017, SIAM Journal on Applied Dynamical Systems, 16, 882
  * Leighton (1969) Leighton R. B., 1969, The Astrophysical Journal, 156, 1
  * Lumley (2007) Lumley J. L., 2007, Stochastic tools in turbulence. Courier Corporation
  * Moffatt (1978) Moffatt H. K., 1978, Cambridge Monographs on Mechanics and Applied Mathematics
  * Moffatt & Dormy (2019) Moffatt K., Dormy E., 2019, Self-exciting fluid dynamos. Vol. 59, Cambridge University Press
  * Muñoz-Jaramillo & Vaquero (2019) Muñoz-Jaramillo A., Vaquero J. M., 2019, Nature Astronomy, 3, 205
  * Nauman & Pessah (2016) Nauman F., Pessah M. E., 2016, The Astrophysical Journal, 833, 187
  * Nigro & Veltri (2011) Nigro G., Veltri P., 2011, The Astrophysical Journal Letters, 740, L37
  * Nigro et al. (2017) Nigro G., Pongkitiwanichakul P., Cattaneo F., Tobias S., 2017, Monthly Notices of the Royal Astronomical Society: Letters, 464, L119
  * Parker (1955) Parker E. N., 1955, The Astrophysical Journal, 122, 293
  * Parker (1979) Parker E. N., 1979, Cosmical magnetic fields: Their origin and their activity. Oxford University Press
  * Pongkitiwanichakul et al. (2016) Pongkitiwanichakul P., Nigro G., Cattaneo F., Tobias S., 2016, The Astrophysical Journal, 825, 23
  * Press et al. (1992) Press W. H., Vetterling W. T., Teukolsky S. A., Flannery B. P., 1992, Numerical Recipes Example Book (FORTRAN). Cambridge University Press Cambridge
  * Proctor (2007) Proctor M. R., 2007, Monthly Notices of the Royal Astronomical Society: Letters, 382, L39
  * Racine et al. (2011) Racine É., Charbonneau P., Ghizaru M., Bouchat A., Smolarkiewicz P. K., 2011, The Astrophysical Journal, 735, 46
  * Richardson & Proctor (2010) Richardson K., Proctor M., 2010, Geophysical and Astrophysical Fluid Dynamics, 104, 601
  * Rincon (2019) Rincon F., 2019, Journal of Plasma Physics, 85, 205850401
  * Saikia et al. (2016) Saikia S. B., et al., 2016, Astronomy & Astrophysics, 594, A29
  * Schmid (2010) Schmid P. J., 2010, Journal of fluid mechanics, 656, 5
  * Schmid (2022) Schmid P. J., 2022, Annual Review of Fluid Mechanics, 54, 225
  * Schmidt & Colonius (2020) Schmidt O. T., Colonius T., 2020, Aiaa journal, 58, 1023
  * Schrinner et al. (2011) Schrinner M., Petitdemange L., Dormy E., 2011, Astronomy & Astrophysics, 530, A140
  * Sirovich (1987) Sirovich L., 1987, Quarterly of applied mathematics, 45, 561
  * Taira et al. (2020) Taira K., Hemati M. S., Brunton S. L., Sun Y., Duraisamy K., Bagheri S., Dawson S. T., Yeh C.-A., 2020, AIAA journal, 58, 998
  * Takens (1981) Takens F., 1981, in , Dynamical systems and turbulence, Warwick 1980. Springer, pp 366–381
  * Tobias (1997) Tobias S., 1997, Astronomy and Astrophysics, 322, 1007
  * Tobias (2021) Tobias S., 2021, Journal of fluid mechanics, 912
  * Tobias & Cattaneo (2013) Tobias S. M., Cattaneo F., 2013, Nature, 497, 463
  * Tu et al. (2014) Tu J. H., Rowley C. W., Luchtenburg D. M., Brunton S. L., Kutz J. N., 2014, Journal of Computational Dynamics, 1, 391
  * Vainshtein & Cattaneo (1992) Vainshtein S. I., Cattaneo F., 1992, The Astrophysical Journal, 393, 165

##  Appendix A Influence of DMD rank and delay

(a)

(b)

Figure 6:  (a) Spectrum of DMD eigenvalues of the multiscale dynamo model
without imposed $\alpha$-effect; only eigenvalues with positive imaginary part
(frequency) are shown. Filled circles, standard DMD of $r=14$; filled squares,
Hankel DMD with rank $r=14$ and delay $d=9$; empty squares, Hankel DMD, $r=6$,
$d=9$. Dashed, neutral stability line. Superscipt $H$ denotes frequencies
obtained with Hankel DMD. (b) The model error
$||B_{s}-B_{s}^{model}||_{2}/||B_{s}||_{2}$ in this model as a function of
rank and delay. The colours are spaced logarithmically. Symbols represent the
analysis parameters in panel (a). The empty circle corresponds to standard DMD
with $r=5$ according to $\sigma^{2}$-criterion. Solid line,
$\sigma^{2}$-criterion; dotted, $\sigma$-criterion.

In this section, we show how varying the rank and delay influences accuracy of
the linear DMD model, on the example of augmented dynamo with imposed
$\alpha$-effect (17-19), (22-23). As the dynamics in this system is
predominantly periodic, one would expect the DMD modes computed in the data
set to be nearly neutral (i.e. with $\Re(\omega)\approx 0$). However, DMD
eigenvalues $\omega_{2}$, $\omega_{3}$ and $\omega_{6}$, computed with the
standard DMD algorithm with rank $r=14$ according to $\sigma$-criterion, are
dampened and have pronounced negative growth rates (figure 6a). This will
contribute to overall error (12) of the linear DMD approximation (8), since
heavily dampened components will decay rapidly. On the other hand, Hankel DMD
with delay $d=9$ results in all the modes being close to neutral stability,
reducing the error from 12% to 2%.

Overall, there is no analytical expression for optimal $d$ and $r$ which tend
to increase with the data complexity. So far, we performed a point-by-point
analysis in the parameter space of $(r,d)$; it can be generalized by spanning
a range of $r$ and $d$ and tracking error (12) of linear DMD approximation
(8). Figure 6(b) presents this error as a function of both $r$ and $d$, with
the DMD models corresponding to figures 6(a) denoted accordingly. With or
without delay, increasing rank reduces the error of the linear approximation.
Increasing delay improves accuracy at first, but when $d$ is large, the
dependence of error on $d$ saturates. To get a better understanding of this
behaviour, we evaluated how delay affects the rank of Hankel data matrix
$Q_{H}$ (13). We evaluate this in terms of rank $r^{\sigma(\sigma^{2})}$ of
$Q_{H}$ according to either $\sigma$ or $\sigma^{2}$ criteria. These ranks are
denoted in figure 6(b) by dotted and solid lines, respectively. They estimate
the number of independent vector components required to describe certain
percentage of the information contained in $Q_{H}$. At first, their values
increase with $d$, indicating increase in data rank as newly added delayed
system states introduce new dynamical information. However, this trend
saturates around $d=5$. Further increasing delay does not introduce new
dynamical components, and the representation error becomes predominantly a
function of $r$ when $d$ is large. A reasonable minimum bound for $d$
therefore is that the rank of matrix $Q_{H}$, defined either through $\sigma$,
$\sigma^{2}$ criteria or numerical estimation of the matrix rank, should
saturate.

Another insight from figure 6(b) is that increasing delay does not recover the
dynamical components that were entirely removed by truncating the POD modal
basis at rank $r$, as demonstrated by the abrupt decrease in error between
$r=13$ and $r=14$, from $\epsilon\sim 0.1$ to $\epsilon\sim 0.01$. At $r=14$,
the weak DMD harmonic with frequency $\omega_{2}$ first appears, improving
accuracy. When $r<14$, this dynamical component is not identified by DMD,
independently of $d$. Note that FFT of $B_{s}$ shows other, even weaker
frequencies of $\Im(7\omega_{0})$, $\Im(\omega_{u}\pm\omega_{2})$ in figure
4(b), not captured by standard or delayed DMD with $r=14$. As a final note,
using increasingly high delays like those used in figure 4(a) can result in
eigenvalues moving up in the complex plane of figure 6(a) so that
$\Re(\omega)>0$ and therefore some DMD modes become exponentially growing. In
this case, increased accuracy of frequency detection will be accompanied by
higher errors in the linear model (8).