Scaling of the Hosking integral in decaying magnetically dominated turbulence

Abstract

The Saffman helicity invariant of Hosking & Schekochihin (Phys. Rev. X, vol. 11, issue 4, 2021, 041005), which we here call the Hosking integral, has emerged as an important quantity that may govern the decay properties of magnetically dominated non-helical turbulence. Using a range of different computational methods, we confirm that this quantity is indeed gauge invariant and nearly perfectly conserved in the limit of large Lundquist numbers. For direct numerical simulations with ordinary viscosity and magnetic diffusivity operators, we find that the solution develops in a nearly self-similar fashion. In a diagram quantifying the instantaneous decay coefficients of magnetic energy and integral scale, we find that the solution evolves along a line that is indeed suggestive of the governing role of the Hosking integral. The solution settles near a line in this diagram that is expected for a self-similar evolution of the magnetic energy spectrum. The solution will settle in a slightly different position when the magnetic diffusivity decreases with time, which would be compatible with the decay being governed by the reconnection time scale rather than the Alfvén time.

1. Introduction

The subject of decaying turbulence plays important roles in laboratory and engineering applications (Proudman & Reid 1954; Stalp, Skrbek & Donnelly 1999), including superfluid (Nore, Abid & Brachet 1997) and supersonic ones (Kitsionas et al. 2009), as well as those where decaying temperature fluctuations are of interest (Warhaft & Lumley 1978), and in many areas of astrophysics ranging from star formation (Mac Low et al. 1998) to solar physics (Krause & Rüdiger 1975) and especially the early Universe (Christensson, Hindmarsh & Brandenburg 2001; Campanelli 2007; Kahniashvili et al. 2010; Brandenburg, Kahniashvili & Tevzadze 2015). The application to the early Universe focusses particularly on the decay of magnetic fields and their associated increase of typical length scales during the radiation-dominated epoch from microphysical to galactic scales (Brandenburg, Enqvist & Olesen 1996); see also Durrer & Neronov (2013), Subramanian (2016) and Vachaspati (2021) for reviews.

The properties of stationary Kolmogorov turbulence are governed by a constant flux of energy to progressively smaller scales. In decaying turbulence, however, this energy flux is time dependent. The rate of energy decay is governed by certain conservation laws, such as the Loitsiansky integral (Proudman & Reid 1954) or the Saffman integral (Saffman 1967); see Davidson (2000) for a review. As the energy decreases, the kinetic energy spectrum declines, primarily at large wavenumbers, causing the peak of the spectrum to move toward smaller wavenumbers or larger length scales.

The presence of a conservation law can also cause an inverse cascade. A famous example in magnetohydrodynamic (MHD) turbulence is magnetic helicity conservation, which can lead to a very pronounced increase of the typical length scale (Frisch et al. 1975; Pouquet, Frisch & Leorat 1976). In an alternative approach by Olesen (1997), it has been argued that the slope of the initial energy spectrum determines the temporal evolution of the spectrum. Olesen (1997) found the possibility of an inverse cascade for a wide range of initial spectral slopes. His argument was based on the observation that the MHD equations are invariant under rescaling space and time coordinates, ${\boldsymbol {x}}$ and $t$, respectively. Here, the relation between spatial and temporal rescalings depends on the initial spectrum. Brandenburg & Kahniashvili (2017) found that this relation can, instead, also be governed by the presence of a conservation law. In both cases, the formalism of Olesen (1997) predicts a self-similar behaviour of the energy spectrum. This can imply an increase of spectral energy at large length scales, i.e. an inverse cascade, or at least inverse transfer.¹

