2.1. Equations and definitions

The equations for MHD in the incompressible case are written as usual as

(2.1) \begin{align} [ \partial _t + \boldsymbol{u} \cdot \nabla ]\boldsymbol{u} \equiv D_t \boldsymbol{u} &= -\nabla P^\star +\boldsymbol{b} \cdot \nabla \boldsymbol{b} + \nu \Delta \boldsymbol{u} +\boldsymbol{F}_V, \\[-4pt] \nonumber \end{align}

(2.2) \begin{align} [\partial _t + \boldsymbol{u} \cdot \nabla ] \boldsymbol{b} \equiv D_t \boldsymbol{b} &= \boldsymbol{b} \cdot \nabla \boldsymbol{u} + \eta \Delta \boldsymbol{b}, \\[8pt]\nonumber\end{align}

together with $\nabla \cdot \boldsymbol{u}=0, \nabla \cdot \boldsymbol{b}=0$ ; here, $\boldsymbol{u}, \boldsymbol{b}$ are the velocity and magnetic field (in Alfvén velocity units), $P^\star$ is the total pressure, $\nu, \eta$ are the viscosity and magnetic diffusivity and $\boldsymbol{F}_V$ is a forcing term. We will use in this paper two types of forcing. The first one is the ABC (Arnold–Beltrami–Childress) forcing, which is a superposition of Beltrami vortices and thus fully helical; it is defined as

(2.3) \begin{equation} \boldsymbol{f}_{\mathrm {ABC}} = \left [ \cos y + \sin z \right ] {\hat{\boldsymbol{x}}} + \left [\sin x + \cos z \right ] {\hat{\boldsymbol{y}}} + \left [ \cos x + \sin y \right ] {\hat{\boldsymbol{z}}} . \end{equation}

The ABC flow is an eigenfunction of the curl, and it is an exact solution of the Euler equations; thus, for large enough viscosity, it is stable but turbulence develops as the Reynolds number $R_V$ increases. We also take the Taylor–Green (TG) forcing written as

(2.4) \begin{equation} \boldsymbol{f}_{\mathrm {TG}} = f_0^{tg} \left \{ \left [ \sin x \ \cos y \ \cos z \right ] {\hat{\boldsymbol{x}}} - \left [\cos x \ \sin y \ \cos z \right ] {\hat{\boldsymbol{y}}} + 0 \ {\hat{\boldsymbol{z}}} \right \}\!, \end{equation}

with $f_0^{tg}=3$ . This forcing is globally non-helical, but it can be viewed in fact as an assembly of helical patches of both signs and varied intensities.

We solved numerically our MHD system in a fully periodic box using a classic pseudo-spectral solver, involving a $2/3$ dealiasing technique with a parallel CPU-MPI code (CUBBY; Ponty et al. Reference Ponty, Mininni, Montgomery, Pinton, Politano and Pouquet2005). With these two forcings, we compute four simulations altogether, for up to hundreds of thousands of eddy turn-over times, and we record the same number of snapshots for the three-dimensional field components for $\textbf{v}$ and $\textbf{b}$ . Some of the data and a few of the statistical properties of these runs are given in table 1.

Table 1. Characteristics of the runs, with the linear resolution $n_p$ of the cubic grid, $\nu$ the viscosity, the Reynolds number $R_V= L_{int} V_{rms}/\nu$ and Taylor Reynolds number $R_\lambda =\lambda V_{rms}/\nu$ ; $R_N= \eta _v n_p/3$ is the so-called Kaneda criterion based on the resolution in terms of the Kolmogorov dissipation length $\eta _v$ . TG denotes Taylor–Green runs, and ABC denotes ABC runs (see text). For each run, $T_{max}$ is its total duration, with $\tau _{nl}$ in these units between 1 and 2. For the magnetic variables, we give the magnetic Reynolds number $R_M= L_{int} V_{rms}/\eta$ , the magnetic Prandtl number $P_m=\nu /\eta$ and $r=R_M/R_M^C$ is the ratio of the magnetic Reynolds number to the (approximate) critical value for the threshold of the dynamo. All these non-dimensional numbers need different definition of scales, like the integral scale $L_{int}=2 \pi \sum E(k)/k/\sum E(k)$ defined using the isotropic energy spectrum computed along the simulation at each wavenumber $k$ , the Taylor scale $\lambda =\sqrt {10}\eta _v R_V^{1/4}=\sqrt {10}L_{int}R_V^{-1/2}$ in the inertial range and the Kolmogorov scale $\eta _v= [\frac {\epsilon _v}{\nu ^3}]^{-1/4} = L_{int}R_V^{-3/4}$ at the onset of the dissipation range. We need also one characteristic velocity which is usually taken as the root mean square of the kinetic energy $V_{rms}= \sqrt {2 \sum E(k)}$ . The nonlinear time is taken as $\tau _{nl}= L_{int}/V_{rms}$ . All the scales and velocity are averaged in time as the simulations develop.

