Data-driven scale identification in oscillatory dynamos

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 (DMD), 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 multiscale 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 non-linear self-interaction, and the ‘fast’ cluster, centred 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 non-linear 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.

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 (Saikia et al. 2016; Muñoz-Jaramillo & Vaquero 2019). 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′ and fluctuating field b′ correlate well, so that the mean electromotive force, or emf ε = 〈u′ × b′〉, 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 ε_i = αδ_ijB_j. In a nutshell, the dynamo can be described with the one-dimensional α − Ω model,

where η 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. Ω-effect corresponds to stretching of poloidal magnetic field with potential A into toroidal field B by large-scale differential rotation Ω′; α-effect parametrizes generation of A from B by helical turbulent motions. Depending on the form of α, equations (1) and (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 (equations 1 and 2) is kinematic, i.e. magnetic field grows linearly from a prescribed flow parametrized by Ω′ and α. 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 (Racine et al. 2011; Schrinner, Petitdemange & Dormy 2011; Käpylä et al. 2013).

Symmetry-preserving non-linear term B³ was introduced into equation (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 (α-effect) when magnetic energy becomes comparable to kinetic energy of the flow (e.g. Jones, Thompson & Tobias 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 non-linear interactions, it is necessary to unambiguously identify coherent dynamo structures in multiscale dynamo models like state-of-art DNS (with about 10⁶ 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 non-linear 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 (Sirovich 1987; Lumley 2007; 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 (Bagheri 2010; Schmid 2010, 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 DMD, 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 α−Ω model (equations 1 and 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

Consider the original non-linear system (equations 1 and 2) or an analogous one that evolves the magnetic potential A and the magnetic field B and depends on a set |$\boldsymbol {\mu }$| of dimensionless functions and parameters,

The functional form of the α-effect and the boundary conditions are

Dynamo number D = Ω′α₀/η²l³ is the dimensionless parameter of the system defining frequency of waves; if η = l = Ω′ = 1, then D = α₀. We solve system (equation 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 travelling waves for D > 20, as shown in Fig. 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, |$\boldsymbol {q}_k (x) \in R^N$|⁠, at time steps k = 1, 2, …, K, with frequency sufficient to resolve the fastest flow dynamics, and assemble them into the data matrix Q,

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

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 q, defined on the domain x, POD modes are eigenfunctions of the integral equation

with autocorrelation function 〈q(x)q*(x′)〉 as a kernel, and 〈…〉 denoting temporal average. The eigenvalues σ² ≥ 0 correspond to the average energy content in each mode ϕ, and the first r modes represent more of the system energy than 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 (equation 5), modes ϕ are eigenvectors of positive semi-definite correlation matrix Q*Q. In practice, it is convenient to find ϕ by factorizing the data matrix Q with singular value decomposition (SVD),

which generalizes the eigenvalue decomposition to non-square matrices. The data set is thus represented as a product of the matrix Φ of the POD modes ϕ_i, the matrix Σ of their singular values σ_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 (equations 6 and 7) (see Schmidt & Colonius 2020).

2.2 Dynamic mode decomposition

Unlike statistics-based POD, DMD harnesses temporal dynamical behaviour of the flow. It is closely related to Koopman operator theory, proposing that the temporal dynamics of non-linear system (equation 3) can be described with a linear operator acting on suitable flow observables (Koopman 1931). For an arbitrary non-linear 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 non-linear system (equation 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},

where b_i(0) is the initial magnitude of the eigenvector ψ_i. Discrete eigenvalues λ_i define whether ψ_i are decaying, neutral, or growing in time. We augment the data matrix (equation 5) with one more snapshot of the system, and construct another matrix Q′, such that

A least-squares approximation to |$\mathcal {A}$| can be obtained by calculating the matrix product of the pseudo-inverse of Q with Q′. However, for low-rank matrices like Q, pseudo-inverse 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 SVD (equation 7) at rank r. After this operation, equation (9) is re-arranged as

where |$\mathcal {A}_r$| is a projection of |$\mathcal {A}$| on the POD modes Φ_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

The second equation in equation (11) relates discrete-time eigenvalues λ to continuous-time eigenvalues ω, both complex. In the following, we will refer to ψ and ω as DMD modes and DMD eigenvalues. If the real part of the eigenvalue ℜ(ω) < 0, the mode is dampened; if ℜ(ω) > 0, the mode will be growing. In a steady-state system, its significant modes are expected to be nearly neutral, ℜ(ω) ≈ 0. The imaginary part ℑ(ω) 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 ∑_rσ_r/∑_totσ = 0.99, retaining 99 per cent 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 σ², and not σ, represents the energy of POD modes in equation (6). The latter criterion means that 99 per cent of energy is retained and results in a sparser spatial basis. The relative error of the DMD model can be estimated as

where || · || is L₂-norm and Q^model is the approximation to Q calculated according to equation (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 (equation 5) as

where d is the delay parameter, and augment (equation 9) in the same way. This is done because our initial choice of system observables, |$\boldsymbol {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 than the rank of Q. Arbabi & Mezic (2017) proved that in the limit d → ∞ DMD eigenvalues and eigenmodes approximating the system,

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 equation (9) (Le Clainche & Vega 2017).

Construction of equation (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 equation (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₂ norm (Fujii et al. 2019).

3 PARKER’S DYNAMO WAVES

As a benchmark problem, we take the system (equations 1 and 2) augmented with (equation 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 Fig. 1.

3.1 Proper orthogonal decomposition

Fig. 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 per cent of the data set according to σ criterion, compared to 2 modes for A. This difference can be attributed to the non-linear term B³ in equation (2), and therefore more complex dynamics. Potential A is influenced by this non-linearity only indirectly through α effect. According to σ² criterion, 99 per cent 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

If the solution of equation (3) is invariant under the action of symmetry group G, and ϕ(x) is an eigenfunction of equation (6), then for every element g ∈ G such as equation (15), g°ϕ is also an eigenfunction of equation (6) with the same eigenvalue (Aubry, Lian & Titi 1993, Proposition 3).

The POD modes, corresponding to the singular value pairs, are shown in Fig. 2(b). Following the fields parity in Fig. 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 Fig. 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 (Fig. 2c, upper and middle panel). They enhance each other at some spatial locations, and weaken at others, producing travelling wave dynamics observed in Fig. 2(c), bottom panel. The relative error between this two-mode model and toroidal field B in Fig. 1, evaluated as L₂ norm, is only 3.7 per cent.

3.2 Dynamic mode decomposition

We calculate the DMD eigenmodes and eigenvalues of the dynamo waves in Fig. 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 equation (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, ℑ(ω_{0, 1}) = ±3.558, corresponds to the frequency of dynamo linear instability at D = 30, and is the main frequency observed in Fig. 1. The second frequency is three times larger than the first one, ℑ(ω_{2, 3}) ≈ 3ℑ(ω_{0, 1}), and can be interpreted as an interaction of the first mode with itself through non-linear damping B³ in equation (2), as the periodic temporal behaviour implies sin ³(ω) ∼ sin (3ω). Note the implication of choosing σ or σ² cut-off: if only 2 POD modes are retained according to σ² 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 Fig. 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 Fig. 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 Fig. 2(c) with just one complex DMD mode ψ₀ and its temporal coefficient. Indeed, consider a DMD mode ψ₀(x) = ψ^r + ψⁱi, and the corresponding projection of the data set on this mode, c₀(t) = c^r + cⁱi. One-mode approximation for the data matrix Q would be

When the corresponding contribution from the conjugate mode ψ₁ will be added, the imaginary part of equation (16) will be cancelled, as the data matrix Q is real; the real part will be double of that of equation (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.

4 MULTISCALE DYNAMO MODEL

Parker’s dynamo model (equations 1 and 2) parametrizes interaction of small-scale dynamo and turbulence with α-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  per cent of the total magnetic energy, as can be seen in Fig. 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 (equations 1 and 2) with extra terms, mimicking the multiscale 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

Temporal and spatial evolution of the small-scale flow field u is prescribed by equation (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 non-linear in fluctuations due to the terms |$\nabla \times (u^{\prime } \times b^{\prime }-\overline{u^{\prime } \times b^{\prime }})$|⁠, non-linearity 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 non-linearity γb³ into equation (19), which allows to control the saturation amplitude of b. The dissipative term in equation (19) is proportional to b_xx; μu_xB models generation of b through interaction of the flow and the large-scale field. Parameters γ and μ are introduced to control the strength of non-linearity 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 ε = ε₀g(ub) on the scale of the dynamo waves in equation (17) is calculated by filtering product ub with a Gaussian filter (equation 20) and ensuring ε = 0 at the boundaries x = 0, L. Alternatively, time averaging or envelope smoothing can be employed to model action of emf on large scales.

Dynamo waves were consistently observed when ε₀ > 800 and the rest of parameters set to Ω′ = η = l = 1, μ = 5, u₀ = 0.5. We set γ = 0, allowing the dynamo to saturate only through large scales. Large value of σ = 30 removes completely small-scale fluctuations from ε and thus inhibits non-oscillatory dynamo solutions (Richardson & Proctor 2010). Frequency ω_u and wave number k_u of small-scale u were set to ω_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. Fig. 3(a) and (c) shows results of numerical integration of equations (17)–(21) with these parameters and ε₀ = 1500. The large-scale field B is again dipolar, antisymmetric 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 fig. 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 (Fig. 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 ω_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 α-effect as in equation (4).

4.1 Imposing α-effect

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

so that magnetic fluctuations scale as b ∼ cos (πx/L)B. If electromotive force is defined as

then α-effect is similar to equation (4), yet dynamo waves appear from interaction of the flow and small-scale magnetic field. Formally, the extra component in equation (22) can be interpreted as a large-scale one, although it is not constant in x compared to mean differential rotation Ω′ in equation (18). Magnetic fluctuations will thus contain direct interaction of this component in equation (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 (equations 17–21) without imposed α-effect, where large-scale shear interacts with B indirectly through Ω′A_x and u_xB. Controlling α-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 (equations 1–2).

We perform a numerical simulation of augmented model (equations 17–19), together with (equations 22 and 23), with dipolar initial conditions. Model parameters were set to ε₀ = 30, Ω′ = l = μ = 1, η = 0.1, u₀ = 0.5, σ = 10, ω_u = 13, k_u = 13, resulting in ‘classic’ large-scale dynamo waves with a period T ≈ 20, about ten times longer (Fig. 3b). The small-scale non-linearity was set non-zero, γ = 10, in order to reduce the magnitude of b to values below those of large-scale B (Fig. 3d), as in dynamo without imposed α-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 antisymmetric u_x and antisymmetric B. Note the difference with multiscale dynamo system (equations 17–21): the sign of the large-scale footprint in b is periodic in time, following the sign of u_xB ∼ cos (πx/L)B. This is less apparent in Fig. 3(c) because the product of B and u_x remains small scale. The resulting multiscale 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 (Fig. 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 α-effect

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 Fig. 4(c) denote the cumulative contribution of the first i POD modes to the full data set, both in terms of σ and σ². As expected, the cumulative energy of the modes |$\sum _i \sigma _i^2$| grows faster than ∑_iσ_i: only 5 modes are necessary to reproduce 99  per cent of energy in the flow, compared to 14 according to σ-criterion. However, with such a small r only the principal dynamo wave with ω₀ 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 ω₀₋₂ with ℑ(ω) < 3 and high-frequency modes ω₃₋₆ with ℑ(ω) > 10 (Fig. 4b). The mode ψ₀ with the slowest oscillation frequency ℑ(ω₀) ≈ 0.282 corresponds to the dynamo wave, with a shape similar to the Parker’s dynamo wave (Fig. 2c). It is the large-scale oscillation visible in Fig. 3. The other frequencies of the slow-evolution cluster, ℑ(ω₁) ≈ ±3ℑ(ω₀), ℑ(ω₂) ≈ ±5ℑ(ω₀), correspond to the non-linear interactions of the mode ψ₀ with itself through the non-linearities b³, B³ in equations (18) and (19), as shown in Section 3.2. The weaker frequency component ℑ(ω₂) appears through secondary interaction of ω₀ and ω₁ in this non-linearity, consider, for example [sin (ω₀) + sin (3ω₀)]³. Despite being dynamically related, the characteristic spatial length scale of the modes ψ₁, ψ₂ is smaller than the length scale of the large-scale dynamo wave ψ₀ (Fig. 4d). The eigenvalues in the second, ‘fast’, cluster are located around the forcing frequency ω_u = 13, but none of them coincides with it exactly. Instead, they are result of interaction of the modes in the slow-evolution cluster, ψ_0,1, with the forcing u. Indeed, the non-linear terms in model (equations 17 and 19) imply b_t ∝ u_xB ∝ cos (± ω_ut)exp (ω_{0, 1}t) ∝ ± ω_u ± ℑ(ω_{0, 1}). The time dependence of the DMD modes is exponential according to equations (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 Fig. 4d).

Comparison of DMD frequencies and amplitudes with FFT results in Fig. 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 on to the DMD modes (Jovanović, Schmid & Nichols 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 ω₂₋₆ are larger compared to the FFT amplitudes; weaker DMD eigenvalues of the ‘fast’ cluster, ω₃ and ω₆, 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 ω₂, ω₃, and ω₆ are now in a better agreement with FFT results (Fig. 4b). The real parts of both weak and strong modes become considerably more neutral (see Fig. A1a). 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 σ² cut-off, the principal dynamical modes ω₀, ω₄, ω₅ with the largest amplitudes are identified correctly (Fig. 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 α-effect

After identifying the dominant interactions in the system with a well-defined α-effect, it becomes straightforward to perform the same analysis on the original multiscale dynamo system (equations 17–21). We analyse the modified dynamo waves in Fig. 3(a), discarding transients. As mentioned before, FFT of these data (Fig. 4a) also results in a ‘slow’ frequency cluster around the dynamo wave with |$\Im (\omega ^{\prime }_0) = 0.209$| and the ‘fast’ one around ω_u. Here, primes are used to distinguish these frequencies from those with imposed α-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 ω_u with magnetic field. Without a priori defined large-scale α-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ω_u = 26 (Table 2). As this frequency cluster has much lower energy, it was not shown in Fig. 4(a).

The singular value spectra of this system are similar to those with imposed α-effect (solid lines in Fig. 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 (Fig. 3a and c). It is thus sufficient to use 4 POD modes according to σ² and 14 modes according to σ-criterion for subsequent DMD decomposition. In Fig. 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 σ-criterion reliably identifies only the large-scale dynamo wave with ℑ(ω₀) and two strongest modal pairs of its quadratic interaction with the flow, ℑ(ω₄), ℑ(ω₅). 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ℑ(ω₀) and its interactions with the flow from the ‘slow’ cluster, as well as to approximate second-order quadratic interactions proportional to 2ω_u. As previously, Hankel DMD with d = 20 allows to use low, σ² 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 Fig. 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 α-effect, DMD modes are slightly more symmetric with respect to x/L = 0.5, as b has much less pronounced systematic time-periodic component (Fig. 3). The first two complex–conjugate DMD modes with the dominant frequency ℑ(ω₀) are able to reproduce the spatial modulation in Fig. 3(a), indicating that it is related to the spatial shape of filtered α-effect. As the non-linearity γb³ was set to zero in this model, the ‘slow’ modes ω₀, ω₁, and ω₂ remain large scale (Fig. 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 α-effect have qualitatively similar dynamics.

5 VALIDITY OF LINEAR APPROXIMATION

In the following, we analyse how Hankel DMD improves scale detection on the data from multiscale dynamo system with imposed α-effect (equations 17–20 and 22–23). As this system has fewer dynamical components (Fig. 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 multiscale dynamo without imposed α-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 equation (8). According to equations (8) and (11), temporal coefficients |$c_i^L(t)$| of the DMD modes ψ_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 Ψ,

where Ψ^† denotes its pseudo-inverse. By comparing |$c_i^I(t)$| and |$c_i^L(t)$|⁠, we can assess how linear is the dynamics of the ith DMD mode.

Fig. 5(a) compares |$c_0^I(t)$| and |$c_0^L(t)$| of the dynamo wave mode ψ₀. Excellent agreement between the two indicates that the dynamo wave component evolves essentially linearly. The energetic modes from the ‘fast’ cluster, ψ₄ and ψ₅, 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 ψ₂, ψ₃, and ψ₆ in both ‘slow’ and ‘fast’ clusters is not entirely linear. In Fig. 5(b), we illustrate this on the temporal coefficient |$c_2^I(t)$| of the mode with ω₂. This mode, besides the principal oscillating frequency, exhibits a slower modulation on the time scale comparable to the one of ψ₀. 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 Fig. A1a), and rapidly becomes out-of-phase with the actual coefficient |$c_2^I(t)$|⁠. The ψ₂ component is thus quickly lost in the linear DMD model.

Callaham, Brunton & Loiseau (2022) showed that POD and DMD modes of a flow with two principal frequencies are non-linearly correlated through triadic interactions, creating mixed-frequency content. In our augmented system, the non-linearity 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) and (19), the mode with eigenvalue ω₂ appears through cubic interaction between the dynamo wave ω₀ and its first non-linearity with ℑ(ω₁) = 3ℑ(ω₀) as [exp (ω₀) + exp (3ω₀)]³. When expanded, this gives the following set of harmonics: [exp (3ω₀) + 3exp (5ω₀) + 3exp (7ω₀) + exp (9ω₀)]. In the FFT of |$c_2^I(t)$| (not shown for brevity), all these harmonics are present, and the frequency content of 7ℑ(ω₀) is only 10 times weaker than the principal signal frequency, ℑ(ω₂) ≈ 5ℑ(ω₀). The fast-oscillating contribution to the time series of |$c_2^I(t)$| corresponds to the frequencies of ℑ(ω_u ± ω₂), i.e. interaction of ψ₂ 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 ψ₃, ψ₆ from the ‘fast’ modal cluster are also modulated, although less then ψ₂ (Fig. 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 ψ₂(x) were minor; but its eigenvalue ω₂ was identified as less dampened (Fig. A1a). The mode now decays much slower and its linear temporal evolution represents the real system well at early times (Fig. 5d). Eventually, it gets out of phase with the instantaneous modal coefficient, indicating that the temporal signal of ψ₂-mode is still mixed-frequency. Hankel DMD also recovers modes with ω₃ and ω₆ as nearly non-decaying and in-phase with the actual signal (Fig. 5c), when standard DMD would result in rapidly decaying, out-of-phase behaviour similar to Fig. 5(b).

Further improvement is achieved by increasing rank r beyond the one given by σ-criterion, as shown by high-rank DMD results in Fig. 4(b) (red dots). The rank of r = 20 and delay d = 6 are large enough to separate the weakest dynamical components with 7ℑ(ω₀), ℑ(ω_u ± ω₂) from the mode ψ₂. Then, the temporal coefficient |$c_2^I(t)$| becomes uncorrelated with them (Fig. 5e), the corresponding mode ψ₂ behaves linearly and according to equation (8). With smaller delays, these additional dynamical components are not recovered and spurious eigenvalues appear (crossed-out circles in Fig. 4b).

6 CONCLUSIONS

In this work, we used POD and DMD for scale identification in one-dimensional dynamo models. Our first benchmark model belongs to a family of α − Ω 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, Brunton & Loiseau 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 non-linear damping term B³ 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 α-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 ℑ(ω₀), and the ‘fast’ cluster, centred about the small-scale flow frequency ω_u. Analysis of the frequency distribution suggests that the ‘slow’ cluster results from the non-linear interactions of the dynamo wave with itself through the cubic terms B³, b³ in equations (18) and (19), while the ‘fast’ cluster appears through quadratic non-linearity proportional to the product of u and B (or b) in equations (17) and (19). As non-linearities in the models with and without imposed α-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 Fig. 4(d) are slightly asymmetric and take into account the systematic large-scale footprint of the dynamo wave in b through term u_xB. In the dynamo without imposed α-effect, the modes corresponding to cubic non-linear interactions remain uncontaminated with small scales, because cubic non-linearity for b was suppressed in this model, γ = 0 compared to γ = 10 in the case with imposed α-effect. The model without imposed α-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 σ²-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 α-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 α-effect, the errors of linear model (equation 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, Tao & Kutz (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 centre 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

APPENDIX A: INFLUENCE OF DMD RANK AND DELAY

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 α-effect (equations 17–19 and 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 ℜ(ω) ≈ 0). However, DMD eigenvalues ω₂, ω₃, and ω₆, computed with the standard DMD algorithm with rank r = 14 according to σ-criterion, are dampened and have pronounced negative growth rates (Fig. A1a). This will contribute to overall error (equation 12) of the linear DMD approximation (equation 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 per cent.

Overall, there is no analytical expression for optimal d and r that 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 (equation 12) of linear DMD approximation (equation 8). Fig. A1(b) presents this error as a function of both r and d, with the DMD models corresponding to Fig. A1(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 (equation 13). We evaluate this in terms of rank |$r^{\sigma (\sigma ^2)}$| of Q_H according to either σ or σ² criteria. These ranks are denoted in Fig. A1(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 σ, σ² criteria or numerical estimation of the matrix rank, should saturate.

Another insight from Fig. A1(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 ϵ ∼ 0.1 to ϵ ∼ 0.01. At r = 14, the weak DMD harmonic with frequency ω₂ 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 ℑ(7ω₀), ℑ(ω_u ± ω₂) in Fig. 4(b), not captured by standard or delayed DMD with r = 14. As a final note, using increasingly high delays like those used in Fig. 4(a) can result in eigenvalues moving up in the complex plane of Fig. A1(a) so that ℜ(ω) > 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 (equation 8).