The presence of a conservation law can only affect the behaviour of the system if the conserved quantity is actually finite. Thus, even though magnetic helicity is always conserved at large magnetic Reynolds numbers, it may not play a role if the magnetic helicity is zero. However, MHD turbulence always has non-vanishing fluctuations of magnetic helicity. Hosking & Schekochihin (2021) have shown that the asymptotic limit, $I_H$, of the integral of the two-point correlation function of the local magnetic helicity density $h({\boldsymbol {x}},t)={\boldsymbol {A}}\boldsymbol {\cdot }{\boldsymbol {B}}$, is invariant in the ideal (non-resistive) limit and also independent of the gauge ${\boldsymbol {A}}\to {\boldsymbol {A}}'={\boldsymbol {A}}-{\boldsymbol {\nabla }}\varLambda$ for any scalar $\varLambda$. Here, ${\boldsymbol {B}}={\boldsymbol {\nabla }}\times {\boldsymbol {A}}$ is the magnetic field expressed in terms of the magnetic vector potential ${\boldsymbol {A}}$. Hosking & Schekochihin (2021) argued that the conservation of $I_H$ determines the decay of non-helical MHD turbulence, and together with self-similarity, it leads to the inverse cascading in non-helical magnetically dominated decaying turbulence (Brandenburg et al. 2015), which was found independently for relativistic turbulence (Zrake 2014). Brandenburg et al. (2015) argued that the decay was compatible with the conservation of anastrophy, i.e. the mean squared vector potential, but this explanation remained problematic owing to the gauge dependence of ${\boldsymbol {A}}$. Subsequent work by Brandenburg et al. (2017) with a different initial condition and lower resolution ($1152^3$ instead of $2304^3$) found a somewhat faster decay compatible with the conservation of the standard Saffman integral, which is the two-point correlation function of momentum, which is proportional to the velocity $\boldsymbol {u}$. This discrepancy suggests a possible dependence on the magnetic Reynolds and Lundquist numbers.

Hosking & Schekochihin (2021) studied the conservation properties of $I_H$ using also hyperviscosity and magnetic hyperdiffusivity to different degrees to study the dependence on the magnetic Reynolds and Lundquist numbers of the simulation, as well as justifying the role of the Sweet–Parker regime of magnetic reconnection. The latter arises because topological constraints on the magnetic field prevent ideal relaxation, as was first suggested by Zhou et al. (2019), Zhou, Loureiro & Uzdensky (2020) and Bhat, Zhou & Loureiro (2021). Hosking & Schekochihin (2021) have shown that the two physical models, Alfvén- vs reconnection-controlled decay, can be unambiguously distinguished by studying how the energy decay law scales with the hyper-diffusion order, $n$, and they found the reconnection time scale to be the relevant one. Hosking & Schekochihin (2022) argue that it prolongs the effective magnetic decay time and makes primordial magnetogenesis models consistent with the bound from GeV observations of blazars (Neronov & Vovk 2010).

The goal of the present paper is to provide independent evidence for the conservation properties of $I_H$ due to magnetic helicity fluctuations, i.e. in the absence of net magnetic helicity. We use a range of different methods to assess the reliability of the results at different stages during the decay. A particular difficulty is to determine the relevant length scale at which the Hosking integral is to be evaluated. We begin by reviewing its definition in § 2, discuss then our numerical stimulation set-up in § 3 and present our results in § 4. We conclude in § 5.

2. Expressions for the Hosking integral

In this section we introduce different expressions for the Hosking integral,² which we compare later in § 4.1.

2.1. Definition of the Hosking integral

We recall that the Hosking integral, $I_H$, is defined analogously to the Saffman integral in hydrodynamics, and emerges as the asymptotic limit of the integral of the two-point correlation function of the local magnetic helicity density $h({\boldsymbol {x}},t)={\boldsymbol {A}}\boldsymbol {\cdot }{\boldsymbol {B}}$. The latter, which we call $\mathcal {I}_H(R)$, is given by (Hosking & Schekochihin 2021)

(2.1)\begin{equation} \mathcal{I}_H(R)=\int_{V_R}\text{d}^3r\left\langle h({\boldsymbol{x}})h({\boldsymbol{x}}+\boldsymbol{r}) \right\rangle. \end{equation}

Here, the angle brackets denote an ensemble average, and $V_R$ is some volume with length scale $R$. We define the length scale of magnetic fluctuations as

(2.2)\begin{equation} \xi_{M}=\frac{\displaystyle\int\text{d}k\,k^{{-}1}E_{M}(k)} {\displaystyle\int\text{d}k\,E_{M}(k)}, \end{equation}

