Reconnection-controlled decay of magnetohydrodynamic turbulence and the role of invariants

We present a new theoretical picture of magnetically dominated, decaying turbulence in the absence of a mean magnetic field. We demonstrate that such turbulence is governed by the reconnection of magnetic structures, and not by ideal dynamics, as has previously been assumed. We obtain predictions for the magnetic-energy decay laws by proposing that turbulence decays on reconnection timescales, while respecting the conservation of certain integral invariants representing topological constraints satisfied by the reconnecting magnetic field. As is well known, the magnetic helicity is such an invariant for initially helical field configurations, but does not constrain non-helical decay, where the volume-averaged magnetic-helicity density vanishes. For such a decay, we propose a new integral invariant, analogous to the Loitsyansky and Saffman invariants of hydrodynamic turbulence, that expresses the conservation of the random (scaling as $\mathrm{volume}^{1/2}$) magnetic helicity contained in any sufficiently large volume. Our treatment leads to novel predictions for the magnetic-energy decay laws: in particular, while we expect the canonical $t^{-2/3}$ power law for helical turbulence when reconnection is fast (i.e., plasmoid-dominated or stochastic), we find a shallower $t^{-4/7}$ decay in the slow `Sweet-Parker' reconnection regime, in better agreement with existing numerical simulations. For non-helical fields, for which there currently exists no definitive theory, we predict power laws of $t^{-10/9}$ and $t^{-20/17}$ in the fast- and slow-reconnection regimes, respectively. We formulate a general principle of decay of turbulent systems subject to conservation of Saffman-like invariants, and propose how it may be applied to MHD turbulence with a strong mean magnetic field and to isotropic MHD turbulence with initial equipartition between the magnetic and kinetic energies.


I. INTRODUCTION
The nature of the decay of magnetohydrodynamic (MHD) turbulence is an important outstanding problem in fluid dynamics, with far-reaching consequences in astrophysics, from the evolution of primordial magnetic fields in cosmology [1][2][3] to the dynamics of the solar wind [4]. Naturally, the subject of decaying turbulence is one with a long history. In the hydrodynamic case, the basic problem of determining the exponent of the energy-decay power law was solved by Kolmogorov, in the third of his seminal 1941 papers on turbulence [5]. Kolmogorov's approach can be summarised as follows: (i ) identify an ideal invariant that is better conserved than the kinetic energy, then (ii ) posit a decay of the kinetic energy, occurring on the dynamical timescale, that conserves that invariant. In the case of hydrodynamic turbulence, Kolmogorov identified the Loitsyansky integral, as the relevant invariant. Physically, the conservation of I L represents the net conservation of angular momen- * david.hosking@physics.ox.ac.uk tum of the turbulent eddies [6,7]. Note, though, that the invariant controlling the decay is not simply the mean angular momentum, L = r × u , which, while conserved, is zero by isotropy. The Loitsyansky integral is therefore an example of an invariant that encodes the conservation of a quantity that individual eddies are expected to possess, but that has vanishing mean value due to its randomly directed nature. Eq. (1) implies the scaling U 2 L 5 ∼ const, where L is the correlation scale of the turbulence, and U is the typical velocity at that scale. Together with the identification of the dynamical timescale as τ ∼ L/U , this is enough to fix the decay rate of the kinetic energy E = U 2 /2, as which results in Kolmogorov's famous decay law This result has been confirmed to excellent precision numerically [7,8].
In this paper, we show how Kolmogorov's philosophy can be adapted to MHD turbulence. Despite extensive studies of the self-similar decaying solutions admitted by the MHD equations [9][10][11], two problems have hindered past attempts to predict specific decay laws. First, there appear to exist a number of different regimes, depending on properties of the initial conditions [11][12][13][14]. Restricting attention, for the moment, to the magnetically dominated case, where magnetic energy is much greater than kinetic, there are two canonical possibilities. First, there are helical field configurations, where the volumeaveraged magnetic-helicity density, h = A·B , is nonzero. Then, conservation of h (a phenomenon sometimes referred to as 'selective decay' [14][15][16][17][18][19][20]) provides a scaling that can be used to constrain the decay laws, B 2 L ∼ const [21][22][23], where B is the typical size of the magnetic field at the correlation scale L. However, magnetic helicity is not sign-definite, so there also exist nonhelical field configurations, for which h B 2 L. For such fields, the conservation of h does not impose a constraint on the decay. In this work, we show that even in such cases, the decay is still controlled by helicity conservation, in a manner formally analogous to the control of the decay of hydrodynamic turbulence by angular-momentum conservation.
Second, since, besides velocity, MHD has an additional field, B, there is no longer a dimensional inevitability in the identification of the decay timescale, as there was in hydrodynamics. Previous treatments [21][22][23][24] have assumed the Alfvénic scaling U ∼ B in order to determine the ideal timescale uniquely [25], though this is not well reproduced in numerics, where B U appears to be maintained if it was true initially, and, furthermore, for helical magnetic fields, a faster decay of the kinetic energy than the magnetic energy is often observed [1,22,[26][27][28][29][30].
In fact, it is intuitively clear that relaxation on ideal timescales may not be possible for a strong initial magnetic field, because of the topological constraints imposed by magnetic-flux freezing. As was hypothesised by J.B. Taylor, magnetic fields with non-trivial topologies relax via the reconnection of magnetic field lines [15,20], which transfers magnetic energy to larger scales. Reconnection, therefore, provides a physical explanation for the inverse transfer of magnetic energy observed in both helical and non-helical decaying MHD turbulence [26,[31][32][33][34][35][36][37][38]. However, unless the values of the dissipation coefficients are sufficiently small for reconnection to occur in the plasmoid-dominated [39] or stochastic [40] regime, magnetic reconnection occurs in the so-called 'Sweet-Parker' regime [41,42], and is slow, i.e., has a rate that is proportional to a negative fractional power of the Lundquist number, S = BL/η, where η is the fluid resistivity. It should then be expected that the decay will proceed on the Sweet-Parker reconnection timescale, not the ideal one. The critical Lundquist number at which magnetic reconnection becomes fast (i.e., independent of S) is very large, ∼ 10 4 [39], much larger than the typical Reynolds numbers at which Eq. (3) becomes a good description of the decay of hydrodynamic turbulence (Re 10 2 [8]). The requirement of such a large Lundquist number (corresponding to very thin current sheets), together with the large scale separation between the box size and the energy-containing scale needed to eliminate finite-box-size effects, mean that direct numerical simulations aiming to measure decay laws will generally be in the slow-reconnection regime [43]. In [44], it was demonstrated that two-dimensional MHD turbulence indeed decays on the Sweet-Parker timescale, and similar evidence has been presented for three-dimensional, nonhelical turbulence too [38], though it was interpreted as arising from the two-dimensional decay mechanism put forward by [44]. One of the main goals of the present work will be to verify the reconnection-controlled nature of the three-dimensional decay in both the helical and non-helical cases, and to establish the corresponding decay laws for both energies.
Sweet-Parker reconnection is defined by the conditions of (i ) efficient conversion of magnetic energy to kinetic energy of reconnection outflows, and (ii ) a balance between the inductive term, ∇ × (u × B), and the resistive dissipation term in the MHD induction equation, so that reconnection occurs in a time-invariant manner. This last requirement means that reconnection-controlled decaying turbulence in the Sweet-Parker regime is very different from decaying hydrodynamic turbulence, in that it is sensitive to the precise form of the dissipation term. Importantly, this means that different decay power laws are expected in numerical simulations depending on whether Laplacian dissipation, ∝ η∇ 2 B ≡ η 2 ∇ 2 B, or hyperdissipation, ∝ η n ∇ n B, is employed. This fact has not been widely appreciated, and leads to a simple test of whether reconnection indeed governs the decay timescale: simulations at moderately large (but not so large so as to be in the 'fast' reconnection regime) Lundquist numbers should exhibit different decay power laws depending on the order of hyper-dissipation, n. In this paper, we present such simulations, conducted with the incompressible, spectral MHD code Snoopy [45]. They turn out to be in excellent agreement with these expectations.
An outline of the rest of this paper is as follows. In Section II, we consider the decay of helical MHD turbulence from a magnetically dominated state. Since magnetic helicity is better conserved than magnetic energy in the limit of vanishing resistivity, helicity is precisely a quantity that can control the evolution in the same way as the Loitsyansky integral does in hydrodynamics. Applying a Kolmogorov-style argument to helicity conservation for a decay occuring on the ideal timescale, L/B, yields a power law of t −2/3 for both the magnetic and kinetic energies [21][22][23], though this decay law is not well-supported by numerics, where a shallower decay law for magnetic energy, and a steeper decay law for kinetic energy are typically observed [1,22,[26][27][28][29][30]46]. We show that both of these results are expected for a decay occurring via magnetic reconnection in the Sweet-Parker regime. In particular, we show that the magnetic energy should decay as t −4/7 , with t −2/3 only achieved by fast reconnection. We also find that, provided the dominant flows are contained within Sweet-Parker sheets, the faster decay of kinetic energy, as t −5/7 , is a natural consequence of their changing aspect ratio.
In Section III, we consider the decay of non-helical MHD turbulence from a magnetically dominated state. The mechanism controlling this type of decay has so far remained unknown: because the mean helicity density vanishes, its conservation cannot be used to derive a scaling relation relating B, the characteristic magnetic field, to its characteristic length scale L. Numerically, decay laws for both the magnetic and kinetic energies of close to t −1 have been observed [30,[33][34][35]38], though there is no definitive theoretical explanation for this behaviour. An influential idea is that in the absence of an integral invariant, the decay might satisfy the well-known scaling symmetry of the MHD equations [47], including dissipative terms. Such a decay would have a t −1 power law for the magnetic energy [9,10,48,49]. Another suggestion is that the non-helical decay is effectively two-dimensional. In this case, it is the conservation of 'anastrophy', or the square of the magnetic vector potential, that should control the decay [14,17]. A Kolmogorov-style argument then leads to a t −1 power law independently of whether the decay occurs on the ideal [21,24,30,34], or the Sweet-Parker [37,38,44] timescale (see Appendix A). That both treatments predict the same power law is a coincidence related to the fact that anastrophy conservation implies constant Lundquist number for n = 2 resistive dissipation [44] (incidentally, the scaling argument is essentially this same point with the direction of implication reversed). However, it is not clear why fully three-dimensional, isotropic turbulence should twodimensionalise in this way, nor why any special significance should be given to constant Lundquist number. Indeed, we show in Section III that both are inconsistent with numerical evidence.
Instead, we propose a treatment of the non-helical decay controlled by the conservation of fluctuations in magnetic helicity. The key point is that the vanishing of the total magnetic helicity does not necessarily imply that any given magnetic field structure is non-helical. Indeed, non-helical magnetic structures generally relax on ideal timescales to zero magnetic energy, assuming they are not constrained by higher-order topological invariants. We therefore expect that the natural state of the turbulence will be to contain a collection of helical structures, though there will be equal abundances of positive-and negativehelicity structures so that the zero-overall-helicity constraint is satisfied. For such a turbulence, we identify a new integral invariant whose relation to magnetic helicity is precisely analogous to that of the Loitsyansky integral to angular momentum. We refer to this invariant as the 'Saffman helicity invariant', due to the close analogy between it and the integral invariant proposed by Saffman for hydrodynamic turbulence [50], and also between our arguments and the arguments usually associated with the Saffman invariant and its various anisotropic generalisations [7,51]. As we will show, the conservation of our new invariant implies the scaling B 4 L 5 ∼ const. This scaling implies a magnetic-energy-decay power law of t −20/17 if reconnection occurs on the Sweet-Parker timescale, or t −10/9 if reconnection is fast. These power laws are different from t −1 , but are still in excellent agreement with published numerical results, and with our own numerical results presented below. We also find that, for non-helical magnetic fields, the rate at which the aspect ratio of the current sheets changes is much smaller than in the helical case, explaining why a faster decay of kinetic energy is not observed for a non-helical decay from initial states that have small kinetic energy.
Finally, in Section IV, we discuss the behaviour of systems with fractional helicity, which we show will ultimately always transition to the fully helical regime as long as the system size is sufficiently large. We also discuss possible applications of the Saffman formalism to wider classes of turbulent decays. As an example, we suggest the existence of a Saffman-type cross-helicity invariant that may control the critically balanced decay of MHD turbulence in the presence of a mean magnetic field, recently studied in [37]. We also conjecture that the simultaneous conservation of both the Saffman-type cross-helicity invariant and the magnetic helicity might govern the initial period of decay of an MHD state starting with U ∼ B.

II. DECAY OF HELICAL TURBULENCE
We first consider the decay of MHD turbulence from an initial state where magnetic energy dominates kinetic (B U ), and the magnetic field is helical, i.e., the volume-averaged magnetic-helicity density, is large, ∼ B 2 L. The case of initial parity between the two energies (U ∼ B) will be discussed in Section IV E. Magnetic helicity is a topological invariant -the total magnetic helicity of a collection of flux tubes is equal to the amount of (signed) flux linked by these tubes. As such, its conservation is related to Alfvén's theorem, which states that for η n = 0, the magnetic field is frozen into fluid motions, ensuring that all topological invariants are precisely conserved. The special significance of magnetic helicity is that, unlike other topological invariants, it remains approximately conserved (i.e., better conserved than energy) for small but non-vanishing η n [52]. This statement is true independently of the reconnection regime. A proof, adapted from [53], is as follows.
A. Helicity is conserved by (hyper-resistive) reconnection In hyper-dissipative MHD, the evolution of the magnetic helicity in a closed volume whose surface is everywhere normal to B satisfies Note that setting n = 2 here yields the familiar expression in terms of the current helicity. After integrating the integral on the right-hand side by parts, applying the Cauchy-Schwarz inequality, and then integrating by parts again, one obtains where dE M /dt = η n d 3 r B · (∇ n B) is the rate of magnetic-energy decay due to Ohmic heating. For n = 2, the other integral in Eq. (6) is just twice the magnetic energy [53]. More generally, we can write where δ η is the resistive dissipation scale. Eq. (6) then implies Eq. (8) states that the rate of change of magnetic helicity is smaller than the rate of the energy decay due to Ohmic heating (which will be even smaller than the true magnetic-energy-decay rate, because magnetic energy can also be converted to the kinetic energy of reconnection outflows) by a factor equal to the ratio of the integral scale to the resistive dissipation length scale, which becomes arbitrarily small as η n → 0 + . Q.E.D. It may appear counter-intuitive that reconnection, a process that, by definition, changes the topology of magnetic field lines, can conserve helicity, a topological invariant. The resolution of this apparent paradox is that self-linkages, i.e., twists of the magnetic flux tube are also associated with helicity. For example, during the unlinking of two linked tori by reconnection to form a single torus, the resulting torus ends up twisted, and the total helicity of the configuration is conserved [52].

B. Theory of helical decay
The conservation of the volume-averaged magnetic helicity, Eq. (4), implies the scaling The remaining ingredient required to compute the magnetic-energy decay law is the decay timescale, as a function of U , B and L. The simplest possible treatment is to assume Alfvénic dynamics, with U ∼ B, similarity between the integral scales of the magnetic and kinetic energies, L, and, therefore, a decay timescale L/B ∼ L/U . For future reference, we compute the expected decay power law under such an assumption for a scaling more general than Eq. (9), viz., Eq. (2) becomes d dt with solution Now setting α = 2 for helicity conservation, the canonical t −2/3 power law is recovered [21]. However, this prediction has proved to be in poor agreement with numerics, which have found a decay of the magnetic energy closer to t −1/2 , and, unexpectedly, a decay of the kinetic energy faster than t −2/3 [1,22,[26][27][28][29][30]46]. These discrepancies are readily resolved by assuming the decay to occur on the reconnection, rather than ideal, timescale. Naturally, the reconnection timescale depends on the reconnection regime, i.e., on whether the reconnection is slow (Sweet-Parker) or fast (plasmoid-dominated [39] or stochastic [40]). Fast reconnection, by definition, occurs on dynamical timescales, so will again produce a t −2/3 decay. However, as explained in the Introduction, extant numerical simulations of decaying MHD turbulence mostly probe the slow regime, owing to the large Lundquist numbers and hence large resolutions required for fast reconnection to take place. In this case, the reconnection timescale is where S n = BL n−1 /η n is the hyper-Lundquist number. Using Eq. (13) for the decay timescale, Eq. (11) becomes d dt the solution of which is Again, for a helicity-conserving decay, we set α = 2 to obtain a power law of For n = 2 (Laplacian dissipation), the power law is which is indeed shallower than p M = 2/3. Indeed, some recent studies at large resolution have reported p M 0.58, in remarkable agreement with this prediction [30,46]. However, we caution against direct comparison with those simulations, because they employed time-dependent dissipation coefficients. For numerical studies using n = 4 hyper-dissipation [22], Eq. (15) predicts an even slower magnetic-energy-decay exponent of p M = 8/17 0.47. This is in excellent agreement with the p M 0.5 found numerically in [22].
It may appear counter-intuitive that reconnection can result in a faster decay of kinetic energy than magnetic energy, as reconnection outflows are typically Alfvénic, i.e., the outflow velocity is approximately equal to the upstream Alfvén speed of the magnetic field prior to reconnection -a condition that is hard-wired into the Sweet-Parker scalings. However, the current sheets, where reconnection occurs, are not volume-filling. Denoting the current sheet width by δ, the volume occupied by the current sheet formed when two structures of volume L 3 reconnect and merge is L 2 δ. Therefore, we expect for Alfvénic outflows, where E K and E M should be understood as the total kinetic and magnetic energies in the system, respectively. This result allows for the possibility of different decay rates for the kinetic and magnetic energy, because the ratio δ/L need not be constant in time. For example, with hyper-dissipative Sweet-Parker Thus, kinetic energy is indeed expected to decay more quickly than magnetic energy, simply due to the changing aspect ratio of the Sweet-Parker sheets. For n = 2, this effect is relatively modest: the kinetic-energy decay exponent is 5/4 times greater than the magnetic one, so the decay exponents for E M and E K are p M = 4/7 0.57 and p K = 5/7 0.71, respectively. For n = 4, though, the kinetic-energy decay exponent is 13/8 times greater than the magnetic one, so these decay exponents become p M = 8/17 0.47 and p K = 13/17 0. 76. We note that the validity of Eq. (18) requires the initial kinetic energy to be much smaller than the magnetic energy, by a factor δ/L, as otherwise the kinetic energy contained in outflows may be subdominant to the rest. This condition will be true for our simulations, which are initialised with U = 0, but may not be true in all situations of interest. However, as we shall discuss in Section IV E, there is reason to believe that decaying helical MHD turbulence should be driven towards Eq. (18), even if E K is larger than (δ/L)E M initially.

C. Numerical results
In this section, we compare the theoretical predictions of the previous section with numerical simulations of MHD turbulence, initialised as a Gaussian random magnetic field with characteristic scale 1/33 of the box size. These simulations employ the incompressible, spectral MHD code Snoopy [45], with Prandtl number Pm ≡ ν n /η n = 1. We provide further information about Snoopy and describe the details of our numerical setup in Appendix D.
In simulations with n = 2 dissipation, resolution constraints prevent the use of Lundquist numbers that are sufficiently large to achieve good conservation of the magnetic helicity. Nonetheless, supposing that helicity decays as For small but non-vanishing η n , therefore, we expect to find an α somewhat smaller than 2 such that B α L ∼ const. We can determine this value of α numerically by measuring E M (t) and where E M (k) is the spectral magnetic-energy density.
As long as η n is not too large, the decay should still occur on the Sweet-Parker timescale. We can then use Eq. (15) to determine the expected magnetic-energy decay exponent based on the empirically determined value of α, and compare it with the value measured in our simulations. While we do not expect that this procedure should yield exact agreement between the predicted and empirical decay exponents, as it neglects the role of Ohmic diffusion (which, when S n is small, will ultimately become more important to the decay than reconnection), we expect approximate agreement that becomes better as S n increases and α, therefore, becomes closer to 2.
In Fig. 1, we present the results of such a comparison, plotting the empirically measured values of p M against those of α, for simulations with n = 2 and n = 4 (how the error bars are determined is explained in Appendix D). This figure shows remarkable agreement between our simulations and the Sweet-Parker decay curves (coloured), despite the fact that we do not reach the asymptotic value of α = 2 with n = 2 dissipation. For n = 4, we also find excellent agreement, and do reach α = 2. As other authors have noted [22], this asymptotic scaling is reached much more rapidly in the hyper-dissipative case, and indeed we were forced to choose relatively small values of S 1/n n in order to populate the part of Fig. 1 with α < 2. To illustrate this point further, we plot two simulations with n = 6 hyper-dissipation, which also exhibit the α = 2 scaling (and are almost coincident in Fig. 1). The faster attainment of the correct asymptotic scaling with hyper-dissipation will be important in establishing the correct value of α for non-helical turbulence, previously unknown, in the next section.
Interestingly, we note that for α 2, the decay exponent p M is consistently somewhat larger than our theoretical prediction based on Sweet-Parker reconnection. This suggests these simulations may be at the start of the transition to the fast-reconnection regime. We do not find the same transition in two-dimensional simulations (see Appendix A), despite employing even larger Lundquist numbers, which is consistent with the intuitive expectation that fast reconnection should 'turn on' more quickly in three dimensions, due to turbulence in the reconnection region.  15), for decays occurring on the Sweet-Parker timescale, with n = 2, n = 4 and n = 6 shown in blue, red, and magenta, respectively. The grey solid curve depicts the 'ideal' scaling given by (12). Simulation results are in excellent agreement with the coloured curves, and not with the grey curve. This confirms that the decay takes place on the Sweet-Parker, rather than ideal, timescale. The full set of decay curves from which this plot was derived can be found in Appendix D.
For the decay of the kinetic energy also, we find that the predictions of the previous section are in good agreement with our simulations, as shown in Fig. 2. The upper panel shows the evolution of the magnetic-and kineticenergy spectra for a run with n = 4, η 4 = 2 × 10 −8 , confirming that at any given time, kinetic energy is contained at much smaller scales than the magnetic energy, consistent with the expectation that reconnection outflows should have a width δ L. The inset shows the relative sizes of the total magnetic and kinetic energies in the same run, which are in excellent agreement with Eq. (18). The lower panel shows the kinetic-energy decay exponents, p K , plotted against the corresponding exponents for the magnetic energy, p M , confirming that a faster decay of kinetic energy is realised in our simulations, in reasonably good agreement with our theoretical prediction, Eq. (19).
We note that, while we do observe faster decay of kinetic energy than magnetic in our simulations, the difference is not as stark as in the numerical study by [22], who found E K ∝ E 2 M . As we will discuss in Section IV, this discrepancy may arise because the initial state in their study was one with equipartition between the magnetic and kinetic energies, U ∼ B, unlike the U B we have employed here. We will argue that when U ∼ B, the conservation of cross-helicity, even though the latter is , where time is in code units based on normalising the box size and initial mean-square magnetic field to 2π and 1, respectively, so that 1 time unit is approximately the initial Alfvén crossing time of the box. The peak of the magnetic energy is at much larger scales (smaller k) than the peak of the kinetic energy, consistent with the expectation that the kinetic energy should be contained within the Sweet-Parker sheets. Inset: decay of the total magnetic and kinetic energies. The total kinetic energy curve is much below the magnetic energy curve, and coincides with (δ/L)EM (black), in agreement with Eq. (18), with δ/L computed as the ratio of the wavenumbers at which kEM (k) and kEK (k) peak. Bottom panel: Plot of the kinetic-energy decay exponent, pM , against the magnetic-energy decay exponent, pK , as measured in simulations with S 1/n n, 0 > 9.0. Results are in reasonable agreement with the theoretical prediction, Eq. (19) (coloured lines), and are inconsistent with EK ∝ EM (grey line). a sign-indefinite quantity, might play an important role in governing the decay. Simultaneously conserving magnetic helicity and fluctuations in the cross-helicity, under the same formalism as we are about to propose for helicity conservation in non-helical turbulence, does imply the scaling E K ∝ E 2 M conjectured in [22]. To conclude, the results of this section represent what we consider compelling evidence that helical, magnetically dominated MHD turbulence relaxes by the selfsimilar coalescence of magnetic structures via magnetic reconnection, as was suggested by J.B. Taylor [20]. Physically, this implies that the correct way to think about the system is as consisting of a collection of magnetic structures that are unable to relax under ideal dynamics due to the topological constraints imposed by the flux freezing, and, therefore, relax via coalescence on a timescale at which these constraints can be broken, i.e., on the reconnection timescale.

III. DECAY OF NON-HELICAL TURBULENCE
Having established a theory of the reconnectioncontrolled decay of helical turbulence from a magneti-cally dominated state, we now consider the case of nonhelical turbulence, i.e., turbulence for which the volumeaveraged magnetic-helicity density [Eq. (4)] vanishes. Here, as in Section II, we consider the decay from an initial state with predominant magnetic energy (B U ), postponing the discussion of the case with U ∼ B until Section IV E. As we have already noted, the mechanisms controlling the evolution of such turbulence are not well understood. Numerically, a power law close to t −1 has been measured [30,[33][34][35]38], prompting comparisons with the two-dimensional decay [21,24,44], which conserves 'anastrophy', or the square of the magnetic vector potential (see Appendix A), resulting in a t −1 decay law independently of the reconnection regime. The evidence that has been presented for this picture in three dimensions relies on demonstrating that the meansquare magnetic vector potential (defined according to some particular, necessarily non-unique gauge choice) Cartoon of a typical merger of two helical structures ('blobs'). As explained in the main text, B 2 L = B 2 L for helical turbulence in which this is the only allowed process. (b) This additional 'annihilation' process is also possible in non-helical turbulence, and should occur equally frequently. The presence of this process modifies the previous scaling to changes more slowly with time than does the magnetic energy [38]. However, this will be true for any decay satisfying B α L ∼ const for any α > 0. Here, we propose a different theory of non-helical decay.

A. Qualitative theory of non-helical decay
The key point, already made in Section I, is that a vanishing mean helicity density does not imply that individual magnetic structures are non-helical, because helicity is not a sign-definite quantity. Moreover, if not constrained by higher-order topological invariants, nonhelical magnetic structures will relax to zero energy on the ideal timescale (as required by J.B. Taylor relaxation [15,20]). For example, consider a toroidal structure without a twist -such a structure will shrink under the magnetic tension force (driving outflows along its axis) to zero magnetic energy. In contrast, for a toroidal structure with a net twist relative to the poloidal axis, such a relaxation is topologically impossible. Thus, we expect a collection of helical structures of different signs to be the natural state of non-helical MHD turbulence. Indeed, visualisations of our simulations do appear to support this intuition, see Fig. 3.
To motivate the more formal approach to follow, we first present an informal 'cartoon' of the expected dy-namics ( Fig. 4). Consider a volume filled with helical magnetic structures, which we shall refer to as 'blobs', as we wish to remain agnostic about their precise morphology. As in the helical case, we expect that topological constraints will hinder their relaxation on ideal timescales, and that instead the blobs will evolve via coalescences with other blobs on the prevailing reconnection timescale.
For simplicity, suppose that all blobs have helicity ∼ H, and for the moment, that they all have the same sign of helicity. When any two blobs merge, the resulting structure will have helicity H = 2H, implying that the characteristic magnetic field and length scale, B and L, will satisfy B 2 L 4 ∼ 2B 2 L 4 . If every blob in the system undergoes such a pairwise merger, then the total number of blobs, N , will decrease by a factor of 2: N = N/2. Assuming the blobs fill all space, their characteristic size must then increase as L = 2 1/3 L. Together, these relations imply B 2 L ∼ const, which, of course, is precisely the condition obtained from the conservation of the volume-averaged helicity density, Eq. (9).
In contrast, for a system with vanishing total helicity, there will be equal numbers of blobs with helicities −H and +H. When two blobs with opposite helicities merge, the resulting structure will be non-helical and, assuming that no higher-order topological invariants constrain its subsequent evolution, it will relax to zero magnetic en-ergy on the ideal timescale [54]. In other words, when blobs of opposite helicity merge, they mutually 'annihilate'. Otherwise, the signs of the helicities of the blobs should not modify the dynamics, so the system will have no preference between like-helicity and opposite-helicity mergers. This implies that, after one merger timescale, we will have N = N/4, because of any four randomly chosen blobs, on average, two will annihilate, and two will merge to form a single blob. Again, assuming blobs fill all space, we have L = 2 2/3 L, which, together with This new scaling is the non-helical analogue of Eq. (9). This picture, while conceptually simple, does rely on significant assumptions about the nature of the dynamics that are not obviously justifiable, e.g., that all structures have the same length scale and helicity, that the only dynamical process is mergers (rather than, say, fragmentation due to MHD instabilities) and that all mergers are pairwise processes. We now discuss how to formalise it.
B. Invariant-based theory of non-helical decay

The Saffman helicity invariant
We propose that the general evolution of a collection of localised structures of mixed magnetic helicity should conserve the integral where, as before, h = A · B is the helicity density, and angle brackets denote an ensemble average. The form of this integral is immediately reminiscent of the integrals that govern hydrodynamic decay: the already-introduced Loitsyansky integral, Eq. (1), and the Saffman integral [50], The Saffman integral is finite for, and conserved by, hydrodynamic turbulence initialised with strong long-range spatial correlations, corresponding to a kinetic-energy spectrum ∝ k 2 at the largest scales [50] (we shall return to the connection with small-k part of the spectrum below). Such turbulence is known as 'Saffman turbulence' [7]. When strong correlations are absent, I P = 0, and the Loitsyansky integral, Eq. (1), is conserved instead [7,8,50]. This case is known as 'Batchelor turbulence', after Batchelor and Proudman, who explored its properties in [55]. Physically, the conservation of I P is related to the conservation of linear momentum, P ≡ d 3 x u, much as the conservation of the Loitsyansky integral is related to angular-momentum conservation, and, as we are about to show, I H is related to helicity conservation. Owing to this analogy, of which we shall make further profitable use in Section IV, we will refer to I H as the 'Saffman helicity invariant'. We proceed by establishing the claims that, in non-helical turbulence, I H is gauge invariant, finite, and conserved. These arguments are in most respects analogous to those originally made by Saffman for I P [50] (see [7] for a review).
Assuming that volume and ensemble averages are the same, so I H is the density of magnetic helicity squared. Gauge invariance is then guaranteed in the same manner as for the magnetic helicity, by arranging that the surface of the volume V is always normal to the magnetic-field direction. This can always be achieved because of our assumption that the magnetic field forms localised structures, which are arbitrarily small compared to V in the limit V → ∞.

I H is finite.
Returning to the definition of I H , Eq. (22), and assuming that the system has no preference for accumulation of like-or opposite-helicity structures, h(x)h(x + r) is zero if r extends beyond the characteristic size of a helical structure, L, and of size ∼ h 2 ∼ B 4 L 2 otherwise. Integrating over r then gives I H ∼ B 4 L 5 . Formally, this requires the two-point magnetic-helicity-density correlation function to decay faster than r −3 as r → ∞. This is the condition for the magnetic structures to be 'localised'.
Another way to obtain this scaling is to consider the volume integral in Eq. (24) as a random walk in the net helicity contained within the volume V -the number of 'steps' is V /L 3 , so Then there is cancellation of V in Eq. (24), and the scaling I H ∼ B 4 L 5 is recovered.

I H is conserved.
According to Eq. (25), the expectation value of the square helicity in a large volume of non-helical turbulence is H 2 V ∝ V . As η → 0 + , the magnetic helicity is conserved during all processes that occur locally inside the volume V , so H 2 V ∝ V can only be changed by processes occurring at the surface of V , S = ∂V . These are fluxes of helicity in or out of the volume, or else reconnection with magnetic structures not contained within the volume. However, both are random processes and hence the net rate of change of square helicity associated with them scales as ∼ S ∼ V 2/3 . Ultimately, then, we find with the result that in the limit V → ∞, H 2 V is a conserved quantity. Therefore, so is I H .
An alternative proof of the invariance of I H , which follows directly from the MHD induction equation under the assumption of sufficiently rapid decay of long-range correlations, is detailed in Appendix B.
In the context of this work, the primary importance of the conservation of I H is that it implies precisely the same scaling as we found from our qualitative theory, Eq. (21). The more formal argument proposed here is much more general, however: while it requires that helical structures be localised (i.e., that local correlations decay sufficiently quickly with distance), and that there be no preference for accumulation of like-helicity or oppositehelicity structures, it does not require that all structures be of the same size or magnitude of helicity at any given instant, and neither does it require that the only relevant dynamics be pairwise mergers. The cartoon presented in Section III A, should, therefore, be considered as a particular example of dynamics that would conserve I H , and that, therefore, must produce the correct scaling, but not as the only dynamics allowed in the system, or required to be prevalent in order for the scaling (21) to hold.

Permanence and impermanence of the large scales
We now discuss some important consequences of the invariance of I H for the small-k asymptotics of the spectra. In particular, we shall find that the phenomenon of 'inverse transfer' of magnetic energy [34,35], whose explanation has so far been unclear, follows naturally from the conservation of I H . We shall also find that the invariance of I H implies an invariant small-k asymptotic of the helicity-variance spectrum, a fact that we shall utilise in Section III D to provide a numerical test of our theory.
As is well known, the Saffman and Loitsyansky integrals are, respectively, proportional to the coefficients of k 2 and k 4 in the small-k asymptotic expansion of the kinetic-energy spectrum. To see why, note that if correlations between distant points decay sufficiently quickly, viz., if u(x) · u(x + r) < O(r −5 ) as r → ∞, then the kinetic-energy spectrum, may be Taylor-expanded in kL 1, which yields (under the assumptions of statistical isotropy and homogeneity), Thus, Saffman turbulence, with I P = 0, has E(k → 0) ∝ k 2 [50], while Batchelor turbulence, with [55]. Owing to the invariance of I P and I L , Eq. (28) leads to a phenomenon known as the 'permanence of the large-scale eddies'as hydrodynamic turbulence decays, the small-k part of its energy spectrum remains unchanged. In particular, this means that non-helical hydrodynamic turbulence supports no inverse energy transfer. This is unlike nonhelical MHD, which, if initialised with E(k → 0) ∝ k 4 (or steeper, like in our simulations -see Appendix D), has been found in simulations to increase its energy content at large scales [34,35]. This may be interpreted as a consequence of the non-invariance of the magnetic equivalent of the Loitsyansky integral, which is related to the magnetic energy spectrum via the expansion analogous to Eq. (28), where I B will be defined and discussed in Section IV C, but for now can be assumed to be zero. Of course, there was no reason to suspect that I L M should have been a dynamical invariant, as angular momentum does not have a magnetic equivalent -L M ≡ r × B is not a conserved quantity in MHD. In fact, the growth of I L M , and hence the inverse transfer, can be recovered immediately from the conservation of the Saffman helicity invariant, as so if I H is conserved while B 2 decays, I L M must grow.
Owing to the presence of the inverse transfer, it would appear that there is no 'permanence of the large scales' principle for non-helical MHD. However, a modified version of this principle may be obtained by noting that I H is proportional to the coefficient of k 2 in the small-k expansion of the helicity-variance spectrum, Namely, at small k, provided that h(x) h(x+r) < O(r −3 ) as r → ∞. Thus, non-helical MHD does have a kind of 'permanence of the large scales' phenomenon, though it manifests itself in the helicity-variance spectrum, not in the energy spectrum. Detecting this phenomenon numerically will be a useful test of our theory and a way of confirming the conservation of I H (see Section III D).

C. Decay laws
Let us now compute the decay laws associated with the conservation of I H . Eq. (15) with α = 4/5, as demanded by Eq. (21), implies a magnetic-energy decay if the dynamics occur on the Sweet-Parker timescale (with Laplacian viscosity), or Eq. (12) gives if reconnection is fast, i.e., either stochastic or plasmoiddominated. These exponents are both close to −1, so are consistent with previous numerical results that have reported a power law decay ∝ t −1 [30,[33][34][35]38]. For the kinetic energy, we again expect Eq. (18) to hold, provided the kinetic energy is dominated by reconnection outflows. Under the Sweet-Parker scalings and B 4 L 5 ∼ const, this becomes For n = 2 and n = 4, this gives respectively. The closeness of these exponents to 1 indicates that the current-sheet aspect ratio changes more slowly in the non-helical case than the helical case, and explains why no significant difference in the kinetic-and magnetic-energy decay laws has been reported in the previous numerical studies cited above.

D. Numerical results
To test the theory proposed in Sections III B and III C, we now present results from simulations of decaying nonhelical turbulence, analogous to those presented for helical fields in Section II C.
We first address the question of the scaling and conservation of I H in our simulations. In a periodic box, rather than infinite space, one should interpret the limit V → ∞ in Eq. (24) as requiring V to be large compared to the energy-containing scales, but small compared to the box size, where the assumption of isotropy fails, as does the random-walk scaling of magnetic helicity, Eq. (25), if the total helicity in the box is constrained to be exactly zero by the initial condition. Fig. 5 shows plots of H 2 V /V vs. R, where V is a cube of width 2R, and we take an ensemble average over many different cube positions in the simulation box (employing the Coulomb gauge, ∇ · A = 0, for numerical convenience). As R → 0,  15), for decays occurring on the Sweet-Parker timescale, with n = 2, n = 4 and n = 6 shown in blue, red, and magenta respectively. The grey solid curve depicts the scaling given by Eq. (12). Simulation results are in excellent agreement with the coloured curves, with better agreement as α increases towards the limiting value of 4/5. As in the helical case, simulations at the largest Lundquist numbers appear to be on the brink of the transition to the fast-reconnection regime. The full set of decay curves from which this plot was derived can be found in Appendix D.
ing the strong scaling of I H with B and L -I H ∼ B 4 L 5 -this decay is very slow, i.e., Eq. (21) holds well. For comparison, if anastrophy were conserved, as in two dimensions, then, according to the robust Sweet-Parker scaling for n = 4, viz., BL ∼ const, B 2 ∼ t −3/4 (see Appendix A), I H would grow : B 4 L 5 ∼ t 3/8 (and, indeed, even faster growth should be expected under the oftenassumed B 2 ∼ t −1 , L ∼ t 1/2 , since then B 4 L 5 ∼ t 1/2 ). Thus, the distinction between our theory of non-helical decay and the conjecture of quasi-two-dimensional dynamics [30,34,38] is measurable in numerical simulations, and there is strong evidence in support of the former over the latter. The conservation of I H may also be demonstrated from the invariance of the small-k part of the helicity-variance spectrum, Θ(k), as explained in Section III B 2. The evolution of this spectrum is shown in Fig. 6. The small-k part exhibits a k 2 power law, indicating that the expansion Eq. (33) is valid, and that I H is finite. As the turbulence decays, the small-k part is preserved to good approximation, though there is a small amount of decay, owing to the finite scale separation between the box scale (close to which, Θ(k) ∝ k 2 should fail) and the energycontaining scales. Nonetheless, the behaviour is once again markedly different from what should be expected  Fig. 2. As in the helical case (Fig. 2), the peak of the magnetic energy is at much larger scales (smaller k) than the peak of the kinetic-energy spectrum, consistent with the expectation that the kinetic energy should be contained within the Sweet-Parker sheets. Inset: decay of the total magnetic and kinetic energies. The total kinetic energy curve is much below the magnetic energy curve, and coincides with (δ/L)EM (black), in agreement with Eq. (18). Here, we compute δ/L as the ratio of the wavenumbers at which kEM (k) and kEK (k) peak. Bottom panel: The kinetic-energy decay exponent, pK , against the magnetic-energy decay exponent, pM , as measured in simulations with S 1/n n, 0 > 9.0. Results are in reasonable agreement with the theoretical prediction, Eq. (36) (coloured lines), though as noted in the main text, this prediction is very similar to EK ∝ EM (grey).
under the often-assumed B 2 ∼ t −1 , L ∼ t 1/2 laws. In that case, an inverse transfer of helicity variance should be expected, with the small-k part of Θ(k) growing like B 4 L 5 ∼ t 1/2 .
Turning to the measured decay laws, we note that, as in the helical case, the asymptotic scaling B 4 L 5 ∼ const will not necessarily be satisfied for a decay at finite η n .
where LV is the angular momentum contained in a sphere of radius R, about the centre of the sphere. We find |LV | 2 ∼ R 5 , at late times, though the scaling is closer to R 4 initially. (c) |PV | 2 , where PV is the linear momentum contained in a cube with V = (2R) 3 . We find |PV | 2 ∼ R 2 , though the scaling appears to become somewhat stronger with time.
Nevertheless, we can still measure a value of α for which B α L ∼ const is satisfied, and compare the measured value of the magnetic-energy decay exponent, p M , with the one expected for a decay on Sweet-Parker or ideal timescales. The results of this comparison are shown in Fig. 7. While agreement with the Sweet-Parker curves (coloured) for α < 4/5 is not quite as good as in the helical case, we still observe (i ) a clear preference for the Sweet-Parker prediction over the prediction of a decay on ideal timescales (grey), and (ii ) convergence to B 4/5 L ∼ const for the hyper-dissipative simulations.
Finally, we describe the decay of kinetic energy. Fig. 8 again shows that, like in the helical case, the kinetic energy is peaked at smaller scales than the magnetic energy is, consistent with the expectation of Alfvénic outflows in current sheets. We also find that Eq. (18) is very well satisfied, and a reasonable agreement with Eq. (36), although the magnetic-and kinetic-energy decay exponents are very close to each other.

E. The behaviour of other invariants
To conclude the discussion of numerical results, we now address the evolution of the other, better known invariants during our simulations, namely, the cross-helicity, as well as the Loitsyansky and Saffman integrals.
First, we consider the cross-helicity, which is an ideal invariant of the (incompressible) MHD equations, though we find that it is not conserved in our simulations any better than energy is. In Fig. 9(a), we plot the average value of the squared total cross-helicity contained in a cube of volume V = (2R) 3 , in a manner analogous to Fig. 5 for the magnetic helicity. We find that H 2 c V ∼ V for R π/2, as is expected from the randomwalk argument. For R > π/2, H 2 c V /V increases with R but not as fast as R 3 , which indicates that there is no net cross-helicity in the box. Even so, one may consider the Saffman-type 'invariant', I Hc , that is associated with cross-helicity (see Sections IV D and IV E). It is given by the value of H 2 c V /V in the flat region of Fig. 9(a), and, as we see, decays. The reason for this behaviour is that the cross-helicity's decay rate has the same scaling with the dissipation coefficients η n (resistivity) and ν n (viscosity) as the decay rate of the magnetic energy does (and as is inevitable, because cross-helicity and energy have the same physical dimensions).
Likewise, the Loitsyansky integral, which encodes the statistics of angular-momentum fluctuations [6,56], is not conserved in our simulations, as is clear from the small-k part of Fig. 8. In Fig. 9(b), we plot |L V | 2 against R, where L V is the total angular momentum contained in a spherical control volume V of radius R, calculated about its centre. Presumably, this is due to injection of angular-momentum fluctuations by the reconnection outflows. Intriguingly, we observe a shift in the scaling properties of |L V | 2 vs. R -at early times, |L V | 2 ∼ R 4 , which is the expected scaling for 'Batchelor turbulence', i.e., when correlations between distant points are weak [7,55]; subsequently, the system appears to evolve towards a state with |L V | 2 ∼ R 5 . The latter is the 'Saffman-turbulence' scaling, and suggests strong long-range correlations in the velocity field [7,50]. A corresponding shift in the analogous quantity for linear momentum, |P V | 2 , vs. R is suggested by Fig. 9(c), which shows a stronger scaling than R 2 (the Batchelor scaling) for R < π/2 at later times, closer to the Saffman-turbulence scaling of R 3 . This is also consistent with Figs. 2 and 8, which appear to show a decreasing slope in the large-scale kinetic-energy spectrum over time, perhaps towards E K ∝ k 2 , the hallmark of Saffman turbulence [7,50]. This spectral behaviour has been noted in other studies [30,34,46], though it was considered an effect of compressibility, owing to the fact that the incompressible simulations of [29] appeared not to see it. However, this may just have been because there was not enough time for the k 2 velocity spectrum to establish itself before the outer scale of the turbulence reached the box size in that study. While not present in decaying hydrodynamic turbulence [7], we suggest the effect might be related to the 'thermalisation' phenomenon that is observed in forced, hydrodynamic turbulence [57][58][59][60]. We shall address this topic specifically in a future publication [61].
Finally, for the reader concerned that the standard scalings for |L V | 2 and |P V | 2 referred to here are not the same as the ∝ R 3 random-walk scalings assumed in Section III B, we show how these scalings may be obtained from the random-walk approach in Appendix C. A more formal derivation of them may be found in [7].

IV. DISCUSSION
A. Case of small, but non-zero, helicity In Section III, we proposed a way to impose the constraint of magnetic-helicity conservation on the decay laws of non-helical MHD turbulence, via the conservation of I H . Of course, no real field configuration will have precisely zero helicity, and therefore it is important to consider the evolution of a field configuration with small, but non-zero, magnetic helicity. In such a case, the system will undergo a transient (though, perhaps, long) period of evolution according to the non-helical decay law B 4 L 5 ∼ const [Eq. (21)], before ultimately entering the fully helical regime, with a corresponding change in the decay law to B 2 L ∼ const [Eq. (9)].
This conclusion is an immediate consequence of the non-helical decay laws, as follows. Suppose that the system starts with some small helicity fraction σ 0 1, defined so that the total helicity is H = σ 0 B 2 0 L 0 V . At a later time, since helicity is conserved, Because σ 0 1, the system is not controlled by its total helicity, at least initially. It therefore evolves according to B 4 L 5 ∼ const. Using this in Eq. (38), we find Thus, the helicity fraction σ will grow with time. This can continue until σ ∼ 1, at which point the system enters 10. Schematic of H 2 V /V as function of V , as a system with small initial fractional helicity, σ0, transitions to the fully helical regime, as explained in the main text. The progression of time is shown by red → blue, with plots at logarithmically spaced time intervals. Note that both axes are plotted on logarithmic scales. the helical regime. This occurs when provided that the energy-containing scale L has not yet reached the system size. This result is intuitive from the cartoon picture presented in Section III: when blobs with one sign of helicity are more populous, the ultimate consequence of random mergers is for the less populous type to be used up (although this can take a long time, which scales with an appropriate negative power of σ 0 , if the initial population imbalance is small).
The same conclusion can be reached from the consideration of the Saffman helicity invariant, Eq. (24), though some care should be taken, as I H is formally infinite in the presence of any net helicity. However, if the helicity fraction σ is small, then we can interpret the limit V → ∞ in Eq. (24) as requiring L 3 V V c , where V c is the critical volume at which H V ceases to be dominated by the net helicity fluctuation owing to its collection of helical structures of random signs, and instead is dominated by the helicity imbalance (see Fig. 10). This condition implies For any V such that L 3 V V c , the arguments presented in Section III B in favour of the conservation of H V remain valid and hence I H still provides the dominant constraint on the decay of magnetic structures. Choosing instead V > V c , H 2 V ∼ σ 2 V B 4 L 2 ∼ const, so we recover the evolution equation for σ, Eq. (38). However, when ultimately σ ∼ 1, there is no longer any possibility of satisfying L 3 V V c , because V c ∼ L 3 . At this point, H 2 V ∼ V B 4 L 2 for any chosen volume, and we are back to the fully helical scaling, Eq. (9).
These arguments suggest that the non-helical decay is ultimately transient for any real system, and there will always eventually be a transition to the helical regime, provided that the growing, energy-containing scale does not reach the system size before this happens. The same conclusion was reached by [1], although their argument was based on different non-helical decay laws than those that we have proposed. A numerical simulation demonstrating the expected transition between the two regimes would be highly desirable, but imposes considerable numerical cost, so is left for future work.

B. General decay principles
Let us now discuss how the principles that guided us in the above might be applied to other types of decaying turbulence. In any type of turbulence with an ideal invariant, F = d 3 x f (x), that is better conserved than the energy, the conservation of that invariant implies This is, in general form, the 'selective decay' principle of [14][15][16][17][18][19][20]. If nonlinear structures (eddies, blobs) are necessarily in possession of F , so that the characteristic size of f can be related to the sizes and correlation scales of the dynamical fields, then Eq. (42) imposes a constraint on these fields that must be satisfied as the turbulence decays. In particular, if the sign of f (x) within each structure is the same, either because it is sign-definite, or because the initial condition stipulates a predominance of structures of one particular sign, then f can be related to the characteristic sizes of dynamical fields and their correlation scales: where ψ is a representative dynamical field, L ψ is its correlation scale, a and b are exponents determined by the functional form of f , and '. . .' denotes the possibility of other fields and scales. Together, Eqs. (42) and (43) imply a constraint on the dynamical fields that must be satisfied during the decay, Alternatively, for sign-indefinite f , there may be no strong predominance of structures associated with either sign of the invariant. In this case, where σ 1 is the fractional imbalance in F , as in Section IV A. While Eq. (42) still holds, its utility is reduced as σ is generally a function of time -the relationship implied by Eqs. (42) and (45), should be considered as an evolution equation for σ, rather than a constraint on ψ and L ψ . However, all is not lost. We propose that when σ 1, conservation of local fluctuations in F imposes a constraint on the decay, through the associated Saffman-type integral (47) where L f is the correlation scale of f . Again, if nonlinear structures are in possession of F , then f 2 L 3 f can be related to the characteristic sizes and correlation lengths of dynamical fields, where c and d are different exponents than those in Eq. (43). Together with Eq. (47), this scaling leads to a different constraint on the decay, which is independent of the fractional imbalance σ. We note that an eventual transition from the balanced regime [Eqs. (46) and (49)] to the imbalanced regime [Eq. (44)] is a general consequence of these results, because, together, Eqs. (46) and (47) imply A necessary and sufficient condition for the fractional imbalance to grow is, therefore, that the scale L f should increase with time. This tends to be the case for any realistic decay problem, even in the absence of inverse transfer of energy in k-space, because small-scale structures generally dissipate faster than large-scale ones. The correspondence of Saffman-type invariants to the small-k asymptotics of spectra, as explained for I H in Section III B, is also a general property: the invariant I F is related to the spectrum of the variance of f (x), provided that f (x) f (x + r) < O(r −3 ) as r → ∞, and that the system is statistically isotropic (though this statement is easily reformulated for anisotropic systems). Therefore, systems that decay while respecting the conservation of a Saffman-type invariant generally have a 'permanence of large scales' principle that applies to the spectrum of the variance of the relevant conserved quantity. This principle provides a convenient way to assess the existence and conservation of candidate invariants in numerical studies (as we did in Fig. 6).
To conclude this section, we note that our discussion has relied on two conditions: (i ) that the invariant, F , is conserved, i.e., that its rate of change is smaller than the energy decay rate, and (ii ) that it can always be considered the case that typical structures possess F , so that f and f 2 may be related to the sizes and correlation lengths of the dynamical fields, via Eqs. (43), (45) and (48). Establishing whether (i ) and (ii ) holds for any given turbulent system requires some physical idea of the decay dynamics. For example, we found in Section III E that the Saffman-type invariant corresponding to crosshelicity was not conserved in our numerical simulations. Nonetheless, different dynamical processes can result in invariants being dissipated at different rates, and indeed there are plausible reasons to argue that cross-helicity (or its Saffman-type-invariant counterpart) might be conserved by the decaying turbulence of interacting nonlinear Alfvén waves, as we will explain in Section IV D.
Dynamical processes may also govern whether structures possess the relevant invariants, condition (ii ). For example, we have argued that in decaying MHD turbulence, non-helical magnetic structures tend to relax to zero energy (as is consistent with J.B. Taylor relaxation), so that at any given time, individual magnetic structures are maximally helical.
As a final general remark, we point out that it is possible for nonlinear structures to possess more than one invariant. If the constraints implied by the conservation of these invariants are not mutually exclusive, then they must be satisfied simultaneously. If, however, these conditions are contradictory, then it will be necessary for some of them to be broken. In that case, the decay cannot take place on the characteristic nonlinear timescale ("eddy-turnover time"), but must instead take place on the timescale on which the weaker constraint (in the sense of quality of conservation) can be broken. This is precisely the situation in magnetically dominated turbulence: magnetic helicity is not the only topological invariant associated with the magnetic field [62], and in principle all higher-order topological invariants might impose constraints on the decay. Conserving all topological invariants is not consistent with a decay in the magnetic energy, however, because it implies B ∼ const. Therefore, the decay can only proceed on the timescale on which the higher-order topological constraints may be violated, i.e., the reconnection timescale [52]. This is illustrative of a general rule that stronger, consistent constraints set the scalings between the integral scales and energies that control the decay, while inconsistent constraints set the decay timescale.

C. When the Saffman helicity invariant fails
Let us now apply the insights of Section IV B to determine the conditions under which the conservation of the Saffman helicity invariant does not impose a constraint on the decay of isotropic, non-helical MHD turbulence.
Naturally, this will be the case if it is zero, which is possible if each magnetic structure is individually nonhelical, i.e., h is pointwise zero, so condition (ii ) in Section IV B is violated. An example would be an ensemble of untwisted, unlinked magnetic tori, such as arise in the final stage of process (b) shown in Fig. (4) [63]. Such a configuration would not be prevented by the invariance of its magnetic topology from relaxing quickly under ideal dynamics, transferring magnetic energy to kinetic. The evolution should then be constrained by invariants pertaining to the velocity field, such as the Loitsyansky integral or cross-helicity -plausibly, the latter may be conserved in a net [cf. Eq. (42)] or fluctuating [cf. Eq. (47)] sense in turbulence with U ∼ B and similar integral scales for both fields, as we shall explain in Sections IV D and IV E.
While we thus acknowledge the possibility that I H can be zero, we view such a field configuration as rather artificial. In any realistic physical scenario, it seems likely that some non-zero fraction of magnetic structures should possess helicity, whether due to twists or to linkages, so I H = 0. Even if this fraction is not 1, we expect that individually non-helical structures should relax on ideal timescales (cf. Fig. 4), and the evolution of those that remain should then be constrained by conservation of I H (also see [64] for a study of the tendency of MHD turbulence to form local helical patches).
The other scenario in which conservation of I H will not impose a constraint on the decay is if it diverges, which it will do if the magnetic-helicity correlation function h(x)h(x + r) decays as r −3 or slower as r → ∞. Physically, I H diverges when our assumption of 'localised magnetic structures' fails. This should be expected if is non-zero, which, as we saw in Eq. (30), corresponds to E M (k → 0) ∝ k 2 . Such a turbulence has long-range (longitudinal) correlations in the magnetic field that decay as r −3 as r → ∞ (the reason for this is the same as the one for which long-range correlations in the velocity field are required for I P = 0 -we refer the interested reader to [7]). Magnetic-helicity correlations should be expected to decay at least as slowly as this (A being formally an integral of B), in which case I H will diverge. So, what principle governs the decay of non-helical turbulence of this sort? We suggest that it should be constrained by the conservation of fluctuations in the magnetic flux, under the formalism described in Section IV B. In this case, I B itself is the relevant invariant, which we might call the 'Saffman flux invariant' [cf. Eq. (47)], as where Φ V = V d 3 x B is the total flux contained within the control volume V . In light of Eq. (30), we note that non-helical turbulence with E M (k → 0) ∝ k 2 should not have an inverse energy transfer; its energy spectrum should instead obey a 'permanence of the large scales' principle. This prediction is consistent with the results of the extensive parameter study by [35]. In decay from a magnetically dominated state, (i.e., α = 2/3) corresponds to a decay law of if the decay occurs on ideal timescales [see Eq. (12)]. This is, inevitably, the same as Saffman's law for the decay of kinetic energy in hydrodynamic turbulence with E(k → 0) ∝ k 2 [50]. Assuming that fast-reconnection outflows are the dominant motions, Eq.
, so E K obeys the same law as E M . If reconnection occurs slowly, in the Sweet-Parker regime, then Eqs. (15) and (19) give for Laplacian (n = 2) dissipation, where the latter law again assumes the dominant motions to be the reconnection outflows. These are interesting predictions to test. However, while it is possible to initialise turbulence in a state with I B = 0, the requirement of strong long-range correlations in B would appear to make it somewhat artificial. For example, in cosmological contexts, causality constraints are expected to rule out spectra shallower than the 'Batchelor' E M (k → 0) ∝ k 4 [30,34,35,65], which, as explained above, corresponds to localised magnetic structures (i.e., rapidly decaying correlations).

D. Decay of MHD turbulence in the presence of strong mean field
In this work, we have so far restricted attention to the decay of isotropic MHD turbulence, i.e., MHD turbulence without a mean magnetic field. The case of turbulence with a strong mean field is fundamentally different, because then the magnetic helicity is not a conserved quantity. Formally, this is because the mean field 'sticks out' of any volume for which one might choose to compute the magnetic helicity. Intuitively, also, this is a different situation compared to isotropic turbulence, because the constraints imposed by topology are much reduced. For example, purely magnetic structures need not persist until they are able to reconnect with each other, instead they can relax by decomposing themselves into Alfvén waves travelling in opposite directions along the mean field.
Aside from energy, MHD turbulence with a strong mean field (described by the 'reduced MHD' equations [66,67]) has only one conserved quantity related to the presence of the magnetic field, the cross-helicity, where u ⊥ is the fluid velocity perpendicular to the mean magnetic field, and b ⊥ is the magnetic-field perturbation. Like magnetic helicity, the cross-helicity is sign-indefinite. In simulations of driven MHD turbulence with a strong mean field, a tendency to develop patches of strong local cross-helicity is observed, even in so-called balanced turbulence where the volume-averaged cross-helicity density is zero [13,14,64,68,69]. The reason for this is that structures of large cross-helicity are also structures containing strong imbalance in the sizes of the two Elsässer fields, Z ± = u ⊥ ± b ⊥ . Importantly, individual Elsässer fields each represent exact nonlinear solutions to the MHD equations. Nonlinearity, and hence turbulent decay, can, therefore, only be present where both fields are present. Since u ⊥ · b ⊥ = (|Z + | 2 − |Z − | 2 )/4, a small cross-helicity density indicates balance between the Elsässer fields, and, therefore, large nonlinearity. Such structures are prone to turbulent decay. In contrast, structures with strong cross-helicity of either sign have reduced nonlinearity and, therefore, are more immune to turbulent decay. Note, however, that these considerations need not apply to isotropic, magnetically dominated MHD turbulence decaying via reconnection, because the possession of cross-helicity does not afford immunity to reconnection. For balanced, Reduced-MHD turbulence, however, they motivate us to conjecture that the decay might be controlled by the 'Saffman cross-helicity invariant' where h c = u ⊥ · b ⊥ [70]. We note that, like I H , I Hc is an example of an invariant that depends on a fourthorder correlation function. The relevance of fourth-order correlators to distinguishing between different species of decaying MHD turbulence has previously been suggested by [11], inspired by the numerical results of [71]. One might expect I Hc to be finite and conserved by precisely the same arguments as we presented for I H , the Saffman helicity invariant, in Section III. By a randomwalk argument analogous to the one in Section III B, we have where l and l ⊥ are the characteristic lengthscales parallel and perpendicular to the mean field, respectively. For Alfvénic motions, b ⊥ ∼ u ⊥ . Note that this scaling is on much firmer ground in the mean-field case than the isotropic case, because of the absence of topological constraints associated with helicity conservation. The parallel and perpendicular length scales can be related to each other by the conjecture of critical balance [37], which states that b ⊥ /l ⊥ ∼ B 0 /l , and is a cornerstone of the theory of strong MHD turbulence [72][73][74]. Critical balance is essentially a statement of causality: the characteristic parallel length scale of an eddy cannot be longer than the distance travelled by an Alfvén wave in one nonlinear timescale. It has been confirmed numerically to great precision in driven RMHD turbulence [75], and appears to be satisfied in decaying turbulence too [37]. Putting these scalings together, we find that I Hc ∝ b 3 ⊥ l 3 ⊥ . If I Hc is indeed conserved in decaying MHD turbulence, this implies Then, since the energies of the Elsässer fields are comparable in balanced turbulence, the turbulent decay would likely be governed by This results in a decay of both magnetic and kinetic energy Intriguingly, this decay law has indeed been observed in simulations of decaying, balanced RMHD turbulence, though was rationalised differently, by assuming local effective conservation of anastrophy, which also implies Eq. (63), as in two-dimensional MHD turbulence [37] (see Appendix A). If indeed the decay of MHD turbulence in the presence of a strong mean field conserves cross-helicity, any initial imbalance will eventually lead to a final state with maximal cross-helicity, i.e., a pure Elsässer state, according to the general argument presented in Section IV B [see Eq. (50)]. Such a state will not decay, since it is an exact solution of the non-linear RMHD equations. Such behaviour has indeed been reported in numerical studies [4,[76][77][78].

E. Decay of isotropic MHD turbulence from an initial state with U ∼ B
Another problem to which the formalism developed here may be usefully applied is to the decay of isotropic MHD turbulence from an initial state with U ∼ B, as opposed to the B U that we have so far considered in this work (apart from in Section IV D). Such a state is the natural final state of the MHD dynamo (see, e.g., [79,80] for reviews). In such a case, we conjecture that the simultaneous conservation of magnetic helicity and cross-helicity might be respected by the decay, as the constraints they imply are not mutually exclusive. Such a prospect has been considered by [9,11,14,16,17] for net-cross-helical initial states -here, however, we shall impose cross-helicity conservation via the Saffman-typeinvariant formalism developed in Section IV B, reflecting the fact that no strong net cross-helicity is generically present in MHD turbulence without a mean magnetic field. Naturally, checking the circumstances under which conservation of Eq. (59) (with h c = u · B in the isotropic case) is valid will require a detailed numerical study, which is left for future work. However, the consequences of this conjecture merit discussion here, because they do appear to be consistent with already-existing numerical studies, and the argument demonstrates a general principle of simultaneous conservation of multiple invariants.
Let us consider a system that, as a result of dynamo or some other process, has reached equipartition between magnetic and kinetic energy, U ∼ B, with the same integral scale L.

Helical magnetic field
First, let us assume that the magnetic field is helical, but that there is no predominance of either sign of crosshelicity. Then the conjecture of simulataneous conservation of magnetic helicity and cross-helicity (the latter as a Saffman-type invariant) implies B 2 L ∼ const [Eq. (9)] and B 2 U 2 L 3 ∼ const, respectively. Together, these conditions imply U ∼ B 2 , or precisely the condition found numerically by [22,26], though without theoretical justification. They conjectured that the decay should take place on the timescale associated with the u · ∇u nonlinearity in the MHD equations, i.e., L/U , which is consistent with the idea that the timescale associated with the magnetic nonlinearity B · ∇B is effectively lengthened by topological constraints on the magnetic field. It is readily verified that such a decay leads to the power laws as found numerically by [22,26,29]. The more rapid decay in kinetic energy will result in a state with U B. The system should then decay in the strong-field regime described in Section IV D. Denoting the small perturbations to the newly established strong magnetic field by δB ∼ U , we should, according to Eq. (63), then have δB 2 ∼ U 2 ∼ t −1 , assuming the decay is critically balanced. Meanwhile, B 2 should decay according to the reconnection-controlled law, either t −2/3 for fast reconnection, or B 2 ∼ t 2n/(5n−3) [Eq. (16)] for Sweet-Parker reconnection. For the n = 4 simulations employed by [22,26], this implies which will persist until the Sweet-Parker outflows dominate the kinetic energy, at which point the kinetic-energy law should change to the one given by Eq. (19), which, independently of the type of reconnection, is always slower than t −1 . The decay laws given by Eq. (66) are very close to those of Eq. (65), and, therefore, either might explain the laws found numerically by [22,26]. Regardless, a more rapid decay of the kinetic than magnetic energy appears to be robust, as does the corresponding establishment of magnetically dominated state. Indeed, a magnetically dominated final state was observed in the simulations of [46], despite being initialised with U B, with decay laws similar to Eq. (65) in the transient U ∼ B regime.

Non-helical magnetic field
Alternatively, let us consider the case of a non-helical magnetic field, initially in equipartition with the kinetic energy. Then, instead of the condition B 2 L ∼ const, we have B 4 L 5 ∼ const [Eq. (21)]. Together with the constraint implied by the conservation of the Saffman cross-helicity invariant, B 2 U 2 L 3 ∼ const, this implies U ∼ B 1/5 , or Unlike Eq. (64), Eq. (67) implies a much faster decay of the magnetic energy than the kinetic energy. However, this decay will be short lived, because the magnetic energy can be maintained at a small, but finite, fraction of the kinetic energy by dynamo. Nonetheless, the magnetic field will always remain just below dynamical strength, because if it were to grow to dynamical levels, cross-and magnetic helicity conservation would force it to decay rapidly.
In the absence of a dynamical-strength magnetic field, the kinetic energy should decay according to the purely hydrodynamic Kolmogorov law, E K ∼ t −10/7 [Eq. (3)]. Because the magnetic energy is tied to a finite fraction of the kinetic energy by the competing effects of dynamo and simultaneous cross-and magnetic-helicity conservation, it must also be the case that E M ∝ t −10/7 , i.e., Such a decay of E M and E K has indeed been observed in non-helical simulations (though initialised with U B), by [38], together with evidence of dynamo action. We note, however, that it is possible that the reason E M never grew to equipartition in these simulations was, rather, that the efficiency of the dynamo was reduced in the absence of the mean-field dynamo effect associated with helical velocity fields (see [79] and references therein).
The arguments presented in this section, if correct, suggest a remarkable general principle: an initially helical velocity field will, due to its tendency to grow a helical magnetic field through mean-field dynamo action, eventually decay in a magnetically dominated state, with E M , E K ∝ t −2/3 (in the fast reconnection regime). In contrast, non-helical velocity fields will always remain in a kinetic-energy-dominated state, with E M , E K ∝ t −10/7 .

F. Conclusion
As Sections IV D and IV E illustrate, the Saffman-typeinvariant approach appears to be an extremely useful general tool that allows consideration of sign-indefinite invariants in the 'selective decay' framework for decaying turbulence [81]. There are many different types of fluid turbulence to which the approach may be usefully applied -a number of them are reviewed in [7], and there are likely to be others, especially in the large variety of plasma systems increasingly of interest in the context of various types of space or astrophysical turbulence (see, e.g., [82][83][84]). The invariant I H , proposed in this work for incompressible magnetohydrodynamic turbulence, should remain an invariant even in the more complex cases of compressible [85,86], relativistic [33] and/or kinetic dynamics [87][88][89][90], since magnetic helicity remains a conserved quantity in such contexts. This means that the invariance of I H (and the physical principle of conservation of local magnetic-helicity fluctuations from which it follows) should provide a constraint on decaying turbulence in a wide variety of magnetised astrophysical systems (provided the dominant motions occur at scales much smaller than the system size). Possible applications include turbulence in star-forming molecular clouds [85,91] and galaxy clusters [3,92]; the generation of seed fields for galactic dynamo [37]; and the evolution of primordial magnetic fields in the early Universe [1][2][3]. With regard to the latter, we note that the non-helical decay laws that we have derived in this work [see Eq. III C] are consistent with observational constraints on magnetic fields in cosmic voids, whereas previously accepted models are not [93]. This point will be addressed specifically in a future publication [94].
In those physical systems where frozen-in magnetic fields are dominant actors in the dynamics, magnetic reconnection is likely to be the key physical process whereby the decay occurs. Marrying this insight to the constraints imposed by invariants appears to be a winning strategy for constructing decay theories, as it has proven to be in the MHD turbulence regimes that we have considered above. code. D.N.H. was supported by a UK STFC studentship. The work of A.A.S. was supported in part by the UK EPSRC grant EP/R034737/1. This work used the ARCHER UK National Supercomputing Service (http://www.archer.ac.uk).
Appendix A: Decay of two-dimensional turbulence In this appendix, we review the problem of decaying two-dimensional MHD turbulence, which, like its threedimensional counterpart, respects the conservation of an invariant associated with the topology of the magnetic field: the square of the magnetic vector potential, 14,17,21,24], sometimes called 'anastrophy'. In two dimensions, anastrophy is well defined (i.e., not gauge dependent), and evolves according to where dE M /dt = η n d 2 r B · ∇ n B is the rate of magnetic-energy decay due to Ohmic heating. Eq. (A1) implies that which is an even slower rate of change than the one we found for the magnetic helicity in three dimensions, Eq. (8). Therefore, like helicity, anastrophy in two dimensions should be conserved as the turbulence decays, for η n → 0 + . Physically, anastrophy conservation is related to the conservation of in-plane magnetic flux [44]. Unlike helicity, though, the anastrophy is manifestly positive definite, so there is only one decay regime, and no Saffman-type invariant. The conservation of anastrophy implies i.e., α = 1 in Eq. (10). According to Eqs. (12) and (15), this implies a power law decay of the magnetic energy if the decay proceeds on the 'ideal' L/B timescale, and for a decay proceeding on the Sweet-Parker timescale, (L/B) S 1/n n . Coincidentally, these laws are the same for n = 2. Remarkably, the Sweet-Parker scaling for the kinetic-energy decay also turns out to be which is independent of n. These results mean that a decay on ideal timescales produces indistinguishable power laws (E K ∼ E M ∼ t −1 ) to a decay on the n = 2 Sweet-Parker timescale. It is perhaps for this reason that over thirty years separates the derivation of the ideal decay law [21,24] and the suggestion of a decay controlled by Sweet-Parker reconnection by [44]. The picture presented in [44] is analogous to (and has inspired) the 'cartoon' model that we proposed in Section III, although we observe that, as for the cartoon model in our study, the formulation based on the conservation of integral invariants, with decay proceeding on the reconnection timescale, is more general, as it does not require that pairwise mergers between structures of equal anastrophy (or helicity) be the only allowed dynamical process. As a consequence of the degeneracy in power laws, it was confirmed in [44] that the Sweet-Parker timescale governed the decay in their simulations by showing that their evolution curves collapsed onto each other when time was renormalised by the initial Sweet-Parker reconnection timescale. An alternative method by which the Sweet-Parker-controlled decay can be established is via the use of hyper-dissipation, in the same manner as we have done in the main part of this work, thanks to hyperdissipation lifting the power-law degeneracy. For example, from Eq. (A5), the magnetic-energy-decay power laws are t −4/5 and t −3/4 for n = 4 and n = 6, respectively.
In Fig. 11(a-c), we show the evolution of the magneticenergy decay rate in our two-dimensional simulations (see Table I in Appendix D for details), normalised by the power of the energy to which it is proportional in the reconnection-based theory [i.e., the power of B 2 on the right hand side of Eq. (14), with α = 1]. As explained in Appendix D (also, see [22]), such plots are preferable to more conventional plots of log E M against log t, because they give an unbiased estimate of the decay law. On these plots, horizontal curves indicate agreement with theoretical expectations. As the resistivity decreases, and hence the Lundquist number increases, Figs. 11(a-c) show increasingly horizontal curves, in agreement with Eq. (A5). The insets to these figures show the scaling that would be expected if the decay proceeded via ideal motions [21,24]. In both cases, we find clearly decreasing curves for the largest Lundquist numbers tested, so these results are inconsistent with the 'ideal' law, Eq. (A4).
In Fig. 11(d) we show E M against t/τ rec,0 , where τ rec,0 is the Sweet-Parker timescale S 1/2 L/B at the start of the self-similar decay period, which we take to occur at t = 1 in all cases for the purposes of this calculation. While such plots are not well-suited to an accurate determination of the decay law, they do show a clear difference in behaviour between the case of Laplacian dissipation (n = 2, blue) and hyper-dissipation (n = 4, red). As previously mentioned, the collapse of the decay curves onto each other under such a normalization was presented as evidence for a reconnection-controlled decay in two-dimensions by [44]. Fig. 11(d) shows that  Table I).
the same behaviour occurs in the n = 4 hyper-dissipative case, and we find the same in the n = 6 case (not shown), so we agree entirely with [44] that the decay is controlled by reconnection.
The remaining panels of Fig. 11 show other relevant quantities besides the magnetic energy. In Fig. 11(e), we show the evolution of the hyper-Lundquist number, which is in excellent agreement with theoretical expectations based on Eqs. (A3) and (A5). Interestingly, the hyper-Lundquist number is expected to grow when n > 2, so even if the simulation starts in the Sweet-Parker reconnection regime, it can ultimately transition to the plasmoid-dominated regime. We reiterate that in such a regime, which is formidable to simulate, the reconnection timescale becomes proportional to the ideal timescale (though longer by a factor of 10 2 [39]), and then we expect a transition to the 'ideal' t −1 decay [Eq. (A4)]. Finally, Fig. 11(f) confirms that BL is indeed constant during the decay, as demanded by the conservation of anastrophy.
Consider now the decay of the kinetic energy. Fig. 12 shows plots of our best numerical estimates of p K vs. p M . While we note that most of our runs do fall close to the p K = 1 line, as predicted by Eq. (A6), the numerical evidence for p K = 1 over p K = p M is not very strong, due to the similarity between the decay laws involved. Another factor that makes this comparison difficult is the high level of noise in the decay curves (see Fig. 15), which is not present in three dimensions, and arises because of the greatly reduced number of magnetic structures in the two-dimensional simulations. Because we initialise both types of simulation with the magneticspectral-energy density peaked at the same wavenumber, k p 33 (see Appendix D), there are initially ∼ 33 times fewer structures in our two-dimensional runs.

Appendix B: Alternative proof of the conservation of the Saffman helicity invariant
In this appendix, we present an alternative proof of the invariance of I H , that follows directly from the MHD induction equation. As in Section III B, we shall rely on the assumption of rapidly decaying spatial correlations.
Restricting to the physical case of Laplacian (n = 2) dissipation, for simplicity, the MHD induction equation reads Relation between the measured magnetic-and kinetic-energy decay exponents in our two-dimensional simulations. Eq. (A6) predicts that all simulations should have pK = 1, independently of pM . While many of the simulations do fall close to this line, the evidence for pK = 1 over pK = pM is not very strong, owing to the closeness of all powers involved, and the large amount of noise in the twodimensional simulations (see end of Appendix A).

Uncurling this equation, we have
where χ is an arbitrary scalar function. From Eqs. (B1) and (B2), it follows that where, as elsewhere, h = A · B is the magnetic helicity density, while is the 'magnetic-helicity flux'. As argued in Section II A, the resistive helicity-dissipation term on the right-hand side of Eq. (B3) is small -its size is ηB 2 /δ η , giving a helicity-dissipation timescale of ∼ Lδ η /η, which is long compared to the magnetic-energy dissipation timescale, δ 2 η /η. Dropping it, Eq. (B3) is a continuity equation for magnetic helicity. Integrating over space, we obtain its conservation law. Alternatively, we may use Eq. (B3), applied at x and x + r, to compute the evolution of the two-point helicity correlation function. After taking an ensemble average, and using statistical isotropy and homogeneity, we have Integrating over r, assuming h(x)F (x + r) < O(r −3 ) as r → ∞, Eq. (B5) gives The conservation of general Saffman-type invariants (see Section IV B) may be shown from their corresponding continuity equations in precisely the same manner.
Appendix C: Random-walk scalings for linear and angular momentum In this appendix, we provide simple arguments for the scalings of |P V | 2 and |L V | 2 vs. R, based on the random-walk argument employed in the main text. These are, respectively, the expectation values of the squared linear and angular momenta contained within a control volume V ∼ R 3 . In the latter case, we take V to be spherical, and compute the angular momentum about its centre. As noted in Section III E, these quantities do not necessarily scale with V as suggested by the naïve random-walk estimates employed in Section III B. However, the correct scalings may be obtained if one accounts for incompressibility, as we now demonstrate.
First, let us consider |P V | 2 for a volume V ∼ R 3 . The random-walk argument put forward in Section III B suggests that |P V | 2 ∝ R 3 . This is the correct scaling for Saffman turbulence, where the velocity field is initialised with long-range correlations and individual eddies can have non-vanishing linear momentum [7,50]. For u initialised with a (longitudinal) correlation function that falls off with distance sufficiently rapidly (Batchelor turbulence), however, the correct scaling turns out to be |P V | 2 ∝ R 2 [7,50]. We can understand this from the random-walk argument, taking into account incompressibility, as follows: where u = ∇ × ψ, and we have assumed that S dS × ψ will sum as a random walk (an assumption that fails if the long-range correlations in u are strong). Therefore, the Saffman integral, for such turbulence. Note also that the linear-momentum fluctuation in V will not be conserved, because it is formally the same size, ∝ R, as the net surface flux that can cause it to change. For the same reason, the conservation of the total magnetic flux d 3 x B also does not provide a constraint on the decay of turbulence without strong long-range correlations in B (see Section IV C for a discussion of turbulence that does have such correlations).
Even if the Saffman integral vanishes, the local fluctuations in the linear momentum dominate |L V | 2 , over the local rotation of the eddies. The reason is that structures further from the origin will contribute more angular momentum than closer structures, leading to |L V | 2 > O(R 3 ) (this effect also means that |L V | 2 depends on the shape of V , which is why it was necessary to assume V to be spherical). In Saffman turbulence, where correlations are long-range and individual flow structures may have finite linear momentum, it turns out that |L V | 2 ∼ R 5 , owing to this effect [7]. This conclusion too can be obtained from a randomwalk argument: the expected square angular momentum in a spherical shell of radius r and width δr satisfying L δr r, is where U t is the typical size of the net translational velocity of a structure. Assuming any two shells are uncorrelated, the total square angular momentum is simply the sum over all shells of Eq. (C3). Integrating over r, we get Like Eq. (C1), the scaling (C4) is also adjusted by incompressibility when correlations fall off quickly with distance. This is because Integrating by parts and expanding the double cross product gives where we have chosen V to be spherical. Since ψ is a random field, the first integral scales as R 2 , so it dominates over the second, which scales as R 3/2 . Therefore, Eq. (C6) implies As above, whether ultimately Eq. (C4) or Eq. (C7) provides the correct scaling depends on the strength of longrange correlations between eddies. A formal derivation of these statements may be found in [7]. In either case, the scaling of |L V | 2 is different to the naïve expectation, |L V | 2 ∼ R 3 , which assumes all eddies to have no translational motion. In that case, the angular momentum of an eddy at a distance r from the origin is ∼ [(r + L)U r − (r − L)U r ] L 3 ∼ U r L 4 , where U r is the typical rotational velocity of a structure. Then Eq. (C3) becomes whence In summary, while the scalings of |P V | 2 and |L V | 2 with R are modified from the naïve R 3 , the correct application of the random-walk argument, taking into account incompressibility, and the greater contribution of more distant structures to |L V | 2 , does result in the correct scalings. We therefore do not consider there to be an essential problem with the application of the randomwalk argument in our treatment of the Saffman helicity invariant in Section III B, or in our discussion of general invariants in Section IV B, though it should be understood that some of these scalings may have to be modified for particular conserved quantities that do not scale with R in the naïve way.

Appendix D: Numerical setup and simulation details
In this work, we have presented numerical simulations conducted with the spectral MHD code Snoopy [45]. The code solves the equations of incompressible MHD with hyper-viscosity and hyper-resistivity both of order n, viz., where p, the thermal pressure, is determined via the incompressibility condition In all our runs, Pm ≡ ν n /η n = 1. The code employs a pseudo-spectral algorithm in a periodic box of size 2π, with a 2/3 dealiasing rule. Snoopy performs time integration of non-dissipative terms using a low-storage, third-order, Runge-Kutta scheme, whereas dissipative terms are solved using an implicit method that preserves the overall third-order accuracy of the numerical scheme (comparisons between the results of Snoopy and other popular MHD codes for various nonlinear problems may be found in [95][96][97]; also see [98] for a test of Snoopy's Hall-MHD module). Units of time are chosen so that the initial magnetic energy is E M = 1/2 (i.e., the unit of time is the initial Alfvén crossing time of the box). We initialise the simulations with a magnetic field whose Fourier representation is where P ij (k) = δ ij − k i k j /k 2 is the projection operator perpendicular to k, G i (k) is the Fourier transform of a two-or three-dimensional Gaussian random field with zero correlation length. The parameter s controls the helicity -s = 1 for a helical field, s = 0 for a nonhelical field -and is related to the fractional helicity, σ, discussed in Section IV A by σ = 2s/(1 + s 2 ) [46]. The function F (k) sets the initial spectrum of the field: where D = 2, 3 is the number of spatial dimensions, a is the initial spectral exponent of the sub-inertial range (small k), and k c sets the initial peak of the magneticenergy spectrum, k p , via k p = (7/4) 1/2 k c . In all runs, we set a = 7, and k c = 25 =⇒ k p 33, so that the magnetic energy is initialised at scales 1/33 of the box size. We note that the initial k 7 sub-inertial-range spectrum is different from the k 4 spectrum (or k 3 in two dimensions) that the system quickly establishes in the subsequent evolution (see Figs. 2 and 8; see also [7] for a discussion of the role of the large-scale spectrum in decaying hydrodynamic turbulence). Our choice to initialise the simulation with this spectrum was motivated by a finding in our exploratory runs that the transient period before the system entered a period of power-law decay was shorter when the large-scale slope was allowed to establish itself organically. Presumably, this is because even if the spectral exponent is the right one, the structure of the synthetic field, Eq. (D4), is not, so it is better not to prejudice the system and let it decide for itself what structures to create at large scales (k < k p ).
We measure the decay exponents p M and p K by plotting |E i 1+1/pi dE i /dt| against time, selecting the parameter p i (i = M, K) so as to obtain a flat curve. As noted in [22], plots of this type give an unbiased estimate of the decay laws, as compared to more conventional logarithmic plots of E against t. The reason for this is that, because the power-law behaviour is only established after a short time t 0 following the initialization of the simulation, a plot of E ∼ (t − t 0 ) −p against t has a bias towards large energies, which decreases over time, giving the false impression of a steeper power law. Furthermore, a logarithmic t-axis exaggerates the importance of the initial times, during which the system has, in fact, not established a steady-state decay. In a similar way, we establish the value of α in Eq. (10) by plotting E α/2 M L against time and selecting the value of α to give a flat curve. The power laws and values of α thus obtained are given for all our runs in Table I   FIG. 14. Same as Fig. 13, but for non-helical simulations.