© The Authors 2024

1. Introduction

Decaying turbulence played an important role in the early Universe during the radiation-dominated era, when the magnetic field is well coupled to the plasma. While turbulence speeds up the decay, it can also lead to a significant increase in the typical length scale, which could then be many times larger than the comoving horizon scale at the time of magnetic field generation (Brandenburg et al. 1996; Christensson et al. 2001; Banerjee & Jedamzik 2004). This is important because magnetogenesis processes during the electroweak era, when the age of the Universe was just a few picoseconds (Vachaspati 1991; Cheng & Olinto 1994; Baym et al. 1996), tend to produce magnetic fields of very small length scales of the order of 1 AU or less (in comoving units).

From the study of decaying hydrodynamic turbulence, it has been known for a long time that turbulent energy density and length scale evolve as power laws (Batchelor 1953; Saffman 1967). The exponents depend on the physics of the decay, specifically on the possibility of a conserved quantity that governs the decay, for example magnetic helicity in the hydromagnetic case (Hatori 1984; Biskamp & Müller 1999). The endpoints of the evolution, however, depend on the relevant timescale, which is traditionally taken to be just the turnover or, in the magnetic case that we consider here, the Alfvén time (Banerjee & Jedamzik 2004). More recently, Hosking & Schekochihin (2023) argued that the turnover time should be replaced by the reconnection time, which could be significantly longer (up to 10^5.5 times). This would result in an endpoint where the magnetic field strength is greater and the turbulent length scale smaller than otherwise, when the decay time is just the Alfvén time.

One of the hallmarks of turbulence is that its large-scale properties are nearly independent of viscosity and resistivity, which act predominantly on the smallest scales of the system. On the other hand, magnetic reconnection is a process that could potentially limit the speed of the inverse cascade. This idea has been invoked by Hosking & Schekochihin (2023) to explain a premature termination of the decay process by the time of recombination of the Universe, when its age was about 400 000 years.

The purpose of our paper is to analyze numerical simulations with respect to their decay times at different values of resistivity. In Sect. 2 we discuss the decay time and its relation to other quantities. In Sect. 3 we present our numerical simulation setup and show the results for a resistivity-dependent decay in Sect. 4. We then make a comparison with a purely hydrodynamic decay in Sect. 5 and with the two-dimensional (2D) hydromagnetic case in Sect. 6, before concluding in Sect. 8. In Appendix A we provide a historical note on the anastrophy, i.e., the mean squared magnetic vector potential, and in Appendix B, we show detailed convergence tests for some of our 2D results.

2. Decay and turnover times

In the following we focus on the decay of magnetic field. Magnetically dominated turbulence is characterized by the turbulent magnetic energy density ℰ_M and the magnetic integral scale ξ_M. Both ℰ_M(t) and ξ_M(t) can be defined in terms of the magnetic energy spectrum E_M(k, t), such that ℰ_M = ∫E_M dk and ξ_M = ∫k⁻¹E_M dk/ℰ_M. In decaying turbulence, both quantities depend algebraically rather than exponentially on time. Therefore, the decay is primarily characterized by power laws,

$$\begin{array}{c}\hfill {\mathcal{E}}_{\mathrm{M}}\propto t^{-p}\mathrm{and}{\xi }_{\mathrm{M}}\propto t^q,\end{array}$$(1)

rather than by exponential laws of the type ℰ_M ∝ e^−t/τ. The algebraic decay is mainly a consequence of nonlinearity. On the other hand, in decaying hydromagnetic turbulence with significant cross-helicity, for example, the nonlinearity in the induction equation is reduced and then the decay is indeed no longer algebraic, but closer to exponential (see Brandenburg & Oughton 2018).

An obvious difference between algebraic and exponential decays is that in the former ℰ_M(t) is characterized by the nondimensional quantity p, while in the latter it is characterized by the dimensionful quantity τ. Following Hosking & Schekochihin (2023), a decay time τ can also be defined for an algebraic decay and is then given by

$$\begin{array}{c}\hfill {\tau }^{-1}=-\mathrm{d}ln{\mathcal{E}}_{\mathrm{M}}/\mathrm{d}t.\end{array}$$(2)

In the present case of a power-law decay, this value of τ = τ(t) is time-dependent and can be related to the instantaneous decay exponent

$$\begin{array}{c}\hfill p(t)=-\mathrm{d}ln{\mathcal{E}}_{\mathrm{M}}/\mathrm{d}lnt\end{array}$$(3)

through τ = t/p(t) (i.e., no new parameter emerges except for t itself). However, a useful way of incorporating new information is by relating τ to the Alfvén time τ_A = ξ_M/v_A through

$$\begin{array}{c}\hfill \tau =C_{\mathrm{M}}{\xi }_{\mathrm{M}}/v_{\mathrm{A}},\end{array}$$(4)

where C_M is a nondimensional parameter, and v_A is the Alfvén velocity, which is related to the magnetic energy density through ${\mathcal{E}}_{\text{M}}=B_{\text{rms}}^2/2{\mu }_0=\rho v_{\text{A}}^2/2$, where ρ is the density, μ₀ the vacuum permeability, and B_rms the root mean square (rms) magnetic field.

As was noted by Hosking & Schekochihin (2023), Eq. (4) can be used to define the endpoints of the evolutionary tracks in a diagram of B_rms versus ξ_M or, equivalently, v_A versus ξ_M (i.e., v_A = v_A(ξ)). They also noted that the location of these endpoints is sensitive to whether or not C_M depends on the resistivity of the plasma. If it does depend on the resistivity, this could be ascribed to the effects of magnetic reconnection, which might slow down the turbulent decay.

Magnetic reconnection refers to a change in magnetic field line connectivity that is subject to topological constraints. A standard example is x-point reconnection (see, e.g., Craig & McClymont 1991; Craig et al. 2005), which becomes slower as the x-point gets degenerated into an extremely elongated structure (Sweet 1958; Parker 1957; see Liu et al. 2022 for a review). It is usually believed that in the presence of turbulence, such structures break up into progressively smaller ones, which makes reconnection eventually fast (i.e., independent of the microphysical resistivity; Galsgaard & Nordlund 1996; Lazarian & Vishniac 1999; Comisso & Sironi 2019). However, whether this would also imply that τ becomes independent of the resistivity remains an unclear issue. For example, Galsgaard & Nordlund (1996) found that resistive heating becomes independent of the value of the resistivity. Another question concerns the speed at which magnetic flux can be processed through a current sheet (Kowal et al. 2009; Loureiro et al. 2012). Also of interest is the timescale on which the topology of the magnetic field changes (Lazarian et al. 2020). These different timescales may not all address the value of C_M that relates the decay time to the Alfvén time.

In magnetically dominated turbulence, the effect of the resistivity is quantified by the Lundquist number. For decaying turbulence, it is time-dependent and defined as