where $E_{M}(k)$ is the magnetic energy spectrum. An alternative length scale can be computed from the spectrum of the magnetic helicity density, $\xi _h$. A comparison between $\xi _{M}$ and $\xi _h$ is made in Appendix A, where we show that, although the ratio $\xi _{M}/\xi _h$ is not exactly a constant in time, it evolves much more slowly than the magnetic energy and is of order unity.³

When $R$ is much smaller than $\xi _{M}$, we have

(2.3)\begin{equation} \mathcal{I}_H(R\ll\xi_{M})\simeq \int_{V_R}\text{d}^3r\left\langle h({\boldsymbol{x}})h({\boldsymbol{x}}) \right\rangle\propto R^3, \end{equation}

whereas when $R\gg \xi _{M}$, but at the same time much smaller than the system scale $L$, $\mathcal {I}_H(R)$ reaches a constant asymptotic value, $I_H$, independent of $R$. Assuming that the volume average approximates the ensemble average (Hosking & Schekochihin 2021), we have

(2.4a,b)\begin{equation} I_H=\mathcal{I}_H(\xi_{M}\ll R\ll L)=\frac{1}{V_R}\left\langle H_{V_R}^2 \right\rangle,\quad H_{V_R}=\int_{V_R}\text{d}^3r\,h(\boldsymbol{r}). \end{equation}

Within this asymptotic range, $\mathcal {I}_H(R)$ is finite and independent of $R$, because the variance of the magnetic helicity $H_{V_R}$ contained in $V_R$ is expected to scale like $\left \langle H_{V_R}^2 \right \rangle \propto V_R^2(V_R/\xi _{M}^3)^{-1}\propto V_R$ (Hosking & Schekochihin 2021).

2.2. The box-counting method

To evaluate (2.4a,b) in numerical simulations, we again replace the ensemble average by a moving average, and hence calculate

(2.5)\begin{equation} \mathcal{I}_H(R)=\frac{1}{VV_R}\int_V\text{d}^3x\left[\int_{V_R} \text{d}^3r\,h({\boldsymbol{x}}+\boldsymbol{r})\right]^2, \end{equation}

where $V$ is the volume of the simulation box. In analogy to the algorithm for calculating the Hausdorff dimension, we call this the box-counting (BC) method. Upon a Fourier transformation,⁴  (2.5) can be recast as a weighted integral

(2.6)\begin{equation} \mathcal{I}_H(R)=\frac{1}{V}\int\frac{\text{d}^3k}{(2{\rm \pi})^3}\,w_R({\boldsymbol{k}}) h^*({\boldsymbol{k}})h({\boldsymbol{k}}). \end{equation}

The weight function $w_R({\boldsymbol {k}})$ depends on the shape of $V_R$. For a cubic region with length $2R$, we have

(2.7)\begin{equation} w_R({\boldsymbol{k}})=w_\text{cube}^\text{BC}({\boldsymbol{k}})\equiv8R^3\prod_{i=1}^3 {\rm j}_0^2(k_{\rm i}R), \end{equation}

whereas for a spherical region with radius $R$,

(2.8)\begin{equation} w_R({\boldsymbol{k}})=w_\text{sph}^\text{BC}(k)\equiv\frac{4{\rm \pi} R^3}{3} \left[\frac{6{\rm j}_1(kR)}{kR}\right]^2. \end{equation}

Here, ${\rm j}_0(x)=\sin x/x$ and ${\rm j}_1(x)=(\sin x-x\cos x)/x^2$ are the first- and second-order spherical Bessel functions, respectively. Note also that in (2.8), $w_\text {sph}^\text {BC}(k)$ depends only on $k=|{\boldsymbol {k}}|$. This allows us to rewrite (2.6) as just a one-dimensional integral

(2.9)\begin{equation} \mathcal{I}_H(R)=\int_0^\infty\text{d}k\,w_\text{sph}^\text{BC}(k)\,\text{Sp}\left( h \right), \end{equation}

where

