4.1. Value of $\mathcal {I}_H(R)$ from different methods

In figure 1 we compare for run K60D1c the temporal dependence of $\mathcal {I}_H(R)$ and $I_H$ computed from the methods introduced in § 2: BC and CI methods with cubic and spherical regions $V_R$ for each, and the fitting method. For the purpose of this comparison, a resolution of $1024^3$ mesh points suffices, although $I_H$ is not as well conserved as for higher resolutions. Our main results are obtained with $2048^3$ mesh points; see table 1. The produced $\mathcal {I}_H(R)$ and $I_H$ curves have different magnitudes, as expected from the fact that $w_R({\boldsymbol {k}})$ are different; see figures 1(a) and 1(b). However, the scaling properties, i.e. the values of $p_H=-\text {d}\ln I_H/\text {d}\ln t$ are unchanged among all the methods; see figure 1(c). In what follows, we use the CI method with cubic regions throughout.

Figure 1. Comparing results for run K60D1c in different methods; (a) $\mathcal {I}_H(R)$ at $t=0$. The vertical line indicates $R=0.115L$ with which we compute $I_H$. (b) Time evolution of $I_H$. (c) Time evolution of the decay exponents $p_H=-\text {d}\ln I _H/\text {d}\ln t$.

4.2. Gauge invariance

In figure 2(a) we compare the $R$-dependence of the correlation integral $\mathcal {I}_H(R)$ at different times under the resistive gauge (black solid) and the Coulomb gauge (red dashed). Noticeable differences appear only at later times and at small $R\ll \xi _{M}$, or when $R$ is close to the system scale $L$. The differences remain negligible in the asymptotic regime for all $t$. Thus, even though $h$ itself, and also its spectrum, are gauge dependent, $I_H$ is not (Hosking & Schekochihin Reference Hosking and Schekochihin2021).

Figure 2. Results for run K60D1c. (a) Comparing $\mathcal {I}_H(R)$ from different gauges. (b) The auto-correlation curves $C_h(R)$. The inset shows $4{\rm \pi} R^2C_h(R)$. Note that the abscissa is normalized by $\xi _{M}$, which is time dependent. For both panels, the pairs of curves are taken at $t=0$, $0.2$, $0.5$, $1.5$, $4.6$, $15$, $46$, $147$ in code units from top to bottom, as indicated by the arrows.

The evolution of the auto-correlation $C_h(R)=\left \langle h({\boldsymbol {x}})h({\boldsymbol {x}}+\boldsymbol {R}) \right \rangle$ can be computed from $(4{\rm \pi} R^2)^{-1}\,\text {d}\mathcal {I}_H(R)/\text {d}R$; see figure 2(b). Using a time-dependent normalization of the abscissa, we see that $h$ is always correlated at approximately $2\xi _{M}(t)$, as also evident from the inset.

Let us point out at this point that $\text {Sp}\left ( h \right )$ is not only gauge dependent, but it can provide some useful insight into the nature of gauge dependence. Candelaresi et al. (Reference Candelaresi, Hubbard, Brandenburg and Mitra2011) used this fact to show that the advective gauge, where $\varphi =\boldsymbol {u}\boldsymbol {\cdot }{\boldsymbol {A}}$, can lead to large gradients in the evolution of ${\boldsymbol {A}}$, which can cause fatal inaccuracies in the nonlinear regime.

4.3. Energy and helicity density spectra