$$\begin{array}{c}\hfill \mathrm{Lu}(t)=v_{\mathrm{A}}(t){\xi }_{\mathrm{M}}(t)/\eta .\end{array}$$(5)

This quantity is similar to the magnetic Reynolds number if we replace v_A by the rms velocity, u_rms. Here, however, the plasma is driven by the Lorentz force, so the Lundquist number is a more direct way of quantifying the resistivity than the magnetic Reynolds number. The Alfvénic Mach number is defined as Ma_A = u_rms/v_A.

In addition to varying η, we also vary ν such that the magnetic Prandtl number Pr_M = ν/η is typically in the range 1 ≤ Pr_M ≤ 5. It is then also convenient to define the Lundquist number based on the reconnection outflow,

$$\begin{array}{c}\hfill {\mathrm{Lu}}_{\nu }(t)=\mathrm{Lu}(t)/\sqrt{1+{\mathrm{Pr}}_{\mathrm{M}}}\end{array}$$(6)

(see Hosking & Schekochihin 2023 for details).

3. Numerical simulations

We performed simulations of the compressible hydromagnetic equations for the magnetic vector potential A, the velocity u, and the logarithmic density lnρ in the presence of viscosity ν and magnetic diffusivity η:

$$\begin{array}{cc}\hfill \frac{\partial \mathit{A}}{\partial t}& =\mathit{u}\times \mathit{B}+\eta {\mathrm{\nabla }}^2\mathit{A},\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill \frac{\mathrm{D}\mathit{u}}{\mathrm{D}t}& =\frac{1}{\rho }[\mathit{J}\times \mathit{B}+\mathbf{\nabla }·(2\rho \nu \mathsf{S})]-c_{\mathrm{s}}^2\mathbf{\nabla }ln\rho ,\hfill \end{array}$$(8)

$$\begin{array}{cc}\hfill \frac{\mathrm{D}ln\rho }{\mathrm{D}t}& =-\mathbf{\nabla }·\mathit{u}.\hfill \end{array}$$(9)

We used a random magnetic field as the initial condition, such that E_M(k, 0) has a k⁴ subinertial range for k < k_p (Durrer & Caprini 2003), and a k⁻² inertial range for k > k_p (Brandenburg et al. 2015). In all cases we chose k_p/k₁ = 60, where k₁ = 2π/L is the smallest wavenumber in our cubic domain of size L³. In Eqs. (7) and (8), B = ∇ × A is the magnetic field, J = ∇ × B/μ₀ is the current density, and S_ij = (∂_iu_j + ∂_ju_i)/2 − δ_ij∇ ⋅ u/3 are the components of the rate-of-strain tensor S.

There is no magnetic helicity on average, but the fluctuations in the local magnetic helicity density h = A ⋅ B lead to a decay behavior where the correlation integral of h, which is also known as the Hosking integral, is conserved (Hosking & Schekochihin 2021; Schekochihin 2022; Zhou et al. 2022). We use the PENCIL CODE (Pencil Code Collaboration 2021), which has also been used for many earlier simulations of decaying hydromagnetic turbulence (Zhou et al. 2022; Brandenburg et al. 2023). All our simulations are in the magnetically dominated regime, because the velocity is just a consequence of and driven by the magnetic field.

Because the magnetic field is initially random, the resulting velocity is also random and it drives a forward turbulent cascade with kinetic and magnetic energy dissipation rates ϵ_K = ⟨2νρS²⟩ and ϵ_M = ⟨ημ₀J²⟩. Their ratio scales with ${\mathrm{Pr}}_{\mathrm{M}}^{1/3}$ (Brandenburg 2014; Galishnikova et al. 2022). If k_p/k₁ is large (we recall that we use the value 60), there is also an inverse cascade (Brandenburg et al. 2015). The inverse cascade was also found in the relativistic regime (Zrake 2014) and is now understood to be a consequence of the conservation of the Hosking integral (Hosking & Schekochihin 2021; Schekochihin 2022; Zhou et al. 2022). However, the role played by the Hosking integral is currently not universally accepted (Armua et al. 2023; Dwivedi et al. 2024). The lack of numerical support could be related to insufficiently large values of k_p/k₁.

We define the kinetic energy spectrum E_K(k, t) analogously to E_M(k, t), using the normalization $\int E_{\text{K}}(k,t)\text{d}k={\mathcal{E}}_{\text{K}}(t)\equiv {\rho }_0{u}_{\text{rms}}^2/2$, where ρ₀ = ⟨ρ⟩=const owing to mass conservation. In magnetically driven turbulence, ℰ_K is about one-tenth of ${\mathcal{E}}_{\text{M}}=B_{\text{rms}}^2/2{\mu }_0$, which seems to be surprisingly independent of the physical input parameters (Brandenburg et al. 2017).

In Table 1 we summarize the results of five simulations, Runs M0–M5, where we vary η and ν and also vary the number of mesh points, N³. We present some relevant output parameters that are defined below. They are all obtained from a statistically steady stretch of our time series data, and the error bars were calculated as the largest departure from any one-third of the time series. A visualization of the z components of B on the periphery of the computational domain for Run M1 is shown in Fig. 1. Figure 2 shows J_z in an xy plane together with a zoom-in on the lower left corner of the domain. The corresponding magnetic and kinetic energy spectra are shown in Fig. 3 for different times in units of [t]=(c_sk₁)⁻¹.

Fig. 1. Visualization of Bz on the periphery of the computational domain for Run M1, where PrM = 5.

Fig. 2. Visualization of Jz at z = π along with a zoom-in on the lower left corner for Run M0, where PrM = 10, at t = 644. The value of ξM at that time is indicated by the length of the short white lines.

Fig. 3. Magnetic and kinetic energy spectra for Run M1 at different times. For large values of k there is a k−2 inertial range. The ∝k4 and ∝k2 slopes are indicated for reference. As elsewhere, time is in units of [t]=(csk1)−1.

Fig. 4. Approach of CM(t), Lu(t), and $C_{\mathrm{L}}^{(1/4)}(t)$ to an approximately constant value toward the end of the simulation for Run M1.

The E_M(k, t)∝k⁻² inertial range can be explained by weak turbulence scaling (Brandenburg et al. 2015), but it becomes a bit shallower near the dissipation subrange. This could follow from some kind of magnetic bottleneck effect associated with reconnection, analogously to the bottleneck effect in hydrodynamic turbulence (Falkovich 1994).

It should be noted that at late times, the subinertial range spectrum of E_M(k, t) becomes shallower than the initial k⁴ slope. This is an artifact of poor scale separation (see Brandenburg et al. 2023 for related numerical evidence). On the other hand, a k² subinertial range for E_K(k, t) has been seen for some time (see Kahniashvili et al. 2013).

In our numerical simulations we use units such that c_s = k₁ = ρ₀ = μ₀ = 1. The resistivity, μ₀η, is therefore the same as the magnetic diffusivity. Nevertheless, most of the results below are expressed in manifestly nondimensional form.

