Scaling of the Hosking integral in decaying magnetically-dominated turbulence

Hongzhe Zhou^{1,2}†, Ramkishor Sharma^{1,3} and Axel Brandenburg^{1,3,4,5}

1Nordita, KTH Royal Institute of Technology and Stockholm University, Hannes Alfvéns väg

12, SE-10691 Stockholm, Sweden

2Tsung-Dao Lee Institute, Shanghai Jiao Tong University, 800 Dongchuan Road, Shanghai

200240, People’s Republic of China 3The Oskar Klein Centre, Department of Astronomy, Stockholm University, AlbaNova, SE-10691 Stockholm, Sweden

4McWilliams Center for Cosmology & Department of Physics, Carnegie Mellon University,

Pittsburgh, PA 15213, USA

5School of Natural Sciences and Medicine, Ilia State University, 3-5 Cholokashvili Avenue, 0194 Tbilisi, Georgia

(21 November 2022)

The Saffman helicity invariant of Hosking and Schekochihin (2021, PRX 11, 041005), which we here call the Hosking integral, has emerged as an important quantity that may govern the decay properties of magnetically dominated nonhelical turbulence. Using a range of different computational methods, we confirm that this quantity is indeed gaugeinvariant 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 et al. 1999), including superfluid (Nore et al. 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 et al. 2001; Campanelli 2007; Kahniashvili et al. 2010; Brandenburg et al. 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 et al. 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

† Email address for correspondence: hongzhe.zhou@su.se

2 H. Zhou, R. Sharma, and A. Brandenburg

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 et al. 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, 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 behavior 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 behavior 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 nonvanishing 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(x, t) = A·B, is invariant in the ideal (non-resistive) limit and also independent of the gauge A → A′ = A − ∇Λ for any scalar Λ. Here, B = ∇ × A is the magnetic field expressed in terms of the magnetic vector potential 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 nonhelical 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 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 twopoint correlation function of momentum, which is proportional to the velocity 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, 2020), and Bhat et al. (2021). Hosking & Schekochihin (2021) have shown that the two physical models, Alfvén- vs. reconnection-controlled decay, can

† An inverse cascade implies a local transfer in wavenumber space to smaller k. If the transfer is nonlocal, or if locality in k space is uncertain, one rather speaks just of inverse transfer.

Scaling of the Hosking integral 3

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, but in the absence of net magnetic helicity, using a range of different methods and 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 section 2, discuss then our numerical stimulation setup in section 3, and present our results in section 4. We conclude in section 5.

2. Expressions for the Hosking integral

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

2.1. The box-counting methodDefinition 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(x, t) = A · B, which we call I_{H}(R). The latter, which we call I_{H}(R), is given by (Hosking & Schekochihin 2021)

I_{H}(R) = ∫_{V_{R}} d^{3}r ⟨h(x)h(x + r)⟩ .  (2.1)

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

ξ_{M} = \frac{∫ dk k^{−1}E_{M}(k)}{∫ dk E_{M}(k)},  (2.2)