In figure 3, we present magnetic energy spectra $E_{M}(k)\equiv \text {Sp}\left ( {\boldsymbol {B}} \right )/2\mu _0$ and magnetic helicity density spectra $\text {Sp}\left ( h \right )$ for run K60D1bt at different times. From the energy spectra, we see that the inertial range shifts both in $k$ and in amplitude. The slope in the inertial range becomes slightly steeper and closer to $k^{-2}$ at later times. This is in agreement with previous works (Brandenburg et al. Reference Brandenburg, Kahniashvili and Tevzadze2015) and may support the reconnection-controlled decay picture (Bhat et al. Reference Bhat, Zhou and Loureiro2021; Zhou et al. Reference Zhou, Wu, Loureiro and Uzdensky2021); see also Zhou et al. (Reference Zhou, Bhat, Loureiro and Uzdensky2019) for two-dimensional MHD, and Zhou et al. (Reference Zhou, Loureiro and Uzdensky2020) for reduced MHD systems. The initial $k^4$ subrange, however, does become slightly shallower at late times when the position of the peak of the spectrum has dropped below $k/k_1\approx 10$. This is presumably a finite-size effect, and so the late time evolution may not be reliable unless the initial $k_\text {peak}$ value was large enough and the total number of mesh points was sufficient to resolve the inertial and sub-inertial ranges reasonably well. Since we do not know a priori what are the quantitative requirements on the initial $k_\text {peak}$ and on the numerical resolution, we must regard our results with care and should discuss the possibility of artifacts as we reach the limits of what can be considered safe. Thus, we can conclude that, except for very late times, ${\mathcal {E}}_{M}$ has remained nearly self-similar. The magnetic helicity density spectrum has, just like $E_{M}(k)$, a $k^{-5/3}$ spectrum at the initial time for $k>k_\text {peak}$. At small $k$, however, $\text {Sp}\left ( h \right )$ has a random noise spectrum $\propto k^2$, which remains unchanged also at late times. This agrees with our expectations; see (2.14), which allows us to determine $I_H$ from the spectral value at small $k$.

Figure 3. Results for run K60D1bt. (a) Magnetic energy spectrum. (b) Spectrum of magnetic helicity density. (c) The non-dimensionalized and compensated $\text {Sp}\left ( h \right )$; see (4.11). (d) Time evolution of the non-dimensionalized $I_H$, (4.12) and (4.13). For the first three panels, the vertical grey lines mark the asymptotic scale chosen to be $k=2{\rm \pi} /(2R)=4.35$ with $R=0.115L$, at which value we have computed $\tilde {I}_H$ in (d).

4.4. Lundquist-number scaling of decay exponents

In figure 4(a), we present the time evolution of the normalized $I_H$ for all runs with $R=0.115L$. With increasing hyper-Lundquist number ${Lu}_n$ (cf. table 1), $I_H$ becomes progressively better conserved. The instantaneous decay exponents $p_H$ vs ${Lu}_n$ are plotted in figure 4(b) for different runs, along with the energy decay exponents $p$. While the latter only has a weak dependence on ${Lu}_n$ and approaches an asymptotic value close to unity (Brandenburg et al. Reference Brandenburg, Kahniashvili and Tevzadze2015; Brandenburg & Kahniashvili Reference Brandenburg and Kahniashvili2017; Reppin & Banerjee Reference Reppin and Banerjee2017; Bhat et al. Reference Bhat, Zhou and Loureiro2021; Zhou et al. Reference Zhou, Wu, Loureiro and Uzdensky2021), at large ${Lu}_n$, $p_H$ decreases and is found to scale approximately as ${Lu}_n^{-1/4}$.

Figure 4. (a) Time evolutions of the normalized $I_H$. (b) The instantaneous decay exponents of $I_H$ ($p_H$) and ${\mathcal {E}}_{M}$ ($p$) vs ${Lu}_n$, and the dash-dotted line indicates $p=10/9$. The size of the symbols increases with time. (c) Same as (b), but plotting for the decay exponent of the mean magnetic helicity ($p_{AB}$), and that of the Hosking cross-helicity integral ($p_C$). (d) Time evolution of $\mathcal {I}_H(R)$ for run K60D3bc. The vertical line indicates the asymptotic scale chosen to be $R/(2{\rm \pi} )=0.115$.