2.2. Brief description of the runs

The dynamo problem, concerning the generation of magnetic fields at both large and small scales, is a long-standing topic (see Brandenburg & Subramanian Reference Brandenburg and Subramanian2005; Rincon Reference Rincon2019 for thorough recent reviews). In the context of this paper, we analyze four simulation runs, focusing on the turbulent dynamo regime. In these simulations, the (fast) dynamo is triggered when the magnetic Reynolds number $R_M$ exceeds a threshold that depends on the magnetic Prandtl number $P_m=\nu /\eta$ (Ponty et al. Reference Ponty, Mininni, Montgomery, Pinton, Politano and Pouquet2005), as seen in run TG3. Sub-critical dynamos can also be observed, where magnetically induced changes to the velocity field play a role, such as in run TG1, which is close to the onset of dynamo action, and living on the sub-critical branch (Ponty et al. Reference Ponty, Laval, Dubrulle, Daviaud and Pinton2007). It is also worth noting that the Lorenz force’s feedback on the velocity field can influence the so-called ‘on–off’ intermittency of the dynamo, as explored in detail by Alexakis & Ponty (Reference Alexakis and Ponty2008) in the context of the ABC runs ABC1 and ABC2. These runs are notable for their short off phases, during which the magnetic field becomes weak enough to revert the system to a hydrodynamic state. The dynamo comes back quickly, with a return to an MHD equilibrium (see figure 1 below). We thus finally have two different types of dynamo turbulence, with a large amount of fluctuations which are analyzed in the next sections, examining the behavior of the third- and fourth-order normalized moment statistics using the full-space temporal field data.

3.1. Field gradient tensors, skewness and kurtosis

Concerning the data points for measurements, which must be statistically independent, they are taken approximately every $\tau _{nl}$ ; we recall that measurement errors go as $\sqrt {6/N_s}$ , with $N_s$ the number of independent data points (see e.g. Sura & Sardeshmukh Reference Sura and Sardeshmukh2008). Large samples are needed also because a parabolic fit is quite sensitive to extreme values. Note that it is shown in (Wan et al. Reference Wan, Oughton, Servidio and Matthaeus2010), in the context of two-dimensional (2-D) MHD turbulence, that an estimate of the kurtosis at small scales requires, in the framework of the DNSs analyzed in that paper, that the (Kolmogorov) dissipation scale be at least twice as large as the cutoff $k_{max}=n_p/3$ ; this condition is well fulfilled by the runs of table 1 (see parameter $R_N$ ).

We analyze the data using the point-wise rates of dissipation of the kinetic and magnetic energy, $ \epsilon _v(\boldsymbol{x}), \epsilon _m(\boldsymbol{x})$ . They are expressed in terms of the symmetric part of the velocity gradient tensor, namely the strain tensor $S^v_{ij}$

(3.1) \begin{equation} S^v_{ij}(\boldsymbol{x}) = \frac {\partial _j u_i (\boldsymbol{x})+ \partial _i u_j(\boldsymbol{x})}{2};\;\;\;\;\; \epsilon _v(\boldsymbol{x}) = {2 \nu } \Sigma _{ij} S^v_{ij}(\boldsymbol{x}) S^v_{ij}(\boldsymbol{x}) \, \end{equation}

and of the magnetic current density, viz. $\epsilon _m(\boldsymbol{x}) = {\eta } \boldsymbol{j}^2(\boldsymbol{x})$ . For completeness, we also define the symmetric part of the magnetic field gradient tensor, namely