(2.10)\begin{equation} \text{Sp}\left( h \right)=\frac{1}{V}\frac{k^2}{(2{\rm \pi})^3}\int_{|{\boldsymbol{k}}|=k}\text{d}\varOmega_k \, \tilde h^*({\boldsymbol{k}})\tilde{h}({\boldsymbol{k}}) \end{equation}

is the Fourier spectrum of the magnetic helicity density and $\varOmega _k$ is the solid angle in Fourier space, normalized such that $\int \text {d}k\,\text {Sp}\left ( h \right )=\left \langle h^2 \right \rangle$.⁵ The magnetic and kinetic energy spectra can then be written analogously as $E_{M}(k)=\text {Sp}\left ( {\boldsymbol {B}} \right )/2\mu _0$ and $E_{K}(k)=\rho _0\,\text {Sp}\left ( \boldsymbol {u} \right )/2$, where $\mu _0$ is the magnetic permeability and $\rho _0$ is the mean density.

2.3. The correlation-integral method

As an alternative to the BC method, (2.1) can be computed straightforwardly by approximating the ensemble average as a volume average, i.e.

(2.11)\begin{equation} \left\langle h({\boldsymbol{x}})h({\boldsymbol{x}}+\boldsymbol{r}) \right\rangle=\frac{1}{V}\int_V\text{d}^3x\,h({\boldsymbol{x}})h({\boldsymbol{x}}+\boldsymbol{r})= \frac{1}{V}\int\frac{\text{d}^3k}{(2{\rm \pi})^3}\,\tilde h^*({\boldsymbol{k}})\tilde{h}({\boldsymbol{k}})\,{\rm e}^{{\rm i}{\boldsymbol{k}}\boldsymbol{\cdot}\boldsymbol{r}}, \end{equation}

which we call the correlation-integral (CI) method. $\mathcal {I}_H(R)$ can then again be recast in the form of (2.6), but the weight functions are now slightly different. For a cubic region $V_R$ we have

(2.12)\begin{equation} w_R({\boldsymbol{k}})=w_\text{cube}^\text{CI}({\boldsymbol{k}})\equiv8R^3\prod_{i=1}^3 {\rm j}_0(k_{\rm i} R), \end{equation}

and for a spherical region we have

(2.13)\begin{equation} w_R({\boldsymbol{k}})=w_\text{sph}^\text{CI}(k)\equiv\frac{4{\rm \pi} R^3}{3} \frac{6{\rm j}_1(kR)}{kR}. \end{equation}

Note that $w_\text {sph}^\text {CI}(k)$ can also be used in (2.9).

2.4. The fitting method

A third way of obtaining $I_H$ is by extracting the coefficient of the leading-order term of the spectrum of $h$. At small $k\ll \xi _{M}^{-1}$ (Hosking & Schekochihin 2021)

(2.14)\begin{equation} \text{Sp}\left( h \right)=\frac{I_H}{2{\rm \pi}^2}k^2+\mathcal{O}(k^4), \end{equation}

given that $h({\boldsymbol {k}})$ is statistically isotropic, and therefore

(2.15)\begin{equation} I_H=\lim_{k\to0}\frac{2{\rm \pi}^2}{k^2}\text{Sp}\left( h \right). \end{equation}

3. Application to decaying MHD turbulence simulations

3.1. Basic equations

We solve the compressible MHD equations in a cubic domain of size $L^3$ using an isothermal equation of state with constant sound speed $c_{s}$, so the gas pressure is $p=\rho c _{s}^2$, where $\rho$ is the density. We allow for the possibility of hyperviscous and hyperdiffusive dissipation of kinetic and magnetic energies and solve the equations for ${\boldsymbol {A}}$ in the resistive gauge. The full set of equations is

(3.1)\begin{gather} \frac{{\rm D}\ln\rho}{{\rm D}t}={-}{\boldsymbol{\nabla}}\boldsymbol{\cdot}\boldsymbol{u}, \end{gather}