Note that the data points of $p_H$ at the largest ${Lu}_n$ values start to level off and even increase with decreasing ${Lu}_n$. We argue that this is an artefact due to the finite size of the computational domain, and not due to entering the fast-reconnection regime where the reconnection time scale becomes independent of ${Lu}_n$. In fact, for run K60D3bc, we have ${Lu}_n^{1/n}=20$ at the end of the simulation, which is still far away from the predicted transition value ${\sim }10^4$ (Loureiro et al. Reference Loureiro, Cowley, Dorland, Haines and Schekochihin2005, Reference Loureiro, Schekochihin and Cowley2007).

To provide evidence of the limitation from a finite-size box, in figure 4(d), we plot the time evolution of the CI $\mathcal {I}_H(R)$. At small $R\lesssim \xi _{M}$, we see an $R^3$ scaling. Thus, $\left \langle H_{V_R}^2 \right \rangle /V_R$ is proportional to $V_R\propto R^3$, and therefore $\left \langle H_{V_R}^2 \right \rangle$ is proportional to $V_R^2$, as expected for a nearly space-filling distribution of magnetic helicity patches of small scale.Footnote ⁹ At the end of the simulation, $\xi _{M}$ has become comparable to the asymptotic scale chosen (vertical line), due to which $I_H$ exhibits an accelerating decay; see the data points of $p_H$ at the largest ${Lu}_n$ values in figure 4(b).

It is worth noting that the mean magnetic helicity density $|H_V|$ is never exactly zero, but it is instead several orders of magnitude below its maximum value. Interestingly, the decline of the instantaneous decay coefficient of $|H_V|$, referred to here as $p_{AB}$, follows a similar decline with ${Lu}_n$ as $p_H$; see figure 4(c).

Following Hosking & Schekochihin (Reference Hosking and Schekochihin2021), we have also performed a similar analysis for the cross-helicity integral, defined as

(4.1)\begin{equation} \mathcal{I}_c(R)=\int_{V_R}\text{d}^3r\left\langle h_{c}({\boldsymbol{x}})h_{c}({\boldsymbol{x}}+\boldsymbol{r}) \right\rangle, \end{equation}

where $h_{c}=\boldsymbol {u}\boldsymbol {\cdot }{\boldsymbol {b}}$ is the cross helicity density. This is motivated by the fact that the cross helicity is also an ideal invariant of MHD (Woltjer Reference Woltjer1958). We found that the asymptotic limit, $I_C$, of $\mathcal {I}_C$, also decays in time, just like $I_H$, but the decay exponent scales with ${Lu}_n$ similarly to that of the magnetic energy density; see figure 4(c), where $p_C=-\text {d}\ln I_C/\text {d}\ln t$. Hence, $I_C$ is not as well conserved as $I_H$.

4.5. Non-dimensionalizing the Hosking integral

Since $I_H$ is conserved, its value can be estimated from the initial condition. Through the definition of $h$ we have