(3.2) \begin{equation} S^m_{ij}(\boldsymbol{x}) = \frac {\partial _j b_i (\boldsymbol{x})+ \partial _i b_j(\boldsymbol{x})}{2};\;\;\;\;\; \sigma _m(\boldsymbol{x}) = \Sigma _{ij} S^m_{ij}(\boldsymbol{x}) S^m_{ij}(\boldsymbol{x}). \end{equation}

Note that $S^m_{ij}$ is a pseudo (axial) tensor. The standard expressions for the integrated (space-averaged) kinetic, magnetic and total energy dissipation rates, can respectively be written also as

(3.3) \begin{equation} \epsilon _V= \nu \big \lt |{\boldsymbol {\omega }}|^2 \big \gt;\;\;\;\;\; \epsilon _M = \ \eta \big \lt |\boldsymbol{j}^2| \big \gt \, \ \epsilon _T=\epsilon _V + \epsilon _M \, \end{equation}

with $\boldsymbol {\omega }=\nabla \times \boldsymbol{u}$ the vorticity. Finally, the skewness and excess kurtosis are written below for a scalar field $f$ , with $S_f=0, K_f=0$ for a Gaussian distribution, with the kurtosis (or flatness) being defined as $F_f =K_f+3$

(3.4) \begin{equation} S_f = \big\lt\, f^3 \big \gt / \big \lt\, f^2 \big \gt ^{3/2},\ K_f = \big \lt \, f^4 \big \gt / \big \lt \, f^2 \big \gt ^2-3. \end{equation}

3.2. Some models for $K(S)$ relations

A $K(S)$ parabolic law has been derived explicitly in Sura & Sardeshmukh (Reference Sura and Sardeshmukh2008) for a model of oceanic sea-surface temperature anomalies, based on the dynamics of a specific linear Langevin model with both additive and multiplicative noises. Analyzing the corresponding Fokker–Planck equation for the stationary PDF, these authors can show analytically that, in the limit of weak multiplicative noise, one has $K(S)= 3S^2/2$ . Multiple other studies show the plausibility of a Langevin model for parabolic $K(S)$ behaviors in different contexts as exemplified e.g. in Hasselmann (Reference Hasselmann1962) and Sattin et al. (Reference Sattin2009) (see also Pouquet et al. (Reference Pouquet, Marino, Politano, Ponty and Rosenberg2025) for nonlinear Langevin models). Note that, in a Langevin equation, in a sense, one is getting rid of the closure problem for turbulent motions since it is linear, with the complex nonlinear small-scale dynamics bundled up in a rapid stochastic forcing with an assumption of (mostly) local interactions among these fast motions.

In fact, models in the framework of fusion plasmas have also been written, for example in the context of magnetically confined experiments. In Garcia (Reference Garcia2012), the dynamics, as for dissipation events, is viewed as a random sequence of bursts as opposed to a quasi-continuum. These bursts are occurring independently and following a Poisson process. When taking for the shape of these bursts a sharp rise and a slow exponential decay, one can compute $K(S)$ relations which, for an exponential distribution of burst amplitudes, becomes $K(S)=3/2(S^2-1)$ .

Figure 1. Run ABC2: temporal evolutions, with kinetic variables in blue and magnetic ones in red. Top, left and right: kinetic energy and its dissipation, $\epsilon _v$ . Middle, left and right: magnetic energy and its dissipation, $\epsilon _m=\eta (\,\boldsymbol{j}^2)$ . Note the different units on axes, and the lin–log scale for magnetic energy. Bottom: probability density functions of the kinetic and magnetic energies (two left plots), and of their dissipation (two right plots), with lin–log plots used for the latter.

Figure 2. Similar plots as in figure 1 but now for the TG3 run.

3.3. Collecting the spatio-temporal data for further analysis