4. Resistivity-dependent decay

We now analyze a collection of runs similar to those of the recent works of Zhou et al. (2022) and Brandenburg et al. (2023), who considered different values of Lu and also included some runs with hyperviscosity and hyperresistivity, unlike what we present in the present work. At variance with those earlier papers, where the focus was always on the decay exponent p(t), here we focus on the evolution of the decay time, τ(t) = t/p(t).

4.1. Decay time

The goal is to determine the prefactor in the scaling relation τ ∝ ξ_M/v_A. Therefore, we write

$$\begin{array}{c}\hfill t/p\equiv \tau =C_{\mathrm{M}}{\xi }_{\mathrm{M}}/v_{\mathrm{A}}\end{array}$$(10)

and determine

$$\begin{array}{c}\hfill C_{\mathrm{M}}(t)=(t/p)v_{\mathrm{A}}(t)/{\xi }_{\mathrm{M}}(t).\end{array}$$(11)

We emphasize here that each of the terms is time-dependent, including τ(t) = t/p(t), as already noted above. Interestingly, it turns out that the quantity C_M(t) eventually settles around a plateau:

$$\begin{array}{c}\hfill C_{\mathrm{M}}=\underset{t\to \infty }{lim}{C}_{\mathrm{M}}(t).\end{array}$$(12)

Here and in the following, when time-dependence is not indicated, we usually mean that the value is obtained as a suitable limit of the corresponding time-dependent function. In numerical simulations with finite domains, the limit t → ∞ needs to be evaluated with some care to prevent the final result from being contaminated by finite size effects. We do this by selecting a suitable time interval during which certain data combinations are approximately statistically stationary. We refer to these results as late-time limits.

In Fig. 4 we show that C_M(t) and Lu(t) approach an approximately constant value toward the end of the run. Time is given both in units of [t]=(c_sk₁)⁻¹ and in initial Alfvén times, (v_A0k_p)⁻¹, where $v_{\mathrm{A}0}=\sqrt{2{\mathcal{E}}_{\mathrm{M}0}}$ is the initial Alfvén velocity and ℰ_M0 the initial magnetic energy density. Allowing for the possibility of power-law scaling, $C_{\mathrm{M}}=C_{\mathrm{L}}^{(n)}{\mathrm{Lu}}^n$, we also plot the prefactor $C_{\mathrm{L}}^{(n)}(t)$ for n = 1/4 in the last panel of Fig. 4 (see details below). Toward the end of the simulation, however, there is an increase in the fluctuations, which follows from the decrease in the magnetic energy.

As discussed below in more detail in connection with 2D simulations, current sheets are underresolved at early times, when ξ_M is small. As we now see from Fig. 4, this underestimates the resulting value of C_M. However, for t > 100, C_M(t) approaches a plateau, suggesting that the simulation now begins to be sufficiently well resolved, at least for the purpose of determining C_M.

In Fig. 5 we show the dependence of C_M on Lu and Lu_ν. We see an approximate scaling ∝Luⁿ with n = 1/4 for Lu < 6000 and a piecewise power-law scaling $\propto {\text{Lu}}_{\nu }^{1/2}$, but with different prefactors for the runs with larger and smaller values of the viscosity. We note that for Pr_M ≫ 1 we have ${\mathrm{Lu}}_{\nu }=v_{\mathrm{A}}{\xi }_{\mathrm{M}}/\sqrt{\eta \nu }$, which explains why the ${\text{Lu}}_{\nu }^{1/2}$ scaling is found to be compatible with being ∝Lu^1/4. We also note, however, that for a higher viscosity, the line ${\text{Lu}}_{\nu }^{1/2}$ is shifted upward (toward larger values of C_M). Owing to the more complicated combined dependence on ν and Lu_ν, we continue to employ the simpler Luⁿ scaling for the following discussion. It is worth noting that the exponent n = 1/4 is discussed in Uzdensky & Loureiro (2016) in connection with the fast growing mode of the tearing instability.

Fig. 5. Dependence of CM on Lu and Luν. An approximate scaling ∝Lu1/4 is found for Lu < 6000 and a piecewise power-law scaling $C_{\text{M}}=C_{\text{L}}^{(n)}{\text{Lu}}^n$ with different prefactors for the runs with higher viscosity (5 × 10−6 red symbols and 4 × 10−6 for the green symbol) and lower viscosity (2 × 10−6 blue symbols).

If it is really true that C_M is proportional to Luⁿ, as found above, we can write $C_{\mathrm{M}}=C_{\mathrm{L}}^{(n)}{\mathrm{Lu}}^n$, and then determine $C_{\mathrm{L}}^{(n)}$ as the late-time limit of

$$\begin{array}{c}\hfill C_{\mathrm{L}}^{(n)}(t)=\frac{t}{p}\frac{v_{\mathrm{A}}^{1-n}}{{{\xi }_{\mathrm{M}}}^{1+n}}{\eta }^n.\end{array}$$(13)

This is the formula that was used to compute $C_{\mathrm{L}}^{(n)}$ and to plot its time-dependence $C_{\mathrm{L}}^{(n)}(t)$ for n = 1/4 (see the last panel of Fig. 4). From Table 1 we find $3.3\le C_{\mathrm{L}}^{(1/4)}\le 3.8$.

4.2. Evolution of ℰ_M(t) and ξ_M(t)

It is of interest to know whether the resistivity dependence of C_M(t) is equally distributed among ξ_M and v_A (or ℰ_M). Figure 6 shows that there are indeed systematic differences in the decay laws for different values of the resistivity, but the differences are small and easily overlooked.

Fig. 6. Dependence of ℰM(t) and ξM(t) for three values of PrM.

To examine the dependence of C_M(t) on Lu, we write the decay laws for ξ_M(t) and ℰ_M(t) in a more detailed form than Eq. (1):

$$\begin{array}{cc}\hfill {\xi }_{\mathrm{M}}(t)& ={\xi }_{\mathrm{M}0}{(1+t/{\tau }_{\xi })}^q,\hfill \end{array}$$(14)

$$\begin{array}{cc}\hfill {\mathcal{E}}_{\mathrm{M}}(t)& ={\mathcal{E}}_{\mathrm{M}0}{(1+t/{\tau }_{\mathcal{E}})}^{-p}.\hfill \end{array}$$(15)

Here the coefficients ξ_M0 and ℰ_M0 just depend on the initial condition, and are thus not dependent on Lu, which means that the Lu-dependence can only enter through the coefficients τ_ξ and τ_ℰ. We can determine them as the limits of the time-dependent functions τ_ξ(t) and τ_ℰ(t), which are obtained by inverting Eqs. (14) and (15), and are given by

$$\begin{array}{c}\hfill {\tau }_{\xi }(t)=\frac{t}{{[{\xi }_{\mathrm{M}}(t)/{\xi }_{\mathrm{M}0}]}^{1/q}-1},\end{array}$$(16)