(4.2)\begin{equation} \langle \tilde{h}^*({\boldsymbol{k}}_1)\tilde{h}({\boldsymbol{k}}_2) \rangle =\int\frac{\text{d}^3k'}{(2{\rm \pi})^3}\,\frac{\text{d}^3k''}{(2{\rm \pi})^3} \epsilon_{ijk}\epsilon_{abc}k'_k k''_c \langle \tilde{A}^*_i({\boldsymbol{k}}_1-{\boldsymbol{k}}')\tilde{A}^*_j({\boldsymbol{k}}') \tilde{A}_a({\boldsymbol{k}}_2-{\boldsymbol{k}}'')\tilde{A}_b({\boldsymbol{k}}'') \rangle. \end{equation}

In the Coulomb gauge, we have, for an isotropic field

(4.3)\begin{equation} \langle \tilde{A}^*_i({\boldsymbol{k}}_1)\tilde{A}_j({\boldsymbol{k}}_2)\rangle =(2{\rm \pi})^3\delta^3({\boldsymbol{k}}_1-{\boldsymbol{k}}_2)P_{ij}({\boldsymbol{k}}_1)M_A(k_1), \end{equation}

where $P_{ij}({\boldsymbol {k}})=\delta _{ij}-k_ik_j/k^2$, and $M_A$ is related to the magnetic energy spectrum via $M_A(k)=2{\rm \pi} ^2 E_{M}(k)/k^4$. Taking advantage of the fact that $\tilde {A}_i({\boldsymbol {k}})$ is a Gaussian field at $t=0$, we can decompose the four-point correlations into products of two-point correlations. For any wave vectors ${\boldsymbol {k}}_1$, ${\boldsymbol {k}}_2$, ${\boldsymbol {k}}_3$ and ${\boldsymbol {k}}_4$

(4.4)\begin{align} & \langle \tilde{A}^*_i({\boldsymbol{k}}_1)\tilde{A}^*_j({\boldsymbol{k}}_2) \tilde{A}_a({\boldsymbol{k}}_3)\tilde{A}_b({\boldsymbol{k}}_4) \rangle\nonumber\\ & \quad =\langle \tilde{A}^*_i({\boldsymbol{k}}_1)\tilde{A}_j(-{\boldsymbol{k}}_2) \rangle\langle \tilde{A}^*_a(-{\boldsymbol{k}}_3)\tilde{A}_b({\boldsymbol{k}}_4) \rangle +\langle \tilde{A}^*_i({\boldsymbol{k}}_1)\tilde{A}_a({\boldsymbol{k}}_3) \rangle\langle \tilde{A}^*_j({\boldsymbol{k}}_2)\tilde{A}_b({\boldsymbol{k}}_4) \rangle\nonumber\\ & \qquad +\langle \tilde{A}^*_i({\boldsymbol{k}}_1)\tilde{A}_b({\boldsymbol{k}}_4) \rangle\langle \tilde{A}^*_j({\boldsymbol{k}}_2)\tilde{A}_a({\boldsymbol{k}}_3) \rangle\nonumber\\ & \quad =(2{\rm \pi})^6[\delta^3({\boldsymbol{k}}_1+{\boldsymbol{k}}_2)\delta^3({\boldsymbol{k}}_3+{\boldsymbol{k}}_4) P_{ij}({\boldsymbol{k}}_1)P_{ab}({\boldsymbol{k}}_3)M_A(k_1)M_A(k_3)\nonumber\\ & \qquad +\delta^3({\boldsymbol{k}}_1-{\boldsymbol{k}}_3)\delta^3({\boldsymbol{k}}_2-{\boldsymbol{k}}_4) P_{ia}({\boldsymbol{k}}_1)P_{jb}({\boldsymbol{k}}_2)M_A(k_1)M_A(k_2)\nonumber\\ & \qquad +\delta^3({\boldsymbol{k}}_1-{\boldsymbol{k}}_4)\delta^3({\boldsymbol{k}}_2-{\boldsymbol{k}}_3) P_{ib}({\boldsymbol{k}}_1)P_{ja}({\boldsymbol{k}}_2)M_A(k_1)M_A(k_2)]. \end{align}

Together with (2.10), we find

(4.5)\begin{equation} \text{Sp}\left( h \right)=\frac{k^2}{2}\int_0^\infty\text{d}k'\int_{{-}1}^1\text{d}\alpha\, \frac{k'^2+\mu^2k'^2-2\mu k'|{\boldsymbol{k}}-{\boldsymbol{k}}'|}{k'^2|{\boldsymbol{k}}-{\boldsymbol{k}}'|^4} E_{M}(|{\boldsymbol{k}}-{\boldsymbol{k}}'|)E_{M}(k'), \end{equation}

where $\alpha =\cos ({\boldsymbol {k}},{\boldsymbol {k}}')$, and

(4.6a,b)\begin{equation} |{\boldsymbol{k}}-{\boldsymbol{k}}'|=\sqrt{k^2+k'^2-2\alpha k k'},\quad \mu=\frac{\alpha k' k-k'^2} {k'\sqrt{k^2+k'^2-2\alpha k k'}}. \end{equation}

Note that (4.5) is similar to the expression for $\text {Sp}\left ( {\boldsymbol {B}}^2 \right )$; see (27) of Brandenburg & Boldyrev (Reference Brandenburg and Boldyrev2020).

Consider a piecewise power-law spectrum,

(4.7)\begin{equation} E_{M}(k)=\left\{\begin{aligned} & E_\text{peak}\left(k/k_\text{peak}\right)^4, & k\leqslant k_\text{kpeak} \\ & E_\text{peak}\left(k/k_\text{peak}\right)^{{-}s}, & k> k_\text{kpeak} \end{aligned}\right\}, \end{equation}

where $s=5/3$ and $s=2$ correspond to the inertial range slopes at early and late times, respectively. For $s=5/3$ we have the relations

(4.8a,b)\begin{equation} E_\text{peak}=\frac{20}{17}{\mathcal{E}}_{M}\xi_{M},\quad k_\text{peak}=\frac{1}{2\xi_{M}}. \end{equation}

The leading-order term ($\propto k^2$) of (4.5) can be readily obtained by taking $k=0$ in the integrand, which gives

(4.9)\begin{equation} \text{Sp}\left( h \right)^\text{early}=\frac{136 \,E_\text{peak}^2}{95\,k_\text{peak}^3}k^2 +\mathcal{O}(k^3) =\frac{5120{\mathcal{E}}_{M}^2\xi_{M}^5}{323}k^2+\mathcal{O}(k^3). \end{equation}

For $s=2$ we obtain

(4.10)\begin{equation} \text{Sp}\left( h \right)^\text{late}=\frac{131072{\mathcal{E}}_{M}^2\xi_{M}^5}{13\,125}k^2+\mathcal{O}(k^3). \end{equation}

In figure 3(c), we show the compensated spectra normalized using $s=5/3$,

(4.11)\begin{equation} \widetilde{\text{Sp}}(h)=\frac{323}{5120\,{\mathcal{E}}_{M}^2 \xi_{M}^5k^2}\text{Sp}\left( h \right), \end{equation}

which is indeed $\sim \mathcal {O}(1)$ initially at small $k$, but then increases to $\mathcal {O}(10)$ at later times. Had we used $s=2$, $\widetilde {\text {Sp}}(h)$ only obtains larger values because the numerical factor in (4.10) is smaller than that in (4.9). Indeed, the increase in $\widetilde {\text {Sp}}(h)$ could be caused by (i) the Lundquist number not being sufficiently high and/or (ii) the magnetic field becoming non-Gaussian at later times. To quantify the latter, we find the excess kurtosis of the magnetic field, defined as $-3+\sum _{i=1}^3 \left \langle B_i^4 \right \rangle /\left \langle B_i^2 \right \rangle ^2/3$, to be $\sim -0.24$ at the end of run K60D1bt. This is slightly larger in magnitude compared with earlier work (Brandenburg & Boldyrev Reference Brandenburg and Boldyrev2020). The influence from non-Gaussianity can be explicitly checked by computing the right-hand side of (4.5) and comparing with $\text {Sp}\left ( h \right )$,Footnote ¹⁰ which is shown in figure 5. The Gaussianity is very well satisfied at $t=0$ since we have initialized the field to be so, but it becomes poor at late times. In particular, the right-hand side of (4.5), which approximates $\text {Sp}\left ( h \right )$ using $E_{M}$, under-estimates the true value of $\text {Sp}\left ( h \right )$ by nearly one order of magnitude, although sharing a similar spectral shape with the latter. Hence, the non-dimensionalized $\widetilde {\text {Sp}}(h)$ would be larger than unity also by approximately one order of magnitude, in agreement with figure 3(c).

Figure 5. For run K60D1bt, comparing the left-hand (solid black) and right-hand (dashed red) sides of (4.5), at two different snapshots, $t=0$ (upper pair) and $t=1$ (bottom pair).

Furthermore, we cannot rule out the possibility that the initial non-self-similar evolution of the magnetic field also contributes to the increase of $\widetilde {\text {Sp}}(h)$. Since the field is initialized to be Gaussian, local patches of magnetic fields with the same sign of magnetic helicity are likely not fully helical. Thus the field configuration close to $t=0$ has an effectively smaller $\xi _{M}$, and our estimation of $\widetilde {\text {Sp}}(h)$ will be lower at early time.

The magnitude of $I_H$ can be estimated from (2.15), and a proper normalization is

(4.12)\begin{equation} \tilde{I}_H=\frac{323}{10240{\rm \pi}^2\,{\mathcal{E}}_{M}^2 \xi_{M}^5}I_H \simeq0.0032\,{\mathcal{E}}_{M}^{{-}2} \xi_{M}^{{-}5}I_H, \end{equation}

with $s=5/3$ for early times, and

(4.13)\begin{equation} \tilde{I}_H=\frac{13\,125}{262144{\rm \pi}^2\,{\mathcal{E}}_{M}^2 \xi_{M}^5}I_H \simeq0.0051\,{\mathcal{E}}_{M}^{{-}2} \xi_{M}^{{-}5}I_H, \end{equation}

with $s=2$ for late times. Both curves are plotted in figure 3(d). The increasing value again highlights the non-Gaussianity of our simulations.

4.6. Evolution in a $pq$ diagram

To put our results into perspective, it is instructive to inspect the evolution in a $pq$ diagram; see § 3.2. In figure 6, we plot $pq$ diagrams for four representative runs. In each panel, the size of the symbols increases with time; the solid line corresponds to the Alfvén relation (3.9), and the dashed line gives the reconnection relation (3.10) for each run.

Figure 6. Panels (a–d) are for runs K200D3t, K60D1c, K60D1bt and K60D1bc, respectively. The symbol size increases with time. The dotted, solid and dashed lines are determined by (3.8), (3.9) and (3.10), respectively.

The most reliable runs are those in figures 6(c) and 6(d), where $N^3=2048^3$ mesh points have been used. We clearly see that the solution evolves along the $\beta =3/2$ line, as expected when the decay is governed by the Hosking integral. For a constant value of $\eta _1$, figure 6(d), the solution also reaches the line $p=2(1-q)$, which is referred to as the scale-invariance line. Note that in some earlier work, it was referred to as the self-similarity line; see the end of Appendix B.3 for a discussion. For a time-dependent $\eta _1(t)$, figure 6(c), the solution settles on a point below the scale-invariance line and is only slightly above the reconnection line. This might be indicative of reconnection playing indeed a certain role in the present simulations.

For our hyperviscous run K200D3t, the solution settles at much smaller values of $p$ and $q$ and is closer to the reconnection line than to the scale-invariance line, but the agreement in this case is not very good either. How conclusive this is in supporting the idea that reconnection plays a decisive role must therefore remain open. Nevertheless, also this solution lies close to the $\beta =3/2$ line, supporting again the governing role of $I_H$.

To compare with the case of a fully helical turbulent magnetic field, we refer to the $pq$ diagram in figure 2(c) of Brandenburg & Kahniashvili (Reference Brandenburg and Kahniashvili2017). There, the system evolved along the $\beta =0$ line toward the point $(p,q)=(2/3,2/3)$. They also considered the non-helical case in their figure 2(b), where the system evolved along the $\beta =1$ line toward the point $(p,q)=(1,1/2)$. In a separate study at lower resolution and with a different initial magnetic field, Brandenburg et al. (Reference Brandenburg, Kahniashvili, Mandal, Pol, Tevzadze and Vachaspati2017) found an evolution along the $\beta =2$ line toward $(p,q)=(6/5,2/5)$. Both results were arguably still consistent with an evolution along $\beta =3/2$ toward $(p,q)=(10/9,4/9)$.