Our methodology is the following. For each field variable at a fixed observation time $T_O$ , as for the vertical component of the velocity $v_z(\boldsymbol{x}, T_O)$ , or the magnetic energy dissipation $\eta \boldsymbol{j}^2(\boldsymbol{x}, T_O)$ , the data are collected roughly every turn-over time $\tau _{nl}$ , for in excess of $5 \times 10^4$ outputs, as shown in column $T_{max}$ in table 1. For each full-cube temporal data output at $T_O$ the spatially dependent 3-D fields we analyze are computed point-wise at each $\boldsymbol{x}$ location in the $n_p^3$ data cube. The PDFs are constructed for these various fields, at that fixed time $T_O$ , and their second, third and fourth moments are evaluated to yield skewness and kurtosis for that time index $T_O$ . These data are then assembled in $K(S)$ plots for the different physical variables of interest. In figures 3 and 4, the time arrow of the data is indicated by the color of the points in the scatter plots, with a rainbow code as given by the color bar at the left of the plots, with purple at early times and red at late times. We note that, having in excess of $5 \times 10^4$ time stamps $T_O^{(i)}$ for all runs, the error on the skewness is less than $0.025$ , and twice that for the kurtosis. As expected, these systems achieve statistical equilibria rapidly, in a few $\tau _{nl}$ (see figures 1and 2 giving the whole time span of the runs for the kinetic and magnetic energies).

Note that this methodology does not correspond to a spatial analysis of intermittency at a fixed time, such as the first maximum of dissipation, given the moderate spatial resolution of the runs analyzed herein. This would lead to a study of the localization and intensity of extreme events at that time, but this is not the purpose of this study. Rather, we instead consider solely here the statistical variations over long times of the spatial intermittency of several field statistics in a global way. We do that in terms of overall variation of non-dimensionalized moments of various functionals of the velocity and magnetic fields as they evolve in time. With samples taken roughly every turn-over time, the statistical data points can be considered as independent. This leads to an investigation of the scatter in values of both $K$ and $S$ for a sample of field functionals, and thus of the $K(S)$ relationships that may emerge for them, with each point on the scatter plots of figures 3 and 4 corresponding to a given $T_O$ . We finally note that these analyses require substantial storage capabilities even at moderate spatial resolutions, so that the various statistical variables of interest and their moments are better computed a priori, and stored as the runs progress.

3.4. Overall data analysis

We give for the ABC2 flow (figure 1) and for the TG3 case (figure 2) the temporal evolution of the kinetic (top) and magnetic (middle) energy as a function of time on the left, and on the right of their respective total dissipations, $\epsilon _{V,M}$ . Note that time is expressed in output counts, with a turn-over time being roughly twice that. At the bottom are given the energy PDFs for the velocity (blue, leftmost) and the magnetic field (red, rightmost plots.) In all cases, there are sustained fluctuations in the amplitude of the fields, and in the case of the ABC run, lapses in both kinetic and magnetic energy corresponding to the on–off mechanism. This is directly related to the fact that the PDFs, in that case, have two relative maxima (with one peak close to zero for the magnetic field), whereas in the case of the TG forcing, the structure of the PDFs is closer to being unimodal. We note that exponential fits in the case of the dissipation fields are plausible, but not yet very strongly marked given the moderate range of Reynolds numbers in these runs.

Figure 3 gives $K(S)$ for various fields; at top, we display the $K(S)$ relationship for the z-component of the velocity, the square vorticity and the kinetic energy dissipation, namely $K_{v_z}(\boldsymbol{x}), K_{\omega ^2}(\boldsymbol{x})$ and $K_{\epsilon _v}(\boldsymbol{x})$ . At bottom we plot the equivalent fields for the magnetic induction, namely $K_{b_z}(\boldsymbol{x}), K_{\sigma _m}(\boldsymbol{x})$ and $K_{\eta \boldsymbol{j}^2}(\boldsymbol{x})$ . The blue lines correspond to $K=[3/2][S^2-1]$ as mentioned before (see § 3.2).

We note the following: whereas for NS turbulence, the three components of the velocity field are Gaussian with the $K(S)$ relation centered on the $(0,0)$ point, here, the peak in values for $K$ at $S\approx 0$ for $v_z$ is up to $K\approx 14$ , and rather narrowly centered around $S_{v_z}\approx 0$ ; high $K$ values are also present for $b_z$ . The hydrodynamic case analyzed in Pouquet et al. (Reference Pouquet, Rosenberg, Marino and Mininni2023) is computed at comparable $R_\lambda$ and $n_p$ (but not $T_{max}$ ), and both $S$ and $K$ for the velocity are close to $0$ . On the other hand, for stratified flows with or without rotation, the vertical component of the velocity, $v_z$ , can have high kurtosis and high values itself; this is associated with the intermittency of dissipation because of the variability of the system dominated by waves with the sudden development of small dissipative scales through shear-related instabilities. For this MHD run, the $x$ and $y$ components of the velocity behave approximately in the same way as $v_z$ (not shown), with a skewness comparable to that of $v_z$ but smaller kurtosis.