$$\begin{array}{c}\hfill {\tau }_{\mathcal{E}}(t)=\frac{t}{{[{\mathcal{E}}_{\mathrm{M}}(t)/{\mathcal{E}}_{\mathrm{M}0}]}^{-1/p}-1}.\end{array}$$(17)

Figure 7 shows the evolution of τ_ξ(t) and τ_ℰ(t) for Runs M3–M5 with fixed and moderate viscosity.

Fig. 7. Evolution of τℰ(t) and τξ(t) for three values of PrM with fixed and moderate viscosity. These times tend to be approximately constant at late times with values approximately consistent with those in Table 1.

Using the late-time limits of Eqs. (14) and (15), and those of Eqs. (16) and (17), the equation for C_M can now be decomposed in the form

$$\begin{array}{c}\hfill C_{\mathrm{M}}=C_{\xi }{C}_{\mathcal{E}},\end{array}$$(18)

such that Eq. (11) is obeyed. Here we can determine C_ξ and C_ℰ as the late-time limits of

$$\begin{array}{c}\hfill C_{\xi }(t)={\tau }^q{{\xi }_{\mathrm{M}}}^{-1}\mathrm{and}{C}_{\mathcal{E}}(t)={\tau }^{p/2}{v}_{\mathrm{A}}.\end{array}$$(19)

In this connection, it it important to remember that for a self-similar evolution (Brandenburg & Kahniashvili 2017), which is here approximately satisfied (see Fig. 3), we have

$$\begin{array}{c}\hfill q+p/2=1,\end{array}$$(20)

so that C_M = (t/p) v_A/ξ_M is obeyed. For t ≫ τ_ℰ, τ_ξ, using again τ = τ(t) = t/p(t), we have

$$\begin{array}{c}\hfill C_{\xi }\approx {({\tau }_{\xi }/p)}^q{{\xi }_{\mathrm{M}0}}^{-1}\mathrm{and}{C}_{\mathcal{E}}\approx {({\tau }_{\mathcal{E}}/p)}^{p/2}{v}_{\mathrm{A}0},\end{array}$$(21)

so that $C_{\mathrm{M}}=(v_{\mathrm{A}0}/p{\xi }_{\mathrm{M}0}){\tau }_{\xi }^q{\tau }_{\mathcal{E}}^{p/2}$. It is then natural to expect that both τ_ξ and τ_ℰ scale in the same way with Lu as C_M itself. Equation (21) also allows us to compute the timescales τ_ξ = p (C_ξξ_M0)^1/q and τ_ℰ = p (C_ℰ/v_A0)^2/p.

In Fig. 8 we show the dependence of τ_ξ and τ_ℰ on Lu. While the two show a similar dependence, approximately ∝Lu^1/4, we note that there is also possible evidence for a leveling off for larger values of Lu.

Fig. 8. Dependence of τξ and τℰ on Lu. The two show a similar dependence on Lu, approximately ∝Lu1/4, but there is possible evidence for a leveling off for larger values of Lu.

5. Comparison with hydrodynamic decay

Hydrodynamic decay is characterized by the kinetic energy density ${\mathcal{E}}_{\text{K}}(t)={\rho }_0{u}_{\text{rms}}^2/2$ and the hydrodynamic integral scale ξ_K(t) = ∫k⁻¹E_K dk/ℰ_K. We define the instantaneous kinetic energy decay exponent p_K(t) = − dlnℰ_K/dlnt and the decay time τ_K(t) = t/p_K(t). We relate τ_K(t) to the turnover time u_rms/ξ_K through τ_K = C_Kξ_K/u_rms, and thus determine C_K as the late-time limit of C_K(t) = [t/p_K(t)] u_rms(t)/ξ_K(t), which is defined analogously to Eq. (11).

Purely hydrodynamic simulations can be performed by just ignoring the magnetic field, or putting B = 0. In Table 2 we summarize such simulations for different values of ν, which is quantified by the Reynolds number: Re = u_rmsξ_K/ν. Figure 9 shows, for Re ≥ 100, that C_K does not change much with Re. This was not expected. It also confirms that the prolonged decay time found in Sect. 4 is indeed a purely magnetic phenomenon. Whether or not the resistivity dependence in the magnetic case must be ascribed to reconnection remains an open question. As discussed in Sect. 2, magnetic reconnection refers to topologically constrained changes of magnetic field lines, but in the present case of a turbulent magnetic field there is no connection with the standard picture of reconnection. An obvious alternative candidate for explaining the resistively prolonged decay time may be related to magnetic helicity conservation in local patches, as described by the conservation of the Hosking integral. While this idea seems more plausible to us, it is not obvious how to distinguish reconnection from magnetic helicity conservation in patches. It is true that magnetic helicity would vanish in two dimensions, but in that case it should be replaced by the anastrophy (see Appendix A for a historical note on this word). One aspect that might be different between the concepts of reconnection and magnetic helicity conservation in patches could be the dependence on Pr_M. Our present results have not yet shown such a dependence, which might support the idea that the dependence on resistivity is related to magnetic helicity conservation in patches.

Fig. 9. Dependence of CK on Re. The line Re0.1 is shown for comparison, but the data are also nearly compatible with being independent of Re.

6. Hydromagnetic decay in two dimensions

We now perform 2D simulations where $\mathit{B}=\mathbf{\nabla }\times (\widehat{\mathit{z}}{A}_z)$ lies entirely in the xy plane (see Table 3 for a summary). Equation (7) then reduces to

$$\begin{array}{c}\hfill \frac{\mathrm{D}{A}_z}{\mathrm{D}t}=\eta {\mathrm{\nabla }}^2{A}_z,\end{array}$$(22)

which obeys conservation of anastrophy, $\langle A_z^2\rangle =\text{const}$ (Fyfe & Montgomery 1976; Pouquet 1978, 1993), and the magnetic helicity density vanishes pointwise, so the Hosking integral is then also zero. These simulations are different from the recent ones by Dwivedi et al. (2024), who performed 2.5D simulations. In their case, there was a magnetic field component out of the plane. The anastrophy was then not conserved and the Hosking integral was finite.

As we stated above, our present measurements of C_M as a function of Lu cannot directly be compared with the reconnection rate determined by Loureiro et al. (2012), Comisso et al. (2015), or Comisso & Bhattacharjee (2016). In addition, it is not obvious how relevant a 2D simulation is in the present context because in 2D the anastrophy is conserved, while the Hosking integral is strictly vanishing.

In our present purely 2D simulations, values of Lu up to about 3 × 10⁵ have been reached. To allow for a longer nearly self-similar evolution, we used k_p/k₁ = 200 instead of 60. The simulation results give the prefactors in the scaling expected from anastrophy conservation as