(3.2)\begin{gather}\frac{{\rm D}\boldsymbol{u}}{{\rm D}t}={-}c_{s}^2{\boldsymbol{\nabla}}\ln\rho +\frac{1}{\rho}[{\boldsymbol{J}}\times{\boldsymbol{B}}+{\boldsymbol{\nabla}}\boldsymbol{\cdot} (2\rho\nu_n\nabla^{2(n-1)}\boldsymbol{\mathsf{S}})], \end{gather}

(3.3)\begin{gather}\frac{\partial{\boldsymbol{A}}}{\partial t}=\boldsymbol{u}\times{\boldsymbol{B}} +\eta_n\nabla^{2n}{\boldsymbol{A}}, \end{gather}

where ${\boldsymbol {J}}={\boldsymbol {\nabla }}\times {\boldsymbol {B}}/\mu _0$ is the current density, $\nu _n$ is the viscosity, $\eta _n$ is the magnetic diffusivity, ${\textsf {{S}}}_{ij}=(\partial _i u_j+\partial _j u_i)/2-\delta _{ij}{\boldsymbol {\nabla }}\boldsymbol {\cdot }\boldsymbol {u}/3$ are the components of the rate-of-strain tensor $\boldsymbol{\mathsf{S}}$ and $n$ denotes the degree of hyperviscosity or hyperdiffusivity, with $n=1$ corresponding to ordinary viscous diffusive operators in (3.2) and (3.3), respectively. Equation (3.3) is here formulated in the resistive gauge, i.e. the scalar potential is $\varphi =-\eta _n\nabla ^{2(n-1)}{\boldsymbol {\nabla }}\boldsymbol {\cdot }{\boldsymbol {A}}$ in the uncurled induction equation $\partial {\boldsymbol {A}}/\partial t=-{\boldsymbol {E}}-{\boldsymbol {\nabla }}\varphi$, where ${\boldsymbol {E}}=-\boldsymbol {u}\times {\boldsymbol {B}}+\eta _n\nabla ^{2(n-1)}\mu _0{\boldsymbol {J}}$ is the electric field.

In some cases, we assume $\nu _n$ and $\eta _n$ to be time dependent, which is indicated by writing $\nu _n(t)$ and $\eta _n(t)$, respectively. In those cases, we assume a power-law variation (for $t>\tau$),

(3.4a,b)\begin{equation} \nu_n(t)=\nu_n^{(0)}[\max(1,t/\tau)]^{r},\quad \eta_n(t)=\eta_n^{(0)}[\max(1,t/\tau)]^{r}, \end{equation}

where $r$ is an exponent, $\tau$ is a decay time scale and $\nu _n^{(0)}$ and $\eta _n^{(0)}$ denote the coefficients at early times. The motivation for time-dependent $\nu _n(t)$ and $\eta _n(t)$ is twofold. On the one hand, perfect self-similarity can only be expected if the value of $r$ is suitably adjusted (Yousef, Haugen & Brandenburg 2004). On the other hand, negative values of $r$ are convenient from a numerical point of view because they allow us to consider Lundquist numbers that gradually increase as the energy of the turbulence at the highest wavenumbers decays.

3.2. Parameters and diagnostic quantities

In the following, we normalize the magnetic energy spectrum $E_{M}(k,t)\equiv \text {Sp}\left ( {\boldsymbol {B}} \right ){/}2\mu _0$ such that $\int \text {d}k\,E_{M}(k,t)=\left \langle {\boldsymbol {B}}^2/2\mu _0 \right \rangle \equiv {\mathcal {E}}_{M}$, where ${\mathcal {E}}_{M}$ denotes the mean magnetic energy density. The magnetic integral scale $\xi _{M}$ has been defined in (2.2), where $\xi _{M}^{-1}(t)$ corresponds approximately to the location where the spectrum peaks, and therefore $\xi _{M}^{-1}(0)\approx k _\text {peak}$. We define the generalized Lundquist number for hyperdiffusion as

(3.5)\begin{equation} {Lu}_n(t)=\frac{v_{A}^\text{rms}\xi_{M}^{2n-1}}{\eta_n}, \end{equation}

where $v_{A}$ is the Alfvén velocity, and $v_A^{\rm rms}$ is its root-mean-square value.