where E_{M}(k) is the magnetic energy spectrum. An alternative length scale can be computed from the spectrum of the magnetic helicity density, ξ_{h}. A comparison between ξ_{M} and ξ_{h} is made in appendix A, where we show that, although the ratio ξ_{M}/ξ_{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 ξ_{M}, we have

I_{H}(R) ≃ ∫_{V_{R}} d^{3}r ⟨h(x)h(x)⟩ ∝ R^{3},  (2.3)

whereas when R ≫ ξ_{M}, but at the same time much smaller than the system scale L, equation (2.1) can be written in an equivalent form by assuming that the volume average approximates the ensemble averageI_{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

† We are grateful to Keith Moffatt for alerting us to the fact that the formerly used term “Saffman helicity invariant” may be misleading, because the term helicity invariant is reserved for integrals which are chiral in character. Moreover, Saffman never considered helicity in his papers. The term “magnetic helicity density correlation integral” may be more appropriate but rather clumsy. Following Schekochihin (2020), we now refer to it as the Hosking integral. We use the term integral instead of invariant as long as we are not in the ideal limit.

‡ We thank the anonymous referee for suggesting this.

I_{H} = \frac{1}{V_{R}} ⟨H_{V_{R}}^{2}⟩ , H_{V_{R}} = ∫_{V_{R}} d^{3}r h(r).  (2.4)

Within this asymptotic range, 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 ⟨H_{V_{R}}^{2}⟩ ∝ V_{R}^{2} (V_{R}/ξ_{M}^{3})^{−1} ∝ V_{R} (Hosking & Schekochihin 2021).

2.2. The box-counting method

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

I_{H}(R) = \frac{1}{V V_{R}} ∫_{V} d^{3}x [∫_{V_{R}} d^{3}r h(x + r)]^{2} ,  (2.5)

where V is the volume of the simulation box. We call this the box-counting (BC) method. Upon a Fourier transformation,† equation (2.5) can be recast as a weighted integral

I_{H}(R) = \frac{1}{V} ∫ \frac{d^{3}k}{(2π)^{3}} w_{R}(k) h^{∗}(k)h(k).  (2.6)

The weight function w_{R}(k) depends on the shape of V_{R}. For a cubic region with length 2R, we have

w_{R}(k) = w_{cube}^{BC}(k) ≡ 8R^{3} ∏_{i=1}^{3} j_{0}^{2}(k_{i}R),  (2.7)

whereas for a spherical region with radius R,

w_{R}(k) = w_{sph}^{BC}(k) ≡ \frac{4πR^{3}}{3} [\frac{6j_{1}(kR)}{kR}]^{2} .  (2.8)

Here j_{0}(x) = sin x/x and j_{1}(x) = (sin x − x cos x)/x^{2} are the first and second order spherical Bessel functions, respectively. Note also that in equation (2.8), w_{sph}^{BC}(k) depends only on k = |k|. This allows us to rewrite equation (2.6) as just a one-dimensional integral,

I_{H}(R) = ∫_{0}^{∞} dk w_{sph}^{BC}(k) Sp (h) ,  (2.9)

where

Sp (h) = \frac{1}{V} \frac{k^{2}}{(2π)^{3}} ∫_{|k|=k} dΩ_{k} h̃^{∗}(k)h̃(k)  (2.10)

is the Fourier spectrum of the magnetic helicity density and Ω_{k} is the solid angle in Fourier space, normalized such that ∫ dk Sp (h) = ⟨h^{2}⟩.‡ The magnetic and kinetic energy spectra can then be written analogously as E_{M}(k) = Sp (B) /2μ_{0} and E_{K}(k) = ρ_{0} Sp (u) /2, where μ_{0} is the magnetic permeability and ρ_{0} is the mean density.

† We define the Fourier transform of f(x) as f̃(k) = ∫ d^{3}x f(x) e^{−ik·x} and the inverse transformations as f(x) = ∫ (2π)^{−3}d^{3}k f̃(k)e^{ik·x}.

‡ Note that Sp (h) is same as Θ in equation (32) of Hosking & Schekochihin (2021).

Scaling of the Hosking integral 5

2.3. The correlation integral method

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

⟨h(x)h(x + r)⟩ = \frac{1}{V} ∫_{V} d^{3}x h(x)h(x + r) = \frac{1}{V} ∫ \frac{d^{3}k}{(2π)^{3}} h̃^{∗}(k)h̃(k)e^{ik·r},  (2.11)

which we call the correlation-integral (CI) method. The Hosking integralI_{H}(R) can then again be recast in the form of equation (2.6), but the weight functions are now slightly different. For a cubic region V_{R} we have

w_{R}(k) = w_{cube}^{CI}(k) ≡ 8R^{3} ∏_{i=1}^{3} j_{0}(k_{i}R),  (2.12)

and for a spherical region we have

w_{R}(k) = w_{sph}^{CI}(k) ≡ \frac{4πR^{3}}{3} \frac{6j_{1}(kR)}{kR}.  (2.13)

Note that w_{sph}^{CI}(k) can also be used in equation (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 ≪ ξ_{M}^{−1} (Hosking & Schekochihin 2021),

Sp (h) = \frac{I_{H}}{2π^{2}} k^{2} + O(k^{4}),  (2.14)

given that h(k) is statistically isotropic, and therefore

I_{H} = lim_{k→0} \frac{2π^{2}}{k^{2}} Sp (h) .  (2.15)

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 = ρc_{s}^{2}, where ρ is the density. We allow for the possibility of hyperviscous and hyperdiffusive dissipation of kinetic and magnetic energies and solve the equations for A in the resistive gauge. The full set of equations is

\frac{D ln ρ}{Dt} = −∇ · u,  (3.1)

\frac{Du}{Dt} = −c_{s}^{2}∇ ln ρ + \frac{1}{ρ} [J × B + ∇ · (2ρν_{n}∇^{2(n−1)}S)] ,  (3.2)
\frac{∂A}{∂t} = u × B + η_{n}∇^{2n}A,  (3.3)

where J = ∇ × B/μ_{0} is the current density, ν_{n} is the viscosity, η_{n} is the magnetic diffusivity, S_{ij} = (∂_{i}u_{j} + ∂_{j}u_{i})/2 − δ_{ij}∇ · u/3 are the components of the rate-of-strain tensor S, and n denotes the degree of hyperviscosity or hyperdiffusivity with n = 1 corresponding to ordinary viscous diffusive operators in equations (3.2) and (3.3), respectively. Equation (3.3) is here formulated in the resistive gauge, i.e., the scalar 6 H. Zhou, R. Sharma, and A. Brandenburg

potential is φ = −η_{n}∇^{2(n−1)}∇·A in the uncurled induction equation ∂A/∂t = −E−∇φ, where E = −u × B + η_{n}∇^{2(n−1)}μ_{0}J is the electric field.

In some cases, we assume ν_{n} and η_{n} to be time-dependent, which is indicated by writing ν_{n}(t) and η_{n}(t), respectively. In those cases, we assume a power law variation (for t > τ),

ν_{n}(t) = ν_{n}^{(0)}[max(1, t/τ)]^{r}, η_{n}(t) = η_{n}^{(0)}[max(1, t/τ)]^{r},  (3.4)

where r is an exponent, τ is a decay time scale, and ν_{n}^{(0)} and η_{n}^{(0)} denote the coefficients at early times. The motivation for time-dependent ν_{n}(t) and η_{n}(t) is two-fold. On the one hand, perfect self-similarity can only be expected if r = (1 − β)/(3 + β) is obeyedthe value of r is suitably adjusted (Yousef et al. 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) ≡ Sp (B) /2μ_{0} such that ∫ dk E_{M}(k, t) = ⟨B^{2}/2μ_{0}⟩ ≡ E_{M}, where E_{M} denotes the mean magnetic energy density. The magnetic integral scale ξ_{M} has been defined in equation (2.2), where ξ_{M}^{−1}(t) corresponds approximately to the location where the spectrum peaks, and therefore ξ_{M}^{−1}(0) ≈ k_{peak}. We define the generalized Lundquist number for hyperdiffusion as

Lu_{n}(t) = \frac{v_{A}^{rms}ξ_{M}^{2n−1}}{η_{n}} ,  (3.5)

where v_{A} is the Alfvén velocity.

To characterize the decay of E_{M} and I_{H}, and the increase of ξ_{M}, we define the instantaneous scaling coefficients

p(t) = −\frac{d ln E_{M}}{d ln t} , p_{H}(t) = −\frac{d ln I_{H}}{d ln t} , q(t) = \frac{d ln ξ_{M}}{d ln t} .  (3.6)

Parametric representations of p(t) vs q(t) are useful in distinguishing different decay behaviors (Brandenburg & Kahniashvili 2017). For purely hydrodynamic turbulence, for example, p(t) and q(t) tend to evolve along a line p/q ≈ 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 ≈ 1 toward p = q = 2/3 (Hatori 1984). In those cases, the spectrum evolves in an approximately self-similar manner of the form

E(k, t) = ξ(t)^{−β}ϕ(kξ(t)),  (3.7)

where ϕ(κ) is a universal function, and β is an exponent that describes the gradual decline of the height of the peak. Self-similarity is not a stringent requirement, and perfect selfsimilarity 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}ξ_{M}^{1+β} ∼ E_{M}ξ_{M}^{1+β} remains constant during the decay, then E_{M}ξ_{M}^{1+β} ∝ t^{−p+q(1+β)} ∝ t^{0}, which gives

p = (1 + β)q.  (3.8)

In particular, β = 0 when magnetic helicity is conserved, and β = 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_{decay} ∼ ξ_{M}/v_{A}^{rms} ∼ ξ_{M}/E_{M}^{1/2},‡ and therefore

1 = q + p/2.  (3.9)

In that case, together with equation (3.8), we get β = 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} ≲ 10^{4} (Loureiro et al. 2005, 2007), we have t_{decay} ∼ Lu_{n}^{1/2n}ξ_{M}/E_{M}^{1/2}, which gives 1 = −p/4n + (1 − 1/2n)q − r/2n + q + p/2, and therefore

(8n − 2)q + (2n − 1)p = 4n + 2r,  (3.10)

where we have also taken into account the time dependence of the hyperresistivity, η_{n} ∝ t^{r}, so that Lu_{n} ∝ E_{M}^{1/2}ξ_{M}^{2n−1}t^{−r} ∝ t^{−p/2+(2n−1)q−r}.¶ Equation (3.8) together with equations (3.9) or (3.10) uniquely determine the values of p and q in terms of β or β, r, and n, respectively. For non-helical decaying MHD turbulence, which is proposed by Hosking & Schekochihin (2021) to conserve I_{H}, and therefore β = 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

p = \frac{10 (2n + r)}{26n − 9}, q = \frac{4 (2n + r)}{26n − 9}.  (3.11)

3.4. Initial conditions

As initial condition we use a magnetic field with a given energy spectrum proportional to k^{4} for k < k_{peak} and proportional to k^{−5/3} for k > k_{peak}, where k_{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_{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 et al. (2022) for the published data sets used to compute each of the figures of the present paper.

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 I_{H}(R) have been implemented and are publicly available since the revision of May 20, 2022. By default, the code computes A in the resistive gauge, but we have also implemented the calculation of A_{C} = A − ∇Λ in the Coulomb gauge by solving ∇^{2}Λ = ∇ · A for the gauge potential Λ.

† Hosking & Schekochihin (2021) introduced the exponent α and wrote the conserved quantity as B^{α}ξ_{M}. Then, α = 2/(1 + β), which is 2 and 4/5 for β = 0 and 3/2, respectively.

‡ This also implies that dE_{M}/dt ∼ −E_{M}/t_{decay} ∼ −E_{M}^{3/2}/ξ_{M}, i.e., the energy dissipation is, as expected, proportional to (v_{A}^{rms})^{3}/ξ_{M}.

¶ Note that for n → ∞, equation (3.10) does not yield equation (3.9). Mathematically, this is because n also enters in the expression for Lu_{n} itself. To see this, one can write instead Lu_{n}^{1/2m} to obtain (4m+4n−2)q +(2m−1)p = 4m+2r, which recovers equation (3.10) in the physically meaningful case m = n, and equation (3.9) when m → ∞ with n < ∞ enforced.

Run  N^{3} E_{M}(0) k_{peak} ν_{n} = η_{n} n r  Lu_{n}
K200D3t 1024^{3} 4 × 10^{−3} 200 1 × 10^{−14} 3 −3/7 1.2 × 10 → 8.5 × 10^{4}
K200D3c 1024^{3} 4 × 10^{−3} 200 1 × 10^{−14} 3 0  1.2 × 10 → 1.4 × 10^{3}
K60D1c 1024^{3} 7 × 10^{−3} 60 5 × 10^{−6} 1 0  5.3 × 10 → 2.5 × 10^{2}
K60D1bt 2048^{3} 3 × 10^{−1} 60 2 × 10^{−6} 1 −3/7 2.1 × 10^{3} → 7.5 × 10^{3}
K60D1bc 2048^{3} 3 × 10^{−1} 60 2 × 10^{−6} 1 0  1.3 × 10^{3} → 4.0 × 10^{3}
K60D3t 1024^{3} 7 × 10^{−3} 60 1 × 10^{−14} 3 −3/7 1.5 × 10^{3} → 6.3 × 10^{6}
K60D3bt 2048^{3} 3 × 10^{−1} 60 4 × 10^{−16} 3 −3/7 1.6 × 10^{5} → 3.2 × 10^{7}
K60D3bc 2048^{3} 3 × 10^{−1} 60 4 × 10^{−16} 3 0  1.7 × 10^{5} → 2.1 × 10^{8}

Table 1. Summary of runs, where N^{3} is the resolution, E_{M}(0) is the initial magnetic energymagnetic energy density, and Lu_{n} is the hyper-Lundquist number at the beginning and the end of the simulation. Runs ending with ‘c’ use a constant (hyper)diffusivity, whereas those with ‘t’ use a time-dependent value ∝ t^{−3/7}.

4. Results

4.1. I_{H}(R) from different methods

In Figure 1 we compare for run K60D1c the temporal dependence of I_{H}(R) and I_{H} computed from the methods introduced in section 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 I_{H}(R) and I_{H} curves have different magnitudes, as expected from the fact that w_{R}(k) are different; see figure 1(a)figures 1(a) and (b). However, the scaling properties, i.e., the values of p_{H} = −d ln I_{H}/d ln t are unchanged among all the methods; see figure 1(bc). In what follows, we use the CI method with cubic regions throughout.

4.2. Gauge invariance

In figure 2(a) we compare the R-dependence of the Hosking integralcorrelation integral 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 ≪ ξ_{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).

The evolution of the auto-correlation C_{h}(R) = ⟨h(x)h(x + R)⟩ can be computed from (4πR^{2})^{−1}dI_{H}(R)/dR; see figure 2(b). Using a time-dependent normalization of the abscissa, we see that h is always correlated at roughly 2ξ_{M}(t), as also evident from the inset.

Let us point out at this point that Sp (h) 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 φ = u · A, can lead to large gradients in the evolution of 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) ≡ Sp (B) /2μ_{0} and magnetic helicity density spectra Sp (h) for run K60D1bt at different times. From the energy

Figure 1. For run K60D1c, comparing I_{H} in different methods, and the decay exponents p_{H} = −d ln I_{H}/d ln t.Comparing results for run K60D1c in different methods. (a) 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} = −d ln I_{H}/d ln t.

Figure 2. Results for run K60D1c. (a) Comparing I_{H}(R) from different gauges. (b) The auto-correlation curves C_{h}(R). The inset shows 4πR^{2}C_{h}(R). Note that the abscissa is normalized by ξ_{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.

spectra, we see that the inertial range shifts both in k and in amplitude , but preserves the same spectral slopes at large and, more or less, also at small k.. 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 twodimensional 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} ≈ 10. This is presumably a finite size effect, and so the late time evolution may not be reliable unless the initial k_{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_{peak} and on the numerical resolution, we must regard our results with care and should discuss the possibility of artifacts as we reach the

Figure 3. Results for run K60D1bt. (a) Magnetic energy spectrum. (b) Spectrum of magnetic helicity density. (c) The non-dimensionalized and compensated Sp (h); see equation (4.11). (d) Time evolution of the non-dimensionalized I_{H}, equation (4.12)equations (4.12) and (4.13). For the first three panels, the vertical gray lines mark the asymptotic scale chosen to be k = 2π/(2R) = 4.35. with R = 0.115L, at which value we have computed Ĩ_{H} in (d).

limits of what can be considered safe. Thus, we can conclude that, except for very late times, 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_{peak}. At small k, however, Sp (h) has a random noise spectrum ∝ k^{2}, which remains unchanged also at late times. This agrees with our expectations; see equation (2.14), which allows us to determine I_{H} from the spectral value at small k.

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 2R/L = 0.23R = 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 the asymptotic value 10/9an 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}.

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 artifact 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 ∼ 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

Figure 4. (a) Time evolutions of the normalized I_{H}. (b) The instantaneous decay exponents of I_{H} (p_{H}) and 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 intetral (p_{C}). (d) Time evolution of I_{H}(R) for run K60D3bc. The vertical line indicates the asymptotic scale chosen to be R/(2π) = 0.115.

time evolution of the Hoskingcorrelation integral I_{H}(R). At small R ≲ ξ_{M}, we see an R^{3} scaling. Thus, ⟨H_{V_{R}}^{2}⟩ /V_{R} is proportional to V_{R} ∝ R^{3}, and therefore ⟨H_{V_{R}}^{2}⟩ 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, ξ_{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

I_{c}(R) = ∫_{V_{R}} d^{3}r ⟨h_{c}(x)h_{c}(x + r)⟩ ,  (4.1)

where h_{c} = u · 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 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} = −d ln I_{C}/d ln t. Hence, I_{C} is not as well conserved as I_{H}.

† For this reason, partially or fully helical magnetic fields would not result in a meaningful definition of I_{H}. In those cases, the conservation of magnetic helicity determines β.

12 H. Zhou, R. Sharma, and A. Brandenburg

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

⟨h̃^{∗}(k_{1})h̃(k_{2})⟩ = ∫ \frac{d^{3}k′}{(2π)^{3}} \frac{d^{3}k′′}{(2π)^{3}} ϵ_{ijk}ϵ_{abc}k′_{k}k′′_{c} ⟨Ã_{i}^{∗}(k_{1} − k′)Ã_{j}^{∗}(k′)Ã_{a}(k_{2} − k′′)Ã_{b}(k′′)⟩ .

(4.2)

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

⟨Ã_{i}^{∗}(k_{1})Ã_{j}(k_{2})⟩ = (2π)^{3}δ^{3}(k_{1} − k_{2})P_{ij}(k_{1})M_{A}(k_{1}),  (4.3)

where P_{ij}(k) = δ_{ij} − k_{i}k_{j}/k^{2}, and M_{A} is related to the magnetic energy spectrum via M_{A}(k) = 2π^{2}E_{M}(k)/k^{4}. Taking advantage of the fact that Ã_{i}(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 k_{1}, k_{2}, k_{3}, and k_{4},

⟨Ã_{i}^{∗}(k_{1})Ã_{j}^{∗}(k_{2})Ã_{a}(k_{3})Ã_{b}(k_{4})⟩
= ⟨Ã_{i}^{∗}(k_{1})Ã_{j}(−k_{2})⟩ ⟨Ã_{a}^{∗}(−k_{3})Ã_{b}(k_{4})⟩ + ⟨Ã_{i}^{∗}(k_{1})Ã_{a}(k_{3})⟩ ⟨Ã_{j}^{∗}(k_{2})Ã_{b}(k_{4})⟩
+ ⟨Ã_{i}^{∗}(k_{1})Ã_{b}(k_{4})⟩ ⟨Ã_{j}^{∗}(k_{2})Ã_{a}(k_{3})⟩
=(2π)^{6} [δ^{3}(k_{1} + k_{2})δ^{3}(k_{3} + k_{4})P_{ij}(k_{1})P_{ab}(k_{3})M_{A}(k_{1})M_{A}(k_{3})
+ δ^{3}(k_{1} − k_{3})δ^{3}(k_{2} − k_{4})P_{ia}(k_{1})P_{jb}(k_{2})M_{A}(k_{1})M_{A}(k_{2})
+δ^{3}(k_{1} − k_{4})δ^{3}(k_{2} − k_{3})P_{ib}(k_{1})P_{ja}(k_{2})M_{A}(k_{1})M_{A}(k_{2})] .  (4.4)

Together with equation (2.10), we find

Sp (h) = \frac{k^{2}}{2} ∫_{0}^{∞} dk′ ∫_{−1}^{1} dα \frac{k^{′2} + μ^{2}k^{′2} − 2μk′|k − k′|}{k^{′2}|k − k′|^{4}}E_{M}(|k − k′|)E_{M}(k′), (4.5)

where α = cos(k, k′), and

|k − k′| = \sqrt{k^{2} + k^{′2} − 2αkk′}, μ = \frac{αk′k − k^{′2}}{k′\sqrt{k^{2} + k^{′2} − 2αkk′}} .  (4.6)

Note that equation (4.5) is similar to the expression for Sp (B^{2}); see equation (27) of Brandenburg & Boldyrev (2020).

Consider a piecewise power-law spectrum,

E_{M}(k) = \{ \begin{matrix}E_{peak} (k/k_{peak})^{4} , & k ⩽ k_{kpeak} \\ E_{peak} (k/k_{peak})^{−s} , & k > k_{kpeak} \end{matrix}\} ,  (4.7)

we have the relationswhere 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

E_{peak} = \frac{20}{17}E_{M}ξ_{M}, k_{peak} = \frac{1}{2ξ_{M}} .  (4.8)

The leading-order term (∝ k^{2}) of equation (4.5) can be readily obtained by taking k = 0 in the integrand, which gives

Sp (h)^{early} = \frac{136 E_{peak}^{2}}{95 k_{peak}^{3}} k^{2} + O(k^{3}) = \frac{5120E_{M}^{2}ξ_{M}^{5}}{323} k^{2} + O(k^{3}).  (4.9)

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

For s = 2 we obtain

Sp (h)^{late} = \frac{131072E_{M}^{2}ξ_{M}^{5}}{13125} k^{2} + O(k^{3}).  (4.10)

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

S̃p(h) = \frac{323}{5120 E_{M}^{2}ξ_{M}^{5}k^{2}} Sp (h) ,  (4.11)

which is indeed ∼ O(1) initially at small k, but then increases to O(10) at later times. Had we used s = 2, S̃p(h) only obtains larger values because the numerical factor in equation (4.10) is smaller than that in equation (4.9). This could be caused byIndeed, the increase in S̃p(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 + ∑_{i=1}^{3} ⟨B_{i}^{4}⟩ / ⟨B_{i}^{2}⟩^{2} /3, to be ∼ −0.24 at the end of run K60D1bt. This is slightly larger in magnitude compared with earlier work (Brandenburg & Boldyrev 2020). Thus,

we conclude that, for the presented runs, the Lundquist number is not yet sufficiently

high and the influence from non-Gaussianity is subdominant.The influence from non-Gaussianity can be explicitly checked by computing the right-hand side of equation (4.5) and comparing with Sp (h),† 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 equation (4.5), which approximates Sp (h) using E_{M}, under-estimates the true value of Sp (h) by nearly one order of magnitude, although sharing a similar spectral shape with the latter. Hence, the non-dimensionalized S̃p(h) would be larger than unity also by approximately one order of magnitude, in agreement with figure 3(c).

Furthermore, we cannot rule out the possibility that the initial non-self-similar evolution of the magnetic field also contributes to the increase of S̃p(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 ξ_{M}, and our estimation of S̃p(h) will be lower at early time.

The magnitude of I_{H} can be estimated from equation (2.15), and a proper normaliza-

† We thank David Hosking for pointing this out.

14 H. Zhou, R. Sharma, and A. Brandenburg

tion is

Ĩ_{H} = \frac{323}{10240π^{2} E_{M}^{2}ξ_{M}^{5}} I_{H} ≃ 0.0032 E_{M}^{−2}ξ_{M}^{−5}I_{H}  (4.12)

which is plotted in figure 3(d)with s = 5/3 for early times, and

Ĩ_{H} = \frac{13125}{262144π^{2} E_{M}^{2}ξ_{M}^{5}} I_{H} ≃ 0.0051 E_{M}^{−2}ξ_{M}^{−5}I_{H}  (4.13)

with s = 2 for late times. Both curves are plotted in figure 3(d). The increasing value again highlights the non-idealnessnon-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 section 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.

The most reliable runs are those in figures 6(c) and (d), where N^{3} = 2048^{3} mesh points have been used. We clearly see that the solution evolves along the β = 3/2 line, as expected when the decay is governed by the Hosking integral. For a constant value of η_{1}, figure 6(d), the solution also reaches the line p = 2(1−q), which is referred to as the scaleinvariance 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 η_{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 β = 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 β = 0 line toward the point (p, q) = (2/3, 2/3). They also considered the nonhelical case in their figure 2(b), where the system evolved along the β = 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 β = 2 line toward (p, q) = (6/5, 2/5). Both results were arguably still consistent with an evolution along β = 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 behavior of MHD turbulence. On dimensional grounds, one would expect β = 3/2, i.e., the solution evolves in a pq diagram along a line where p = (1 + β)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, Sp (h),

Figure 6. Panels (a) to (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 equations (3.8), (3.9), and (3.10), respectively.

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., η_{1} = const and η_{1} = η_{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.

We thank David Hosking, Nuno Loureiro, Kandaswamy Subramanian, and the referees for helpful comments and discussions. We also thank Keith Moffatt and Alex Schekochihin for discussions regarding the problematic naming of the Hosking integral as a Saffman helicity invariant. This work was supported by the Swedish Research Council (Vetenskapsrådet, 2019-04234). Nordita is sponsored by Nordforsk. 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 and Linköping.

Appendix A. Different definitions of the length scale

In section 2.1, we used ξ_{M} to characterize the length scale of magnetic fluctuations.

Figure 7. The ratio between the energy integral scale ξ_{M} and the helicity density integral scale ξ_{h}.

An alternative definition can be

ξ_{h} ≡ \frac{∫ k^{−1}Sp (h) dk}{∫ Sp (h) dk} .  (A 1)

This is plotted here in figure 7 for all the runs. Although the ratio is time-dependent, note that this is on a log-linear scale and thus is a rather weak dependence in comparison with the power-law decay of the magnetic energy.

Appendix B. Self-similarity and invariance under rescaling

In sections 3.2 and 3.3, we discussed pq diagrams. In this appendix, we summarize the essential relations between p, q, and the exponent β, which describes the gradual decline of the spectral peak as it moves toward smaller k.

B.1. Self-similarity

A self-similar spectrum implies that its shape is given by a universal function ϕ = ϕ(κ) such that the spectrum is of the form

E(k, t) = ξ(t)^{−β}ϕ(kξ(t)),  (B 1)

where ϕ(κ) depends on just one argument κ ≡ kξ(t), and ξ = ξ(t) is the temporal dependence of the integral scale, which can be used to describe the gradual shift of the peak of the spectrum towards smaller k.

B.2. Invariance under rescaling

Describing the time dependence of ξ through a power law ξ ∝ t^{q} means that a rescaling of length, x → x′ = xℓ implies a rescaling of time, t → t′ = tℓ^{1/q}. Since E(k, t) has dimensions

[E(k, t)] = [x]^{3}[t]^{−2},  (B 2)

and since ϕ(kξ(t)) must not change under rescaling, we have

E(kℓ^{−1}, tℓ^{1/q}) = ℓ^{3−2/q+β} [ξ(t)ℓ]^{−β} ϕ(kξ),  (B 3)

and therefore β = 2/q − 3 or q = 2/(β + 3); see Olesen (1997).

Scaling of the Hosking integral 17

Conserved quantity dimension q β = 2/q − 3 p = (β + 1)q
1  −1  0 (=0/4)
⟨A · B⟩  [x]^{3}[t]^{−2} 2/3  0  2/3 (=4/6)
⟨A^{2}⟩ (?)  [x]^{4}[t]^{−2} 1/2  1  1 (=8/8)
I_{H}  [x]^{9}[t]^{−4} 4/9  3/2  10/9
I_{u}  [x]^{5}[t]^{−2} 2/5  2  6/5 (=12/10)
1/3  3  4/3 (=16/12)
Loitsiansky  [x]^{7}[t]^{−2} 2/7  4  10/7(=20/14)

Table 2. Summary of coefficients. The question mark on ⟨A^{2}⟩ indicates that the significance of this quantity is questionable.

B.3. Relation to p

The exponent p quantifies the temporal scaling of E_{M} = ∫ dk E_{M}(k, t) ∝ t^{−p}. Inserting (B 1), we have

t^{−p} ∝ E_{M} = ξ_{M}(t)^{−(β+1)} ∫ d(kξ_{M}) E_{M}(kξ_{M}) ∝ ξ_{M}(t)^{−(β+1)} ∝ t^{−(β+1)q}  (B 4)

and therefore p = (β + 1)q. Using invariance under rescaling, we have

p = 2(1 − q).  (B 5)

Table 2 lists various candidates of conserved quantities and the corresponding values of q, β, and p. In a pq diagram, this represents the line on which the instantaneous scaling coefficients p(t) and q(t) tend to settle; see Brandenburg & Kahniashvili (2017), where we called this the self-similarity line. However, if the reconnection time scale really becomes the dominant one, it may be more meaningful to call p = 2(1 − q) the scale-invariance line. It agrees with the assumption of the relevant time scale being the Alfvén time; see equation (3.9). On the other hand, if there really are two distinct scales that evolve differently, the result cannot be self-similar; see also section 11.2.3 of Schekochihin (2020) for a discussion.

REFERENCES

Bhat, Pallavi, Zhou, Muni & Loureiro, Nuno F. 2021 Inverse energy transfer in decaying, three-dimensional, non-helical magnetic turbulence due to magnetic reconnection. MNRAS 501 (2), 3074–3087, arXiv: 2007.07325.

Brandenburg, Axel & Boldyrev, Stanislav 2020 The Turbulent Stress Spectrum in the

Inertial and Subinertial Ranges. ApJ 892 (2), 80, arXiv: 1912.07499. Brandenburg, Axel, Enqvist, Kari & Olesen, Poul 1996 Large-scale magnetic fields

from hydromagnetic turbulence in the very early universe. Phys. Rev. D 54 (2), 1291– 1300, arXiv: astro-ph/9602031.

Brandenburg, Axel & Kahniashvili, Tina 2017 Classes of Hydrodynamic and Magnetohydrodynamic Turbulent Decay. Phys. Rev. Lett. 118 (5), 055102, arXiv: 1607.01360.

Brandenburg, Axel, Kahniashvili, Tina, Mandal, Sayan, Pol, Alberto Roper, Tevzadze, Alexander G. & Vachaspati, Tanmay 2017 Evolution of hydromagnetic turbulence from the electroweak phase transition. Phys. Rev. D 96 (12), 123528, arXiv: 1711.03804.

18 H. Zhou, R. Sharma, and A. Brandenburg

Brandenburg, A., Kahniashvili, T. & Tevzadze, A. G. 2015 Nonhelical Inverse Transfer

of a Decaying Turbulent Magnetic Field. PhRvL 114, 075001, arXiv: 1404.2238. Campanelli, Leonardo 2007 Evolution of Magnetic Fields in Freely Decaying

Magnetohydrodynamic Turbulence. Phys. Rev. Lett. 98 (25), 251302, arXiv: 0705.2308. Candelaresi, Simon, Hubbard, Alexander, Brandenburg, Axel & Mitra,

Dhrubaditya 2011 Magnetic helicity transport in the advective gauge family. Physics of Plasmas 18 (1), 012903–012903, arXiv: 1010.6177.

Christensson, Mattias, Hindmarsh, Mark & Brandenburg, Axel 2001 Inverse cascade in decaying three-dimensional magnetohydrodynamic turbulence. Phys. Rev. E 64 (5), 056405, arXiv: astro-ph/0011321.

Davidson, P. A. 2000 Was Loitsyansky correct? A review of the arguments. Journal of

Turbulence 1 (1), 6. Durrer, Ruth & Neronov, Andrii 2013 Cosmological magnetic fields: their generation,

evolution and observation. Astron. Astrophys. Rev. 21, 62, arXiv: 1303.7121. Frisch, U., Pouquet, A., Leorat, J. & Mazure, A. 1975 Possibility of an inverse cascade

of magnetic helicity in magnetohydrodynamic turbulence. Journal of Fluid Mechanics 68, 769–778.

Hatori, Tadatsugu 1984 Kolmogorov-Style Argument for the Decaying Homogeneous MHD

Turbulence. JPSJ 53 (8), 2539. Hosking, David N. & Schekochihin, Alexander A. 2021 Reconnection-Controlled Decay of

Magnetohydrodynamic Turbulence and the Role of Invariants. Physical Review X 11 (4), 041005, arXiv: 2012.01393.

Hosking, David N. & Schekochihin, Alexander A. 2022 Cosmic-void observations reconciled with primordial magnetogenesis. arXiv e-prints p. arXiv:2203.03573, arXiv: 2203.03573.

Kahniashvili, Tina, Brandenburg, Axel, Tevzadze, Alexander G. & Ratra, Bharat 2010 Numerical simulations of the decay of primordial magnetic turbulence. Phys. Rev. D 81 (12), 123002, arXiv: 1004.3084.

Kitsionas, S., Federrath, C., Klessen, R. S., Schmidt, W., Price, D. J., Dursi, L. J., Gritschneder, M., Walch, S., Piontek, R., Kim, J., Jappsen, A. K., Ciecielag, P. & Mac Low, M. M. 2009 Algorithmic comparisons of decaying, isothermal, supersonic turbulence. A&A 508 (1), 541–560, arXiv: 0810.4599.

Krause, F. & Rüdiger, G. 1975 On the Turbulent Decay of Strong Magnetic Fields and the

Development of Sunspot Areas. Solar Physics 42 (1), 107–119. Loureiro, N. F., Cowley, S. C., Dorland, W. D., Haines, M. G. & Schekochihin,

A. A. 2005 X-Point Collapse and Saturation in the Nonlinear Tearing Mode Reconnection. Phys. Rev. Lett. 95 (23), 235003, arXiv: physics/0507206.

Loureiro, N. F., Schekochihin, A. A. & Cowley, S. C. 2007 Instability of current sheets and formation of plasmoid chains. Physics of Plasmas 14 (10), 100703–100703, arXiv: astro-ph/0703631.

Mac Low, Mordecai-Mark, Klessen, Ralf S., Burkert, Andreas & Smith, Michael D. 1998 Kinetic Energy Decay Rates of Supersonic and Super-Alfvénic Turbulence in Star-Forming Clouds. Phys. Rev. Lett. 80 (13), 2754–2757, arXiv: astroph/9712013.

Neronov, Andrii & Vovk, Ievgen 2010 Evidence for Strong Extragalactic Magnetic Fields

from Fermi Observations of TeV Blazars. Science 328 (5974), 73, arXiv: 1006.3504. Nore, C., Abid, M. & Brachet, M. E. 1997 Decaying Kolmogorov turbulence in a model

of superflow. Physics of Fluids 9 (9), 2644–2669. Olesen, P. 1997 Inverse cascades and primordial magnetic fields. Physics Letters B 398, 321–

325, arXiv: astro-ph/9610154. Pencil Code Collaboration, Brandenburg, Axel, Johansen, Anders, Bourdin,

Philippe, Dobler, Wolfgang, Lyra, Wladimir, Rheinhardt, Matthias, Bingert, Sven, Haugen, Nils, Mee, Antony, Gent, Frederick, Babkovskaia, Natalia, Yang, Chao-Chin, Heinemann, Tobias, Dintrans, Boris, Mitra, Dhrubaditya, Candelaresi, Simon, Warnecke, Jörn, Käpylä, Petri, Schreiber, Andreas, Chatterjee, Piyali, Käpylä, Maarit, Li, Xiang-Yu, Krüger, Jonas, Aarnes, Jørgen, Sarson, Graeme, Oishi, Jeffrey, Schober, Jennifer, Plasson,

Scaling of the Hosking integral 19

Raphaël, Sandin, Christer, Karchniwy, Ewa, Rodrigues, Luiz, Hubbard, Alexander, Guerrero, Gustavo, Snodin, Andrew, Losada, Illa, Pekkilä, Johannes & Qian, Chengeng 2021 The Pencil Code, a modular MPI code for partial differential equations and particles: multipurpose and multiuser-maintained. The Journal of Open Source Software 6 (58), 2807, arXiv: 2009.08231.

Pouquet, A., Frisch, U. & Leorat, J. 1976 Strong MHD helical turbulence and the

nonlinear dynamo effect. Journal of Fluid Mechanics 77, 321–354. Proudman, I. & Reid, W. H. 1954 On the Decay of a Normally Distributed and Homogeneous

Turbulent Velocity Field. Philosophical Transactions of the Royal Society of London Series A 247 (926), 163–189.

Reppin, Johannes & Banerjee, Robi 2017 Nonhelical turbulence and the inverse transfer of

energy: A parameter study. Phys. Rev. E 96 (5), 053105, arXiv: 1708.07717. Saffman, P. G. 1967 The large-scale structure of homogeneous turbulence. Journal of Fluid

Mechanics 27, 581–593. Schekochihin, Alexander A. 2020 MHD Turbulence: A Biased Review. arXiv e-prints p.

arXiv:2010.00699, arXiv: 2010.00699. Stalp, Steven R., Skrbek, L. & Donnelly, Russell J. 1999 Decay of Grid Turbulence in

a Finite Channel. Phys. Rev. Lett. 82 (24), 4831–4834. Subramanian, Kandaswamy 2016 The origin, evolution and signatures of primordial magnetic

fields. RPPh 79 (7), 076901. Vachaspati, Tanmay 2021 Progress on cosmological magnetic fields. Reports on Progress in

Physics 84 (7), 074901, arXiv: 2010.10525. Warhaft, Z. & Lumley, J. L. 1978 An experimental study of the decay of temperature

fluctuations in grid-generated turbulence. Journal of Fluid Mechanics 88, 659–684. Williamson, J. H. 1980 Low-Storage Runge-Kutta Schemes. Journal of Computational Physics

35 (1), 48–56. Woltjer, L. 1958 On Hydromagnetic Equilibrium. Proceedings of the National Academy of

Science 44 (9), 833–841. Yousef, Tarek A., Haugen, Nils Erland L. & Brandenburg, Axel 2004 Self-similar

scaling in decaying numerical turbulence. Phys. Rev. E 69 (5), 056303, arXiv: astroph/0312505.

Zhou, H., Sharma, R. & Brandenburg, A. 2022 Datasets for Scaling of the Saffman helicity integral in decaying magnetically-dominated turbulence, doi:10.5281/zenodo.6640511

(v2022.06.14); see also http://www.nordita.org/~brandenb/projects/Saffman/ for

easier access . Zhou, Muni, Bhat, Pallavi, Loureiro, Nuno F. & Uzdensky, Dmitri A. 2019 Magnetic

island merger as a mechanism for inverse magnetic energy transfer. Physical Review Research 1 (1), 012004, arXiv: 1901.02448.

Zhou, Muni, Loureiro, Nuno F. & Uzdensky, Dmitri A. 2020 Multi-scale dynamics of magnetic flux tubes and inverse magnetic energy transfer. Journal of Plasma Physics 86 (4), 535860401, arXiv: 2001.07291.

Zhou, Muni, Wu, David H., Loureiro, Nuno F. & Uzdensky, Dmitri A. 2021 Statistical description of coalescing magnetic islands via magnetic reconnection. Journal of Plasma Physics 87 (6), 905870620, arXiv: 2104.13757.

Zrake, Jonathan 2014 Inverse Cascade of Nonhelical Magnetic Turbulence in a Relativistic Fluid. ApJL 794 (2), L26, arXiv: 1407.5626.