$$\begin{array}{c}\hfill {\xi }_{\mathrm{M}}(t)\approx 0.13{⟨A_z^2⟩}^{1/4}{t}^{1/2},{\mathcal{E}}_{\mathrm{M}}(t)\approx 15{⟨A_z^2⟩}^{1/2}{t}^{-1},\end{array}$$(23)

where μ₀ = ρ₀ = 1 is used. We also see from Fig. 10 that the spectral peaks evolve underneath an envelope $E_{\text{M}}(k,t)\le 60\langle A_z^2\rangle k$.

Fig. 10. Similar to Fig. 3, but for Run 2M1. The magnetic peaks lie underneath a kβ envelope with β = 1, as expected in the case of anastrophy conservation.

Figure 11 shows magnetic and kinetic energy spectra for Run 2m2 with Lu = 1.8 × 10⁵ collapsed on top of each other by plotting ${\xi }_{\text{M}}{}^{\beta }(t)E_{\text{M}}(k{\xi }_{\text{M}}(t))$ versus kξ_M(t) for β = 1. This plot suggests that the subinertial range scalings of E_M(k, t) and E_K(k, t) are proportional to k⁵ and k³, respectively. Thus, they are steeper than expected in 3D. This behavior is in some ways similar to the steepening observed for helical decaying magnetic fields (see Brandenburg & Kahniashvili 2017).

Fig. 11. Magnetic energy spectra (left) and kinetic energy spectra (right) for Run 2m2 with Lu = 1.8 × 105 and PrM = 100, at times t = 5, 20, 100, and 400, collapsed on top of each other by plotting both vs kξM(t) and scaling them with ξM. This makes their heights agree, as expected.

In Table 3 we also give the ratio ϵ_K/ϵ_M, which is seen to increase with Pr_M (see Fig. 12). This was expected based on earlier results (Brandenburg 2014; Galishnikova et al. 2022), but it was never shown in the decaying case of magnetically dominated turbulence. We see that ${ϵ}_{\mathrm{K}}/{ϵ}_{\mathrm{M}}\propto {\mathrm{Pr}}_{\mathrm{M}}^{0.7}$, which is similar to the previously studied case with large-scale dynamo action rather than the case with just small-scale action where the slope was shallower. The ratio ϵ_K/ϵ_M is also found to increase with Lu, at least for Lu ≪ 10⁶. As discussed in Brandenburg & Rempel (2019), an increase in ϵ_K/ϵ_M with Pr_M may have implications for heating the solar corona, where the possible dominance of kinetic energy dissipation over Joule dissipation is not generally appreciated (see Rappazzo et al. 2007, 2018 for earlier work discussing this ratio).

Fig. 12. Dependence of ϵK/ϵM on PrM. The solid line denotes $1.3{\mathrm{Pr}}_{\mathrm{M}}^{0.7}$, but many of the data points, predominantly those with smaller Lu, are beneath that line. The inset shows that for PrM = 100, ϵK/ϵM increases with Lu like Lu0.6.

Although 2D and 3D runs are in many ways rather different from each other, we now determine the same diagnostics as in the 3D case (see Table 3 and Fig. 13 for a plot of C_M versus Lu). We see that the C_M dependence on Lu is qualitatively similar for 2D and 3D turbulence. Moreover, it becomes shallower for larger values of Lu. There is now evidence that C_M levels off and becomes independent of Lu. It is possible to fit our data to a function of the form

$$\begin{array}{c}\hfill C_{\mathrm{M}}(\mathrm{Lu})\approx C_{\mathrm{L}}^{(1/4)}{\left(\frac{\mathrm{Lu}}{1+\mathrm{Lu}/{\mathrm{Lu}}_{\mathrm{c}}}\right)}^n,\end{array}$$(24)

Fig. 13. Dependence of CM on Lu for the 2D runs. Shallow scaling ∝Lu0.1 is found for 104 < Lu < 105, which is also compatible with a leveling off at Luc = 2.5 × 104, as described by Eq. (24) with n = 1/4. The black (red) data points are for PrM = 5 (PrM = 1). The blue data points denote the 3D results from Sect. 4. The orange symbols are for the runs with PrM = 10 and 20, listed in Table 3. The thin dotted line gives Eq. (24) with n = 1/4 for comparison.

where Lu_c = 2.5 × 10⁴ is a critical Lundquist number characterizing the point where the dependence levels off, n = 1/4, and $C_{\mathrm{L}}^{(1/4)}=3.7$. The value of Lu_c is larger than that found by Loureiro et al. (2012), where the asymptotes for small and large Lundquist numbers cross at a value closer to 5000. However, this difference could simply be related to different definitions of the relevant length scales. We note that our definition of ξ_M does not include a 2π factor. A comparison with the Sweet–Parker value of n = 1/2 results in reasonable agreement for small values of Lu, but there are rather noticeable departures from the data for intermediate values.

Our 2D and 3D results for C_M are seen to be in good agreement with each other, which is similar to the observation made by Bhat et al. (2021). Obviously, larger simulations should still be performed to see whether the agreement between 2D and 3D continues to larger Lundquist numbers.

It should also be noted that some of our data with their nominal error bars do not lie on the fit given by Eq. (24) with n = 1/4. Especially for Pr_M = 100 and smaller values of Lu, we see that C_M lies systematically below the fit. It should be noted, however, that we would expect an increase with Pr_M as (1 + Pr_M)^1/2, which is the opposite trend, if the reconnection phenomenology were applicable (see Eq. (16) in Hosking & Schekochihin 2023).

To check whether our simulations are sufficiently well resolved, we show in Fig. 14 a visualization of J_z(x, y) for Run 2m6 with Pr_M = 10 and 16384² mesh points at t = 464 for a small part of the domain with sizes 2.8ξ_M(t)×0.74ξ_M(t) where a large current sheet breaks up into smaller plasmoids. A comparison between Runs 2m5 and 2m6 with 8192² and 16384² mesh points is shown in Appendix B. Our main conclusions are that higher resolution suppresses the tendency to produce ringing (i.e., the formation of oscillations on the grid scale), but that the results for C_M are not very sensitive to the numerical resolutions, as can be seen by comparing Runs 2m5 and 2m6 in Table 3.

Fig. 14. Visualization of Jz(x, y) of Run 2m6 with PrM = 10, Lu = 1.8 × 105, Luν ≈ 5 × 104, and 163842 mesh points at t = 464 for a small part of the domain with sizes 2.8ξM(t)×0.74ξM(t). The lengths of 100 δc, 2Lc, and ξM are indicated by horizontal white solid, dashed, and dotted lines, respectively. The thickness of the current sheet corresponds to about 3Δx ≈ 21 δc. In its proximity, there are also indications of ringing, indicated by the black circle.