To characterize the decay of ${\mathcal {E}}_{M}$ and $I_H$, and the increase of $\xi _{M}$, we define the instantaneous scaling coefficients

(3.6a,b)\begin{equation} p(t)={-}\frac{\text{d}\ln{\mathcal{E}}_{M}}{\text{d}\ln t},\quad p_H(t)={-}\frac{\text{d}\ln I_H}{\text{d}\ln t},\quad q(t)=\frac{\text{d}\ln\xi_{M}}{\text{d}\ln t}. \end{equation}

Parametric representations of $p(t)$ vs $q(t)$ are useful in distinguishing different decay behaviours (Brandenburg & Kahniashvili 2017). For purely hydrodynamic turbulence, for example, $p(t)$ and $q(t)$ tend to evolve along a line $p/q\approx 5$ toward a point where $p=10/7$ and $q=2/7$ (Proudman & Reid 1954), while for fully helical MHD turbulence, one sees an evolution along $p/q\approx 1$ toward $p=q=2/3$ (Hatori 1984). In those cases, the spectrum evolves in an approximately self-similar manner of the form

(3.7)\begin{equation} E(k,t)=\xi(t)^{-\beta}\phi\left(k\xi(t)\right), \end{equation}

where $\phi (\kappa )$ is a universal function, and $\beta$ is an exponent that describes the gradual decline of the height of the peak. Self-similarity is not a stringent requirement, and perfect self-similarity can also not be expected in a numerical simulation owing to the limited range of scales that can be resolved. As shown by Olesen (1997), the ideal MHD equations are invariant under rescaling; see also Appendix B.2. This causes additional constraints that will be discussed below. To understand the expectations following from self-similarity and invariance under rescaling, we recall the basic relations involving $p$ and $q$ in Appendix B.

3.3. Role of Alfvén and reconnection times

As emphasized above, the ratio $p/q$ is determined by the conserved invariant; see Brandenburg & Kahniashvili (2017) and Appendix B. If $B^2\xi _{M}^{1+\beta }\sim {\mathcal {E}}_{M}\xi _{M}^{1+\beta }$ remains constant during the decay, then ${\mathcal {E}}_{M}\xi _{M}^{1+\beta }\propto t^{-p+q(1+\beta )}\propto t^0$, which gives

(3.8)\begin{equation} p=(1+\beta)q. \end{equation}

In particular, $\beta =0$ when magnetic helicity is conserved, and $\beta =3/2$ when the Hosking invariant is conserved.⁶

The decay time scale yields a second relation between $p$ and $q$. When the decay time is the Alfvén time, we have $t_\text {decay}\sim \xi _{M}/v_{A}^\text {rms}\sim \xi _{M}/{\mathcal {E}}_{M}^{1/2}$,⁷ and therefore

(3.9)\begin{equation} 1=q+p/2. \end{equation}

In that case, together with (3.8), we get $\beta =2/q-3$, which is the relation expected from the invariance under rescaling the MHD equations; see Appendix B.2. Alternatively, if MHD turbulence decay is controlled by slow (Sweet–Parker) reconnection, which occurs for ${Lu}_n^{1/2n}\lesssim 10^4$ (Loureiro et al. 2005; Loureiro, Schekochihin & Cowley 2007), we have $t_\text {decay}\sim {Lu}_n^{1/2n}\xi _{M}/{\mathcal {E}}_{M}^{1/2}$, which gives $1=-p/4n+(1-1/2n)q-r/2n+q+p/2$, and therefore

(3.10)\begin{equation} (8n-2)q+(2n-1)p=4n+2r, \end{equation}

where we have also taken into account the time dependence of the hyperresistivity, $\eta _n\propto t^r$, so that ${Lu}_n\propto {\mathcal {E}}_{M}^{1/2}\xi _{M}^{2n-1}t^{-r}\propto t^{-p/2+(2n-1)q-r}$.⁸ Equation (3.8) together with (3.9) or (3.10) uniquely determine the values of $p$ and $q$ in terms of $\beta$ or $\beta$, $r$, and $n$, respectively. For non-helical decaying MHD turbulence, which is proposed by Hosking & Schekochihin (2021) to conserve $I_H$, and therefore $\beta =3/2$, we have $p=10/9$ and $q=4/9$ with an Alfvén time scale that is independent of $n$ and $r$. With the reconnection time scale we have instead