Figure 3. Top: for run TG3, $K(S)$ for the vertical velocity $v_z$ (left), the square vorticity density $\boldsymbol {\omega }^2$ (middle) and the point-wise kinetic energy dissipation $\epsilon _v$ (right). Bottom: $K(S)$ for $b_z$ (left), $\sigma _m$ (middle) and $\eta \boldsymbol{j}^2$ (right). The color bar at left indicates the temporal clock in units of turn-over times, with early (late) times in blue (red). The blue lines follow $K(S)=3/2[S^2 - 1]$ .

The symmetric and anti-symmetric parts of both the velocity and magnetic field gradient tensors have both high skewness and high kurtosis and for them, a quasi-parabolic fit is appropriate. The magnetic dissipation has higher kurtosis and thus is likely significantly more intermittent than its kinetic counterpart, and similarly for $\sigma _m$ and vorticity. On the other hand, both kinetic and magnetic dissipations have lower skewness and kurtosis than for the other part of their gradient tensors, although their statistics overall are similar. We recall that double exponential (Laplace), or Weibull distributions have small $S$ of either sign and high kurtosis (Bertin & Clusel (Reference Bertin and Clusel2006); Biri, Scharffenberg & Stammer Reference Biri, Scharffenberg and Stammer2015; Reference Aschwanden2016). In that context, we note that it is shown in Sorriso-Valvo et al. (Reference Sorriso-Valvo, Carbone, Perri, Greco, Marino and Bruno2018) that a proxy of energy transfer for the solar wind can be defined based on exact laws for MHD corresponding to the conservation of energy and cross-helicity $H_C=\left \lt \boldsymbol{v} \cdot \boldsymbol{b} \right \gt$ ; these proxy fields display high intermittency in Helios 2 (and Ulysses) data, with plausible stretched exponential fits. A final remark is that data points with $[K,S]\approx 0$ must be dominated by random noise at these times; they could correspond to relaxation intervals following sharp bursts in energy dissipation.

Figure 4. The first two rows are the same as figure 3 for run ABC2 with $R_V\approx 175$ and $T_{max}\approx 3.14 \times 10^5$ . Bottom: $K(S)$ is plotted for the magnetic energy $E_m$ (left) and for the symmetric part of the magnetic field gradient tensor $\sigma _m$ (right), using now log–log coordinates. Approximate fits give for both $\alpha \approx 2.2, \kappa \approx 1.0$ .

We now check whether this behavior is observed as well for another type of forcing. In figure 4 are plotted the $K(S)$ relationships at the top for $v_z$ (left), $\boldsymbol {\omega }^2$ (middle) and $\epsilon _v$ (right) and (middle row), $b_z$ (left), $\sigma _m$ (middle) and $\eta \boldsymbol{j}^2$ (right), as in figure 3 but now for run ABC2 with a fully helical ABC forcing. Values for $(S,K)$ for both runs are comparable except for the vertical component of the velocity due to its specific structure. We added in the last row of figure 4 the $K(S)$ plots for the magnetic energy (left) and $\sigma _m$ (right) for the same run, and using now log–log coordinates. Power laws are conspicuous above some threshold in skewness, with $E_m \approx 0.96S^{2.20}, \ \sigma _m \approx 0.96S^{2.25},\ \eta j^2\approx 1.05 S^{2.17}$ .

Finally, due to the strong symmetries of the initial conditions and forcing of the dynamos analyzed in this paper, one could wonder whether the addition of a small noise would change the results. On the other hand, in view of the length of the computations, well beyond what can be estimated for a reasonable Lyapounov time of separation of trajectories, it is unlikely that the overall results, and in particular the quasi-parabolic law for dissipation, would be altered. Indeed, one can find estimates of the first Lyapounov exponent $\lambda _1$ , in the ABC dynamo for example, with $\lambda _1\sim 0.073$ for a run with $R_V\approx 60, P_M=4$ , comparable to what we have here (Zienicke, Politano & Pouquet Reference Zienicke, Politano and Pouquet1998; Alexakis & Ponty Reference Alexakis and Ponty2008), meaning that, after roughly $14 \tau _{nl}$ , the initial conditions have been forgotten.