Next, we compare the typical length and thickness of current sheets in Run 2m6 with the values defined by Uzdensky et al. (2010) for critical current sheets, L_c = Lu_c η/v_A and ${\delta }_{\mathrm{c}}=L_{\mathrm{c}}/{\mathrm{Lu}}_{\mathrm{c}}^{1/2}$, respectively. These lengths are indicated in Fig. 14. Here we used Lu_c = 2.5 × 10⁴ as the critical Lundquist number, which was found to be representative of all of our cases. We see that the current sheets in Fig. 14 have a length that is comparable to ξ_M ≈ 7 L_c, because Lu ≈ 7 Lu_c. The thickness of the current sheets is about 20 δ_c. They are marginally resolved with about 3Δx, where Δx = 2π/16384 ≈ 3.8 × 10⁻⁴ is the mesh spacing. Thus, the aspect ratio of thickness to length of the current sheets is about $20{\delta }_{\mathrm{c}}/7L_{\mathrm{c}}=(20/7){\mathrm{Lu}}_{\mathrm{c}}^{-1/2}\approx 0.02$. This is about three times the nominal value estimated by Uzdensky et al. (2010) for critical current sheets.

By comparison, Hosking & Schekochihin (2023) estimated for the aspect ratio ${\delta }_{\mathrm{c}}/{\xi }_{\mathrm{M}}={\mathrm{Lu}}_{\mathrm{c}}^{1/2}/{\mathrm{Lu}}_{\nu }$. For Run 2m6 with Lu_ν ≈ 5 × 10⁴, this yields 0.003, which is about six times smaller than our value of 0.02. While these discrepancies could perhaps be explained by the absence of nondimensional prefactors in the definitions of δ_c, we must also consider the possibility that this is simply a consequence of a lack of resolution.

Runs 2m2 and 2m6 have in common that they have the same resistivity and nearly the same value of Lu of about 1.8 × 10⁵, but Run 2m2 has a tenfold larger viscosity than Run 2m6, so Pr_M is increased from 10 to 100. We see that C_M has deceased by only about 15%, which is much less than what is expected if C_M was proportional to ${\mathrm{Pr}}_{\mathrm{M}}^{-1/2}$.

In Fig. 15 we show for Runs 2m2 and 2m6 magnetic field profiles, −B_x(y) and B_y(y), and velocity profiles, −u_x(y) and u_y(y), through a particular current sheet. We see that in Run 2m2 with Pr_M = 100, the profiles are much smoother, even though the value of Lu = 1.8 × 10⁵ is the same in both cases. Thus, the viscosity has a significant effect in smoothing the magnetic field. Even so, the effect on the value of C_M is small.

Fig. 15. Magnetic field profiles [−Bx(y) and By(y) in red] and velocity profiles [−ux(y) and uy(y) in blue] through the current sheet (in black) for Run 2m2 (through x = 0.372) and Run 2m6 (through x = 0.270).

7. Endpoints in the primordial evolutionary diagram

The evolution of primordial magnetic fields is usually displayed in an evolutionary diagram showing the comoving values of B_rms versus ξ_M or, similarly, v_A versus ξ_M. This dependence corresponds to a power law of the form v_A = v_A0(ξ_M/ξ_M0)^κ. Since ξ_M(t)∼t^q, we have $v_{\mathrm{A}}(t)\sim t^{-p/2}\sim {{\xi }_{\mathrm{M}}}^{-p/2q}$, so κ = p/2q = 5/4 for the Hosking scaling with p = 10/9 and q = 4/9.

Following Banerjee & Jedamzik (2004), the time t in Eq. (10) would be replaced by the age of the Universe at the time of recombination, t_rec. As emphasized by Hosking & Schekochihin (2023), this gives an implicit equation for the magnetic field at recombination with the solution B(t_rec)≈10^−8.5 G (ξ_M/1 Mpc) = 10^−14.5 G (ξ_M/1 pc) if C_M = 1. However, B(t_rec) would be much larger and ξ_M much smaller when C_M ≫ 1 is taken into account.

Under the hypothesis of fast reconnection owing to plasmoid instability (Bhattacharjee et al. 2009; Uzdensky et al. 2010), Hosking & Schekochihin (2023) estimated C_M as the square root of an effective cutoff value of about 10⁴ for the Lundquist number and an additional factor of ${\mathrm{Pr}}_{\mathrm{M}}^{1/2}$. Here, they estimated Pr_M ≈ 10⁷, so C_M = 10^5.5 and B(t_rec)≈10⁻³ G (ξ_M/1 Mpc) = 10⁻⁹ G (ξ_M/1 pc).

Our new results challenge the reliability of the anticipated dependence of C_M on Pr_M. With the current results at hand, Fig. 13 suggests that C_M never exceeds the value $3.7{\mathrm{Lu}}_{\mathrm{c}}^{1/4}\approx 47$, even for Pr_M as large as 100. Particularly important is the fact that this result is independent of Pr_M. Given that Hosking & Schekochihin (2023) used Pr_M ≈ 10⁷, which yielded an extra 10^3.5 factor, and therefore C_M = 10^5.5 in their estimate, our new findings imply that the resistivity effect is actually independent of Pr_M. Although this has so far only been verified for values of Pr_M ≤ 5, we suggest that a more accurate formula for the endpoints of the evolution with C_M = 47 ≈ 10^1.7 would be B(t_rec)≈10^−6.8 G (ξ_M/1 Mpc) = 10^−12.8 G (ξ_M/1 pc). These values are still above the lower limits inferred from suppressed GeV photon emission from the halos of blazars (Neronov & Vovk 2010).

8. Conclusions

The present results have shown that, up to the largest Lundquist numbers accessible to our present direct numerical simulations with 2048³ mesh points, the decay times depend on the resistivity. Only for our 2D simulations do we see evidence for a cutoff. The dependence of hydromagnetic turbulence properties on the value of the resistivity is unusual for fully developed turbulence. We wonder whether our results reflect just a peculiar property of decaying turbulence or whether there could also exist aspects of statistically stationary turbulence that depend on the microphysical resistivity. Possible examples of resistively controlled time dependences could include the time that is needed to develop the final saturated magnetic energy spectrum in kinetically forced turbulence, where the magnetic field emerges due to a dynamo action (Haugen et al. 2003; Schekochihin et al. 2004). A resistively slow adjustment phase is reminiscent of what occurs for helical magnetic fields (Brandenburg 2001), where it has also been possible to measure a weak resistivity dependence of the turbulent magnetic diffusivity (Brandenburg et al. 2008).

In the present 3D case, the magnetic helicity vanishes on the average. However, in the spirit of the Hosking phenomenology of a decay controlled by the conservation of the Hosking integral, it is very possible, even in decaying turbulence, that the conservation of magnetic helicity in patches of one sign of magnetic helicity plays an important role in causing the resistively controlled decay speed. Whether or not this is equivalent to talking about reconnection remains an open question. As discussed in Sect. 2, the idea of reconnection in terms of current sheets and plasmoids may not be fully applicable in the context of turbulence, where magnetic structures are more volume filling than in standard reconnection experiments (see Fig. 2). Additional support for a possible mismatch between classical reconnection theory and turbulent decay times comes from our numerical finding that the dependence of C_M on Lu seems to be independent of the value of Pr_M (see Fig. 13). More extensive numerical studies with resolutions of 8192³ meshpoints, which was the resolution needed to see a leveling off in 2D, might suffice to verify our 2D findings in the 3D case.