(3.11a,b)\begin{equation} p=\frac{10(2n+r)}{26n-9},\quad q=\frac{4(2n+r)}{26n-9}. \end{equation}

3.4. Initial conditions

As initial condition we use a magnetic field with a given energy spectrum proportional to $k^4$ for $k< k_\text {peak}$ and proportional to $k^{-5/3}$ for $k>k_\text {peak}$, where $k_\text {peak}$ denotes the position of the peak of the initial spectrum. We assume random phases, which make the field Gaussian distributed. The smallest wavenumber of the domain is ${k}_1$, and $k_\text {peak}$ is taken to be 60 or 200; see table 1 for a summary of the simulations. The corresponding data files for these runs can be found in the online material; see Zhou, Sharma & Brandenburg (2022) for the published data sets used to compute each of the figures of the present paper.

Table 1. Summary of runs, where $N^3$ is the resolution, ${\mathcal {E}}_{M}(0)$ is the initial magnetic energy density, and the last column lists the hyper-Lundquist numbers at the beginning and the end of the simulations. Runs ending with ‘c’ use a constant (hyper)diffusivity, whereas those with ‘t’ use a time-dependent value $\propto t^{-3/7}$.

We use the publicly available Pencil Code (Pencil Code Collaboration et al. 2021). By default, it uses sixth-order accurate finite differences and the third-order Runge–Kutta timestepping scheme of Williamson (1980). The different methods for calculating $\mathcal {I}_H(R)$ have been implemented and are publicly available since the revision of 20 May 2022. By default, the code computes ${\boldsymbol {A}}$ in the resistive gauge, but we have also implemented the calculation of ${\boldsymbol {A}}_{C}={\boldsymbol {A}}-{\boldsymbol {\nabla }}\varLambda$ in the Coulomb gauge by solving $\nabla ^2\varLambda ={\boldsymbol {\nabla }}\boldsymbol {\cdot }{\boldsymbol {A}}$ for the gauge potential $\varLambda$.

4. Results

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 2021).

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. (2011) 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. 2015) and may support the reconnection-controlled decay picture (Bhat et al. 2021; Zhou et al. 2021); see also Zhou et al. (2019) for two-dimensional MHD, and Zhou et al. (2020) 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. 2015; Brandenburg & Kahniashvili 2017; Reppin & Banerjee 2017; Bhat et al. 2021; Zhou et al. 2021), 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. 2005, 2007).

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.⁹ 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 (2021), 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 1958). 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 (2020).

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 2020). 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 )$,¹⁰ 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 (2017). 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. (2017) 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)$.

5. Conclusion

The present results have verified that the Hosking integral is conserved in the limit of large Lundquist numbers. This implies that $I_H$ can indeed control the decay behaviour of MHD turbulence. On dimensional grounds, one would expect $\beta =3/2$, i.e. the solution evolves in a $pq$ diagram along a line where $p=(1+\beta )q=5q/2$. Our highest resolution simulations with $2048^3$ mesh points show that this is compatible with this expectation.

We have also shown that different methods of determining $I_H$ all lead to the same result. The preferred method is based on the magnetic helicity density spectrum, $\text {Sp}\left ( h \right )$, which is also the simplest method. Furthermore, by comparing the resistive and Coulomb gauges, we find that $I_H$ is indeed gauge invariant to high precision.

Whether or not the decay time is governed by the reconnection time rather than the Alfvén time remains uncertain, although our comparison between runs with constant and time-dependent ordinary magnetic diffusivities, i.e. $\eta _1=\text {const.}$ and $\eta _1=\eta _1(t)$, suggest that the reconnection time scale might indeed be the relevant one. Clearly, high resolution simulations are required to obtain meaningful scaling results. A resolution of $2048^3$ is just beginning to yield conclusive results, but higher resolution would be desired to address the role of reconnection more conclusively.