Our present work motivates possible avenues for future research. One choice is to do the same for turbulence with a −αu friction term in the momentum equation. Such calculations were already performed by Banerjee & Jedamzik (2004) and the friction term is also incorporated in the phenomenology of Hosking & Schekochihin (2023). This is to model the drag from photons when their mean-free path begins to exceed the scale ξ_M. This is the case after the time of recombination, when it contributes to dissipating kinetic energy and leads to a decoupling of the magnetic field. The field then becomes static and remains frozen into the plasma. According to the work of Hosking & Schekochihin (2023), this results in a certain reduction of C_M compared to the resistively limited value. Verifying this with simulations would be particularly important.

Another critical aspect to verify is the absence of a dependence of C_M on Pr_M over a broader range of parameter combinations. Given that there are always limitations on the resolution, it may be useful to explore simulations in rectangular domains to cover a larger range of scales. Other possibilities include simulations with time-dependent values of η and ν to obtain a larger separation between ξ_M and the dissipation scale near the end of the simulation. However, there is the danger that artifacts are introduced that need to be carefully examined. Ill-understood artifacts can also be introduced by using hyperviscosity and hyperresistivity, which are used in some simulations. For the time being, however, the possibility of a Pr_M dependence of C_M cannot be confirmed from our simulations. Whether or not this automatically rules out reconnection as the reason for a resistively limited value of C_M, instead of the idea of magnetic helicity conservation in smaller patches, remains uncertain.

As the value of Pr_M is increased, we also see a systematic increase in the kinetic-to-magnetic dissipation ratio, ϵ_K/ϵ_M. Such a dependence has previously been seen for kinetically dominated forced turbulence, but it is shown here, perhaps for the first time, for magnetically dominated decaying turbulence. While such results may be of interest to the problem of coronal heating, it should be remembered that the present simulations have large plasma betas (i.e., the gas pressure dominates over the magnetic pressure). It would therefore be of interest to check whether the obtained Pr_M-dependence also persists for smaller plasma beta. The restriction to two dimensions is another computational simplification that allowed significantly larger Lundquist numbers to be accessed, but it needs to be checked that the results for ϵ_K/ϵ_M are not very sensitive to this restriction. Comparing the 3D Runs M1 and M3 of Table 1 with the 2D Run 2M3 of Table 3, which have Pr_M = 5 and similar Lundquist numbers, we see that the 2D results may overestimate the ratio ϵ_K/ϵ_M by a factor of about two. Future work will need to show whether this can also be confirmed for other parameters.

Acknowledgments

We thank Pallavi Bhat (Bangalore) and David Hosking (Princeton) for detailed comments and suggestions, especially regarding reconnection in the 2D case. We are also acknowledge inspiring discussions with the participants of the program on “Turbulence in Astrophysical Environments” at the Kavli Institute for Theoretical Physics in Santa Barbara. We are also grateful to Sébastien Galtier, Romain Meyrand, Annick Pouquet, and Alex Schekochihin for their interest in shedding light on the origin of the word anastrophy; see Appendix A. This research was supported in part by the Swedish Research Council (Vetenskapsrådet) under Grant No. 2019-04234, the National Science Foundation under Grant No. NSF PHY-1748958, a NASA ATP Award 80NSSC22K0825, and the Munich Institute for Astro-, Particle and BioPhysics (MIAPbP), which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2094 – 390783311. We acknowledge the allocation of computing resources provided by the Swedish National Allocations Committee at the Center for Parallel Computers at the Royal Institute of Technology in Stockholm. F.V. acknowledges partial financial support from the Cariplo “BREAKTHRU” funds Rif: 2022–2088 CUP J33C22004310003. Software and data availability. The source code used for the simulations of this study, the PENCIL CODE (Pencil Code Collaboration 2021), is freely available on https://github.com/pencil-code/. The DOI of the code is https://doi.org/10.5281/zenodo.2315093. The simulation setups and corresponding data are freely available on https://doi.org/10.5531/sd.astro.9 (Brandenburg et al. 2024) and https://doi.org/10.5281/zenodo.10527437.

References

Appendix A: Historical note on anastrophy

In recent years, the term anastrophy for the mean squared magnetic vector potential $\langle A_z^2\rangle =\text{const}$ has become increasingly popular (Tronko et al. 2013; Galtier & Meyrand 2015; Zhou et al. 2021; Hosking & Schekochihin 2021; Schekochihin 2022). In the 1970s, it was referred to as mean square vector potential (Fyfe & Montgomery 1976) or as the variance of the magnetic potential (Pouquet 1978). The term anastrophy was first used in the 1987 Les Houches lecture notes by Pouquet (1993), and it was also used by Vakoulenko (1993), but neither of them provided an explanation of its origin.

Annick Pouquet (2024, private communication) informed us now that the word may have been invented by Uriel Frisch and Nicolas Papanicolaou during a meeting on a Winter Sunday in the late 1970s at Saint-Jean-Cap-Ferrat, while she and Jacques Léorat were also present.

The word has Greek roots, where ‘strophe’ refers to curl or turning, and ‘ana’ therefore hints at the inverted curl of B. This is somewhat reminiscent of the word palinstrophy, which is the mean squared double curl of the velocity, where ‘palin’ means again. This term was also invented by Frisch and Papanicolaou. The palinstrophy is proportional to the rate of change of enstrophy, i.e., the mean squared vorticity.

Appendix B: Resolving the plasmoid instability

The expected thickness of critical layers (i.e., the minimum thickness of current sheets when they become unstable to the plasmoid instability) is expected to be on the order of δ_c (see Sect. 6). To compare this with our 2D simulations, we show in Fig. 10 visualizations of J_z(x, y) in parts of the full domain with sizes 5ξ_M × 10ξ_M at times t = 1, 10, 100, and 464. We note that ξ_M(t) increases with time, so the x and y ranges increase with time too. We see that each panel contains about one pair of current tubes. Furthermore, the other length scales, L_c and δ_c, increase with time by the same factor, as expected for a self-similar evolution.

At late times, for t = 100 and 464, the current sheets are seen to break up into plasmoids. They are obviously better resolved in Run 2m6 with 16384² mesh points than in Run 2m5 with only 8192² mesh points. At lower resolution, there is a higher tendency for ringing (i.e., the formation of oscillations on the grid scale), which indicates that the resolution limit has been reached. Nevertheless, Table 3 and Fig. 13 show that, within the error bars, the values of C_M are similar for the two resolutions.

It is possible that the critical current sheets are not well resolved during a significant fraction of the duration of the simulation. Theoretically, we expect the thinnest sheets to set the reconnection rate, but this is also the place where the value of Pr_M matters because viscosity is unimportant for the larger scales in the plasmoid hierarchy. Therefore, if we do not resolve that sheet, it seems reasonable that we would also not see a dependence on Pr_M. To check this, we now consider a version of Run 2m2, where the viscosity is ten times larger, which increases Pr_M from 10 to 100; see Run 2m6. The result is shown in Fig. B.1. We see that a lower value of Lu for Pr_M = 100 changes the results in an expected fashion, thus rejecting the possibility that the weak dependence of C_M on Pr_M was an artifact of having chosen unreliably large values of Lu; compare Runs 2m1 and 2m2 in Fig. B.2.

Fig. B.1. Comparisons of visualizations of Jz(x, y) for Runs 2m5 and 2m6 with PrM = 10 and 81922 (left) and 163842 (right) mesh points at times t = 1, 10, 100, and 464. In each panel the lengths of the dotted, dashed, and solid lines denote the values of ξM, 5 Lc, and 500 δc. The bottom right panel shows the same current sheet that was presented in Fig. 14 as a blow-up. The thick white line in that panel at (x, y) = (0.27, 0.25) marks the location of the cross section shown in Fig. 15.

Fig. B.2. Similar to Fig. B.1, but for Run 2m1 with Lu = 75, 000 (left) and Run 2m2 with Lu = 182, 000 (right) using PrM = 100 and 163842 mesh points in both cases. The thick white line in the bottom right panel at (x, y) = (0.37, 0.25) marks the location of the cross section shown in Fig. 15.

In the bottom right panel of Fig. B.1, we see the same current sheet that was already presented in Sect. 6 as Fig. 14. However, we also see in Fig. B.1 that there are many other current sheets that are not yet in the process of breaking up.

All Tables

All Figures

	Fig. 1. Visualization of Bz on the periphery of the computational domain for Run M1, where PrM = 5.
In the text

	Fig. 2. Visualization of Jz at z = π along with a zoom-in on the lower left corner for Run M0, where PrM = 10, at t = 644. The value of ξM at that time is indicated by the length of the short white lines.
In the text

	Fig. 3. Magnetic and kinetic energy spectra for Run M1 at different times. For large values of k there is a k−2 inertial range. The ∝k4 and ∝k2 slopes are indicated for reference. As elsewhere, time is in units of [t]=(csk1)−1.
In the text

	Fig. 4. Approach of CM(t), Lu(t), and $C_{\mathrm{L}}^{(1/4)}(t)$ to an approximately constant value toward the end of the simulation for Run M1.
In the text

	Fig. 5. Dependence of CM on Lu and Luν. An approximate scaling ∝Lu1/4 is found for Lu < 6000 and a piecewise power-law scaling $C_{\text{M}}=C_{\text{L}}^{(n)}{\text{Lu}}^n$ with different prefactors for the runs with higher viscosity (5 × 10−6 red symbols and 4 × 10−6 for the green symbol) and lower viscosity (2 × 10−6 blue symbols).
In the text

	Fig. 6. Dependence of ℰM(t) and ξM(t) for three values of PrM.
In the text

	Fig. 7. Evolution of τℰ(t) and τξ(t) for three values of PrM with fixed and moderate viscosity. These times tend to be approximately constant at late times with values approximately consistent with those in Table 1.
In the text

	Fig. 8. Dependence of τξ and τℰ on Lu. The two show a similar dependence on Lu, approximately ∝Lu1/4, but there is possible evidence for a leveling off for larger values of Lu.
In the text

	Fig. 9. Dependence of CK on Re. The line Re0.1 is shown for comparison, but the data are also nearly compatible with being independent of Re.
In the text

	Fig. 10. Similar to Fig. 3, but for Run 2M1. The magnetic peaks lie underneath a kβ envelope with β = 1, as expected in the case of anastrophy conservation.
In the text

	Fig. 11. Magnetic energy spectra (left) and kinetic energy spectra (right) for Run 2m2 with Lu = 1.8 × 105 and PrM = 100, at times t = 5, 20, 100, and 400, collapsed on top of each other by plotting both vs kξM(t) and scaling them with ξM. This makes their heights agree, as expected.
In the text

	Fig. 12. Dependence of ϵK/ϵM on PrM. The solid line denotes $1.3{\mathrm{Pr}}_{\mathrm{M}}^{0.7}$, but many of the data points, predominantly those with smaller Lu, are beneath that line. The inset shows that for PrM = 100, ϵK/ϵM increases with Lu like Lu0.6.
In the text

Fig. 13. Dependence of CM on Lu for the 2D runs. Shallow scaling ∝Lu0.1 is found for 104 < Lu < 105, which is also compatible with a leveling off at Luc = 2.5 × 104, as described by Eq. (24) with n = 1/4. The black (red) data points are for PrM = 5 (PrM = 1). The blue data points denote the 3D results from Sect. 4. The orange symbols are for the runs with PrM = 10 and 20, listed in Table 3. The thin dotted line gives Eq. (24) with n = 1/4 for comparison.

In the text

Fig. 14. Visualization of Jz(x, y) of Run 2m6 with PrM = 10, Lu = 1.8 × 105, Luν ≈ 5 × 104, and 163842 mesh points at t = 464 for a small part of the domain with sizes 2.8ξM(t)×0.74ξM(t). The lengths of 100 δc, 2Lc, and ξM are indicated by horizontal white solid, dashed, and dotted lines, respectively. The thickness of the current sheet corresponds to about 3Δx ≈ 21 δc. In its proximity, there are also indications of ringing, indicated by the black circle.

In the text

	Fig. 15. Magnetic field profiles [−Bx(y) and By(y) in red] and velocity profiles [−ux(y) and uy(y) in blue] through the current sheet (in black) for Run 2m2 (through x = 0.372) and Run 2m6 (through x = 0.270).
In the text

Fig. B.1. Comparisons of visualizations of Jz(x, y) for Runs 2m5 and 2m6 with PrM = 10 and 81922 (left) and 163842 (right) mesh points at times t = 1, 10, 100, and 464. In each panel the lengths of the dotted, dashed, and solid lines denote the values of ξM, 5 Lc, and 500 δc. The bottom right panel shows the same current sheet that was presented in Fig. 14 as a blow-up. The thick white line in that panel at (x, y) = (0.27, 0.25) marks the location of the cross section shown in Fig. 15.

In the text

	Fig. B.2. Similar to Fig. B.1, but for Run 2m1 with Lu = 75, 000 (left) and Run 2m2 with Lu = 182, 000 (right) using PrM = 100 and 163842 mesh points in both cases. The thick white line in the bottom right panel at (x, y) = (0.37, 0.25) marks the location of the cross section shown in Fig. 15.
In the text