Cosmic-void observations reconciled with primordial magnetogenesis

It has been suggested that the weak magnetic field hosted by the intergalactic medium in cosmic voids could be a relic from the early Universe. However, accepted models of turbulent magnetohydrodynamic decay predict that the present-day strength of fields originally generated at the electroweak phase transition (EWPT) without parity violation would be too low to explain the observed scattering of γ-rays from TeV blazars. Here, we propose that the decay is mediated by magnetic reconnection and conserves the mean square fluctuation level of magnetic helicity. We find that the relic fields would be stronger by several orders of magnitude under this theory than was indicated by previous treatments, which restores the consistency of the EWPT-relic hypothesis with the observational constraints. Moreover, efficient EWPT magnetogenesis would produce relics at the strength required to resolve the Hubble tension via magnetic effects at recombination and seed galaxy-cluster fields close to their present-day strength.

It has been suggested that the weak magnetic field hosted by the intergalactic medium in cosmic voids could be a relic from the early Universe.However, accepted models of turbulent magnetohydrodynamic decay predict that the present-day strength of fields originally generated at the electroweak phase transition (EWPT) without parity violation would be too low to explain the observed scattering of γ-rays from TeV blazars.Here, we propose that the decay is mediated by magnetic reconnection and conserves the mean square fluctuation level of magnetic helicity.We find that the relic fields would be stronger by several orders of magnitude under this theory than was indicated by previous treatments, which restores the consistency of the EWPT-relic hypothesis with the observational constraints.Moreover, efficient EWPT magnetogenesis would produce relics at the strength required to resolve the Hubble tension via magnetic effects at recombination and seed galaxy-cluster fields close to their present-day strength.
Where might fields in voids come from?A popular idea (although not the only one, see [23]) is that they could be relics of primordial magnetic fields (PMFs) generated in the early Universe [24], including, prominently, at the electroweak phase transition (EWPT) [25].If so, the physics of the early Universe could be constrained by observations of the fields in voids -a remarkable possibility -provided the magnetohydrodynamic (MHD) decay of the PMFs between their genesis and the present day were understood.However, conventional theory of the decay ( [24]; see [13][14][15] for reviews) appears inconsistent with the EWPT-relic hypothesis: Ref. [26] argue that the lower bound (1) on B is too high to be consistent with PMFs generated at the EWPT without magnetic helicity (a topological quantity that quantifies the number of twists and linkages in the field, which is conserved even as energy decays [27]).Furthermore, they show that the amount of magnetic helicity required for consistency with Eq. ( 1) is greater than can be generated by baryon asymmetry at the EWPT, as estimated by Ref. [28].In principle, other mechanisms of magnetic-helicity generation may have been present in the early Universe; one idea is chiral MHD (see [29] and references therein).Whether enough net helicity can be generated via these mechanisms for PMFs to become maximally helical during their evolution remains an open question [30,31].
On the other hand, Ref. [26] note that their conclusions could be subject to modification by the contemporaneous discovery of "inverse transfer" of magnetic energy in simulations of non-helical MHD turbulence [32,33] (see [30,[34][35][36] for schemes for doing so based on decay laws obtained numerically).Recently, the inverse transfer was explained as a consequence of local fluctuations in the magnetic helicity, which are generically present even when the global helicity vanishes [37].In this paper, we demonstrate how this insight, together with the other key result of Ref. [37], and of Refs.[38][39][40], that the decay timescale is the one on which magnetic fields reconnect, restores consistency of the hypothesis of a nonhelical EWPT-generated PMF with Eq. (1).Intriguingly, we find that reasonably efficient magnetogenesis of non-helical magnetic field at the EWPT could produce relics with the ∼ 10 −11 G comoving strength that, it has been suggested, is sufficient to resolve the Hubble tension [41,42].Relics of this strength would also constitute seed fields for galaxy clusters that would not require much amplification by turbulent dynamo after structure formation to reach their observed present-day strength [43] (although dynamo would still be required to maintain cluster fields at present levels).

Results
We take the metric of the expanding Universe to be where a(t) is the scale factor, normalised to 1 at the present day, t is conformal time (related to cosmic time t by a(t)dt = dt), and x i are comoving coordinates.The expanding-Universe MHD equations can be transformed to those for a static Universe by a simple rescaling [44]: the scaled variables [where ρ, p, B, u, η and ν are the physical values of the total (matter + radiation) density, pressure, magnetic field, velocity, magnetic diffusivity and kinematic viscosity, respectively] evolve according to the MHD equations in Minkowski spacetime.As in previous work (see [13][14][15]), we consider the dynamics of the "tilded" variables in Minkowski spacetime and transform the result to the spacetime (2) of the expanding Universe via Eq.( 3).
Selective decay of small-scale structure.
Historically, it has been believed that statistically isotropic MHD turbulence decays while preserving the small-k asymptotic of the magnetic-energy spectrum E M (k) (see [13,14] and references therein).This idea, sometimes called "selective decay of small-scale structure", amounts to a statement of the invariance in time of the magnetic Loitsyansky integral, which, for isotropic turbulence without long-range spatial correlations, is related to E M (k) by Invariance of I L M implies In writing (6), we have assumed that the magneticenergy spectrum is sufficiently peaked around the energycontaining scale λ B for the latter to be equal to the correlation, or integral, scale of the field.This would not be the case for a scale-invariant magnetic field (often conjectured to be generated by inflationary mechanisms).We exclude such fields from our analysis in this paper, in which we consider "causal" fields -the sort that could be generated at a phase transition -exclusively.
Decay timescale.Eq. ( 6) can be translated into a decay law for magnetic energy by a suitable assumption about   4)], B and λB evolve along lines parallel to the ones shown in blue.The predicted values of modern-day B and λB are given by the intersection of these lines with Eq. (10).We see that even PMFs generated with ρB(t * ) ∼ ργ(t * ) and λB(t * ) ∼ rH (t * ) produce modern-day relics that are inconsistent with Eq. ( 1).
how the energy-decay timescale, depends on B, λ B and t.Regardless of this choice, Eqs. ( 6) and ( 7) have the following important property.Suppose that, after some intermediate time t c , τ (B, λ B , t) can be approximated by some particular product of powers of its arguments.Then, for all t τ (t c ), B2 decays as a power law: B2 ∝ t −p , where p is a number of order unity.Substituting this back into Eq.( 7), one finds which is an implicit equation for B = B(λ B ) that can be solved simultaneously with Eq. ( 6) for B(t) and λ B (t). Eq. ( 8) was first suggested by Ref. [24] on phenomenological grounds.Its great utility, which has perhaps not been spelled out explicitly, is that it implies that one need not know the functional form of τ ( B, λ B , t) during the early stages of the decay in order to compute B and λ B at later times.Thus, the effect of early-Universe physics (e.g., neutrino viscosity) on the decay dynamics can be safely neglected.
Inconsistency with observations.Assuming that the decay satisfies Eq. ( 6) and that its timescale is Alfvénic, viz., when it terminates at the recombination time t recomb [14] [Eq.(27) in Methods], Eq. ( 8) implies [24] B(t recomb ) ∼ 10 −8.5 G λ B (t recomb ) 1 Mpc (10) [see Eq. (31) in Methods].In (9), ρb is the baryon density, which appears because photons do not contribute to the fluid inertia at scale λ B at the time of recombination [45] [13,26].As is shown in Fig. 1, these values and Eq. ( 10) together lead to values of B and λ B at t recomb that violate the observational constraint (1).Note that λ B (t * ) ∼ 10 −2 r H (t * ) is, in fact, a more popular estimate, corresponding to the typical coalescence size of "bubbles of new phase" that form at the phase transition [46]; for this initial correlation scale, the predicted value of B is separated from the allowed values by around three orders of magnitude.A similar calculation led Ref. [26] to conclude that genesis of EGMFs at the EWPT was unlikely (although we note that significant modification of Eq. ( 1) by inclusion of the effects of plasma instabilities in the modelling of the electromagnetic cascade -see the comment below Eq. ( 1) -could alter this conclusion).
Saffman helicity invariant.We argue that the theory outlined above requires revision.First, the idea of "selective decay of small-scale structure" is flawed.This is because the kλ B 1 tail of the magnetic-energy spectrum E M (k) corresponds not to physical structures (as in the Richardson-cascade picture of inertial-range hydrodynamic turbulence) but to cumulative statistical properties of the structures of size λ B [47].Absent a physical principle to support the invariance of I L M (such as angular-momentum conservation for its hydrodynamic equivalent [47,48]), there is therefore no reason to suppose that the small-k asymptotic of E M (k) evolves on a longer timescale than the dynamical one of λ B -scale structures (if this is long compared to the magneticdiffusion timescale at scale λ B , then selective decay is valid, as the simulations of [24,49] confirm, but this is not the regime relevant to PMFs).
Instead, we propose that the decay of PMFs is con-trolled by a different integral invariant [37]: where h = Ã • B is the helicity density ( B = ∇ × Ã).Eq. ( 11) is equivalent to where H V is the total magnetic helicity contained within the control volume V .The invariance of I H can therefore be understood intuitively as expressing the conservation of the net mean square fluctuation level of magnetic helicity per unit volume that arises in any finite volume of non-helical MHD turbulence (see Fig. 2; we refer the reader concerned about the existence of such fluctuations to Section B of the Supplementary Information).Numerical evidence supporting the invariance of I H has been presented by Ref. [37] and independently by Refs.[50,51].From I H = const, we deduce We make two brief remarks.First, growth of I L M , and, therefore, the inverse-transfer effect discovered by Refs.[32,33,52], follows immediately from Eq. ( 13).This is because 5)] grows while B decays.Second, the value of the largescale spectral exponent does not affect the late-time limit of the decay laws in our theory (see Section C of the Supplementary Information), unlike in the "selective-decay" paradigm.Reconnection-controlled decay timescale.The second revision that we propose to the existing theory is that the field's decay timescale τ should be identified not with the Alfvénic timescale (9), but with the magneticreconnection one.This is because relaxation of stochastic magnetic fields via the generation of Alfvénic motions is prohibited by topological constraints, which can only be broken by reconnection.Refs.[37,40,50] have presented numerical evidence for a reconnection-controlled timescale for decays that occur with a dominance of magnetic over kinetic energy (see [38,39] for the same in 2D).Magnetically dominated conditions are relevant to the decay of PMFs because (i) the large neutrino and photon viscosities in the early Universe favour them, and (ii) once established, they are maintained, as reconnection is typically slow compared with the Alfvénic timescale.The identification of τ as the reconnection timescale implies that a number of different decay regimes are possible, as we now explain.
Under resistive-MHD theory, reconnecting structures in a fluid with large conductivity generate a hierarchy of current sheets at increasingly small scales via the plasmoid instability [53].The global reconnection timescale is the one associated with the smallest of these sheets (the "critical sheet"), which is short enough to be marginally stable ( [54,55], see [56] for a review).This timescale is where Pm = ν/η is the magnetic Prandtl number, which appears because viscosity can suppress the outflows that advect reconnected field away from the reconnection site, is the Lundquist number based on the reconnection outflow and S c ∼ 10 4 is the critical value of S for the onset of the plasmoid instability.Eq. ( 14) is a straightforward theoretical generalisation [56] to arbitrary Pm of a prediction for Pm = 1 [54] that has been confirmed numerically [55,57].Pm is given by Spitzer's theory [58] [Pm Sp ∼ 10 7 at recombination, see Eq. ( 37) in Methods] if the plasma is collisional, i.e., if the Larmor radius of protons r L = m i cv th,i /aeB is large compared to their mean free path, λ mfp (m i and v th,i ≡ 2T /m i are the mass and thermal speed of protons respectively).If, on the other hand, r L < λ mfp , which happens if B > B iso ≡ m i cv th,i /eaλ mfp , then the components of the Figure 3. Reconnection-controlled decay of non-helical PMFs.As in Fig. 1, purple regions denote values of B and λB excluded on physical or observational grounds [Eq.( 1)].Under decays that conserve IH [Eq.( 11)], B and λB evolve along lines parallel to the ones shown in blue.The predicted values of modern-day B and λB are given by the intersection of these lines with Eq. ( 8) evaluated at recombination [represented by lines (i-v), which are derived in Methods], with τ the prevailing decay timescale.The blue-gold line shows the locus of possible present-day states resulting from reconnection-controlled decays on the timescales explained in the main text, assuming that the microscopic viscosity of the primordial plasma was controlled by collisions between protons.The effective value of Pm in Eq. ( 14) might have been heavily suppressed when B > Biso if viscosity were then instead governed by plasma microinstabilities -the red-gold line shows the locus of modern-day states corresponding to the extreme choice of Pm 1 for B > Biso.In either case, we see that PMFs generated at the EWPT with a wide range of values of IH produce modern-day relics that are consistent with Eq. ( 1), and even with the stronger version of this constraint [see text below Eq. ( 1)] which is indicated by the pale purple region.
viscosity tensor perpendicular to the magnetic field are reduced by a factor (r L /λ mfp ) 2 , because protons' motions across B are inhibited by their Larmor gyration [59].These are the components that limit reconnection outflows because velocity gradients in reconnection sheets are perpendicular to the mean magnetic field.Therefore, Pm → (r L /λ mfp ) 2 Pm Sp = ( Biso / B) 2 Pm Sp in Eq. ( 14) if B > Biso ≡ a 2 B iso .
The validity of the resistive-MHD treatment that leads to Eq. ( 14) requires the fluid approximation to hold at the scale of the critical sheet: its width must be larger than either r L or the ion inertial length d i = m i c 2 /4πe 2 n i a 2 (n i is the proton number density) [54,60].If δ c < r L , d i , then the physics of the critical sheet is kinetic, not fluid, and the reconnection timescale is rather than (14).Eq. ( 17) is a robust numerical result whose theoretical explanation is an active research topic (see [61] for a recent study, [62,63] for reviews).We shall find in the next section that (17) is not the limiting timescale at recombination for almost any choice of initial condition consistent with EWPT magnetogenesis; our conclusions therefore do not depend sensitively on the validity of (17).
The decay timescale can also be limited by radiation drag due to photons [24]; this imparts a force −α ũ per unit density of fluid [see Eq. ( 54) in Methods].The drag is subdominant to magnetic tension at sufficiently small scales (as it does not depend on gradients of ũ), so does not contribute to Pm in Eq. ( 14).However, it can inhibit inflows to the reconnection layer.Balancing drag with magnetic tension at the integral scale λ B , we find an inflow speed ũ ∼ ṽ2 A /αλ B , so the timescale for magnetic flux to be processed by reconnection is The timescale for energy decay depends on whether largescale drag or small-scale reconnection physics is most restrictive: Comparison with observations.The locus of possible PMF states for different values of I H ∼ B4 λ 5 B under the theory that we have described is represented by the blue-gold line in Fig. 3.We denote the largest value of I H consistent with EWPT magnetogenesis by I H, max ; this corresponds to ρB (t * ) = ργ (t * ) and λ B (t * ) = r H (t * ).For I H 10 −29 I H, max , decays terminate on line (i) in Fig. 3 [Eq.(40) in Methods], which represents Eq. ( 8) with τ = τ rec given by Eq. ( 14) and Pm = Pm Sp .Use of Eq. ( 14) is valid here because δ c r L , d i [see Eqs. ( 41) and (42) in Methods].The Spitzer estimate of Pm is valid at recombination only if B Biso ∼ 10 −13 G [Eq. (44) in Methods], so decays with I H 10 −29 I H, max have a shorter timescale at recombination -they terminate on line (ii) [Eq.(45) in Methods], which represents Eq. ( 8) with τ = τ rec given by Eq. ( 14) and Pm ∼ (r L /λ mfp ) 2 Pm Sp .For I H 10 −2 I H, max , the states on line (ii) have δ c < d i , r L [see Eqs.(46) and (47) in Methods], so Eq. ( 14) is invalid for them.These decays pass through line (ii) at some time before recombination with timescale given by Eq. ( 17).However, they do access the domain of validity of Eq. ( 14) if, before t recomb , B becomes small enough for δ c to be comparable with relevant kinetic scales.When that happens, their timescale becomes much larger than t recomb so further decay is prohibited -these decays all terminate with B ∼ 10 −11 G, which corresponds to δ c ∼ d i at t recomb [see Eq. (46) in Methods].Decays with I H 10 8 I H, max are radiationdrag limited at recombination [line (iv); Eq. ( 55) in Methods] -such decays are inconsistent with EWPT magnetogenesis, but could originate from magnetogenesis at the quantum-chromodynamic (QCD) phase transition, when r H ∼ 10 −6 Mpc [13,26].( The relic of a field with λ B (t * ) ∼ 10 −2 r H (t * ) ∼ 10 −10 Mpc at the EWPT would therefore be consistent with Eq. ( 1) -modulo any modifications for plasma instabilities in voids [17][18][19][20][21][22] -if ρB (t * ) 10 −6.5 ργ (t * ).This confirms the assertion in the title of this paper.Intriguingly, if instead ρB (t * ) ∼ ργ (t * ) and λ B (t * ) 10 −2 r H (t * ), then we find B ∼ 10 −11 G at recombination.PMFs of this strength would provide a seed for magnetic fields in galaxy clusters that would not require significant amplification by turbulent dynamo after structure formation to reach their present day strength of ∼ µG [43], although dynamo would still be required to maintain cluster fields at present levels.We emphasise that a cluster field so maintained by dynamo need not (and, in all likelihood, would not) retain memory of its primordial seed.We also note that PMFs of 10 −11 G strength are considered a promising candidate to resolve the Hubble tension, by modifying the local rate of recombination [41,42].
As an aside, we note that the relevance of reconnection physics is not restricted to non-helical decay [37].Some analogues for maximally helical PMFs of the results of this section (relevant for magnetogenesis mechanisms capable of parity violation) are presented in Section A of the Supplementary Information.Role of plasma microinstabilities.Finally, we note that, for B > Biso , the effective values of ν and η might be dictated by plasma "microinstabilities" rather than by collisions between protons [64] (this is conjectured to happen in galaxy clusters [65]).In Methods, we show that the decay of the integral-scale magnetic energy is too slow to excite the "firehose" instability that is important in the cluster context [see Eq. ( 60)].Nonetheless, we cannot rule out other microinstabilities -for example, the excitation of the "mirror" instability by reconnection has been studied recently by Ref. [66], although its effect on the rate of reconnection remains unclear.The most dramatic effect that microinstabilities in general could plausibly have would be to reduce the effective value of Pm to 1 if B > Biso (see [67,68]).This corresponds to the red-gold line in Fig. 3, which remains consistent with Eq. ( 1) for I H 10 −20 I H, max .Compatibility between the EWPT-magnetogenesis scenario and the observational constraints on EGMFs therefore appears robust.

METHODS
Post-recombination evolution.In the matter-dominated Universe after recombination, the transformation that maps Minkowski-spacetime MHD onto its expanding-Universe equivalent is not Eq.( 3), but [24] ρ = a 3 ρ, p = a 4 p, As a ∝ t 2 in the matter-dominated Universe, t ∝ log t, so a power-law decay in rescaled variables corresponds to only a logarithmic decay in comoving variables [14].Thus, in computing the expected present-day strength of EGMFs, one may assume the decay of B to terminate at recombination with negligible error.
Derivation of Eq. (10).In order to apply Eq. ( 8), we require an expression for the conformal time at recombination, t recomb .From the Friedmann equation, where G is the gravitational constant, the "entropy equation" where g is the number of degrees of freedom of the radiation field and T is the temperature, and Stefan's law for the radiation density where χ = π 2 /90c 5 3 (we work in "energy units" for temperature, with Boltzmann constant kB = 1), it can be shown that dT dt where the subscript 0 refers to quantities evaluated at the present day.Because (g/g0) 1/6 1, one may solve Eq. ( 25) to give an expression for the cosmic temperature as a function of conformal time, With g0 = 2 (for the two photon-polarisation states), one obtains t ∼ 10 16.5 s T 0.3 eV Therefore, Eq. ( 8) becomes τ ∼ 10 16.5 s T 0.3 eV Thus, t recomb ∼ 10 16.5 s.Eq. ( 28) can be used to relate B and λB under the assumption that the decay occurs on the Alfvénic timescale τ ∼ λB/ṽA [Eq.( 9)].As noted in the main text, ṽA should be computed using the baryon density ρb , because the photon mean free path [13] (where σT is the Thompson-scattering cross-section) is large compared with λB at the time of recombination, indicating that photons are not strongly coupled to the fluid [45].However, because ρb ργ at the time of recombination, the decoupling of photons does not affect Eq.(10).The Alfvén speed is where we have used ρb = a 4 ρ b a 4 min b , with mi the proton mass and n b the WMAP value for the baryon number density n b 2.5 × 10 −7 cm −3 a −3 [69], and taken a T0/T [Eq.( 23)].
Comparing Eq. ( 9) and Eq. ( 28), and substituting Eq. ( 30), we have Evaluated at T = T (t recomb ) = 0.3 eV, this is Eq. ( 10).Derivation of line (i) of Fig. 3. Line (i) represents Eq. ( 14) evaluated at the time of recombination t recomb , with Pm = PmSp ≡ νSp/ηSp, where νSp and ηSp are the comoving Spitzer values of kinematic viscosity and magnetic diffusivity respectively [58].We first evaluate PmSp.Under Spitzer theory, the dominant component of the plasma viscosity at the scale of the rate-determining current sheet is due to ion-ion (i.e., proton-proton) collisions.The collision frequency is [58] νii ∼ e 4 ni ln Λii where e is the elementary charge, ni the ion number density, mi the ion mass, Ti the ion temperature, and ln Λii the Coulomb logarithm for ion-ion collisions.Neglecting any anisotropising effect of the magnetic field (see main text), the comoving isotropic kinematic viscosity is [70] νSp where v th,i = 2Ti/mi is the thermal speed of ions, and we have assumed Ti T , used a T0/T [Eq.( 23)], taken ni to be equal to the WMAP value for the baryon number density n b 2.5 × 10 −7 cm −3 a −3 [69], and estimated the Coulomb logarithm ln Λii by ln Λii ln Similarly, the electron-ion collision frequency is [70] νei ∼ e 4 ne ln Λei where ne ni is the electron number density, Te the electron temperature, and ln Λei the Coulomb logarithm for electronion collisions.Eq. ( 35) leads to the Spitzer [58] value for the magnetic diffusivity ηSp ∼ νeimec 2 4πnee 2 a ∼ 10 10.5 cm 2 s −1 T 0.3 eV where we have used ln Λei ln Λii 20, assumed the electron temperature Te T , and again neglected any anisotropy resulting from the magnetic field.From Eqs. ( 33) and ( 36 (37) Let us now evaluate the Lundquist number, Eq. ( 15), in order to compare it with Sc, as Eq. ( 14) requires.Note that, as above, it is the Alfvén speed based on baryon inertia that appears in Eq. ( 15); photons are even more weakly coupled to the cosmic fluid at reconnection scales than at scale λB as the former are typically small compared with the latter.Using Eqs. ( 13), (30) Eq. (38) shows that S Sc ∼ 10 4 [unless B(t * ) or λB(t * ) are very small, in which case their evolution is inconsistent with the observational constraint (1), so we neglect this possibility for simplicity].Substituting Eq. ( 37), we find that the decay timescale ( 14) is τ ∼ 10 5.5 T 0.3 eV Comparing Eqs. ( 28) and ( 39), and again substituting Eq. ( 30), we find Evaluated at T = T (t recomb ) = 0.3 eV, this is line (i) of Fig. 3. Finally, we note that when reconnection occurs under large-Pm conditions with isotropic Spitzer viscosity, the ratio of δc [Eq.( 16)] to rL [defined below Eq. ( 15)] prior to recombination is independent of the magnetic-field strength, temperature and density: where we have used Eqs (36), (37) and (30).Thus, δc > rL always.Furthermore, we find from Eqs. ( 15), ( 16), ( 30), ( 36), (37) and the definition of di [see below Eq. ( 16)] that Therefore, δc > di, rL at recombination for all relevant field strengths, so we are justified in using fluid theory to describe decays with B < Biso [evaluated in Eq. ( 44)].
As described in the main text, Eq. ( 40) is valid when B is small enough for the Larmor radius of ions rL to be larger than their mean free path The critical magnetic field strength above which this condition is no longer satisfied is Derivation of line (ii) of Fig. 3. Line (ii) represents Eq. ( 14) evaluated at the time of recombination t recomb , with magnetic Prandtl number Pm ∼ (rL/λ mfp ) 2 PmSp = ( Biso/ B) 2 PmSp.Note that this suppression of Pm relative to PmSp increases the value of S at any given ṽA and λB relative to the value (38) of S that corresponds to Pm = PmSp.We therefore expect this family of decays also to have S Sc ∼ 10 4 .The inclusion of the factor of ( Biso/ B) 2 in Pm modifies Eq. ( 40) straightforwardly: it becomes Evaluated at T = T (t recomb ) = 0.3 eV, this is line (iv) of Fig. 3.
The analogue of Eq. ( 42) for Pm ∼ ( Biso/ B) 2 PmSp is while the corresponding analogue of Eq. ( 41) is Eq. ( 46) shows that δc di at t recomb if B 10 −11 G, while Eq. ( 47) indicates that δc rL if B 10 −12 G. Following the prescription described in [54], we use the former condition on B as the domain of validity of Eq. ( 14) in Fig. 3, though we note that our results do not depend strongly on this choice -the order-of-magnitude difference between the two critical values of B is comparable to the degree of accuracy to which our scaling arguments are valid.
We also note that the temperature dependence of Eq. ( 46) means that a decaying field that developed δc di before recombination would have done so at a field strength B < 10 −11 G; strictly, therefore, the decay of primordial fields should terminate somewhere below the horizontal part of the blue-gold curve in Fig. 3, not directly on it.However, the difference is order unity and thus negligible for the purposes of our order-of-magnitude estimates.This is because magnetic decay was strongly suppressed by radiative drag at early times [a consequence of the strong temperature dependence of Eq. ( 55)] -i.e., when temperatures exceeded around 10 2 × 0.3 eV.For all relevant values of IH , the magnetic-field strength would therefore have greatly exceeded the critical value required for δc ∼ di until the time that corresponds to this temperature, and by that time the critical field strength indicated by Eq. ( 46) was already within a small factor of its value at recombination.Derivation of line (iii) of Fig. 3. Line (iii) represents Eq. ( 14) at the time of recombination t recomb , with Pm 1.With Pm 1, Eq. ( 38) should be replaced by S ∼ 10 Comparing Eqs. ( 28) and ( 39), and substituting Eq. ( 30), we find Evaluated at T = T (t recomb ) = 0.3 eV, this is line (iii) of Fig. 3.The analogues of Eqs. ( 42) and ( 41) for Pm 1 (but η ∼ ηSp) are and Note that the field strength at which δc ∼ di is approximately equal to Biso at recombination (both are ∼ 10 −13 G), while δc rL.The red-gold line in Fig. 3 therefore extends past line (iii) to line (iv) along the line B ∼ Biso.Radiation drag and the derivation of line (iv) of Fig. 3.As well as by viscosity arising from collisions between ions, the kinetic energy of primordial-plasma flows (after neutrino decoupling) can be dissipated by electron-photon collisions (Thompson scattering).Around the time of recombination, the comoving mean free path of photons, Eq. ( 29), is much larger than the anticipated correlation scale of the magnetic field (and, therefore, of any magnetically driven flows).Under these conditions, the effect of Thompson scattering is to induce a drag on electrons.Owing to the collisional coupling between ions and electrons, this drag can dissipate bulk plasma flows.
The comoving drag force on the fluid per unit baryon density is where [24] α As explained in the main text, the effect of drag is most important at the scale λB (it becomes increasingly subdominant to magnetic tension at smaller scales) where it inhibits inflows to the reconnection layer.When the timescale τα ≡ αλ 2 B /ṽ 2 A on which flux can be delivered to the layer by strongly dragged inflows is larger than the reconnection timescale of the critical sheet τrec [see Eq. ( 19)], τα gives the timescale for energy decay.Eq. ( 28) with τ = τα yields, after substitution of Eqs.(30) and Eq. ( 54) Evaluated at T = T (t recomb ) = 0.3 eV, this is line (iv) of Fig. 3.
Non-excitation of the firehose instability.Plasma with an anisotropic viscosity tensor can, in principle, be unstable to a variety of instabilities that develop at kinetic scales.For a decaying magnetic field, an instability of particular importance is the "firehose", which can generate the growth of small-scale magnetic fields in response to the decay of largescale ones [64,71].This happens if the size of the (negative) pressure anisotropy ∆ exceeds a critical value: where p and p ⊥ are the thermal pressures parallel and perpendicular to the magnetic field, and is the "plasma beta".∆ can be estimated as [64] ∆ ∼ 1 νii where t is cosmic time [defined below Eq. ( 2)].Naturally, the value of βi at any given T depends on the evolution of the magnetic field.A lower bound on the value of B at any given time for a given initial condition is the one that would develop from a decay on the kinetic reconnection timescale, τ ∼ 10λB/ṽA [Eq.( 17)].Solving Eqs. ( 13), ( 17), ( 28) and ( 30) simultaneously, we find that this is Eq. ( 60) suggests that the threshold for instability (56) is never met, unless λB(t * ) and/or B(t * ) are so small as to be inconsistent with the observational constraint (1).Numerical simulation.The numerical simulations visualised in Fig. 2 and described in the Supplementary information were conducted using the spectral MHD code Snoopy [72].The code solves the equations of incompressible MHD in Minkowski spacetime with hyper-viscosity and hyper-resistivity both of order n, viz., where p, the thermal pressure, is determined via the incom-pressibility condition The code uses 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, thirdorder, Runge-Kutta scheme, whereas dissipative terms are solved using an implicit method that preserves the overall third-order accuracy of the numerical scheme.In all runs presented here, we employ νn = ηn = 10 −12 , n = 6 and use a resolution of 512 3 .
For completeness, here we provide the results for maximally helical fields that correspond to those presented in the main text for non-helical fields; these results are relevant for magnetogenesis mechanisms that are capable of parity violation.The decay of such fields conserves the net magnetic helicity, resulting in the self-similar scaling As in the non-helical case, the decay proceeds on reconnection timescales [37]; the possible decay regimes are the same as those described in the main text.Under Eq. ( 64), the PMF evolution in the ( B, λ B ) plane is parallel to Eq. ( 1) [see Fig. 4].Thus, any field satisfying will satisfy the observational constraint (1) at recombination, as is well known (see, e.g., [13,14]).The locus of present-day PMF states for decays that occur on the reconnection timescales explained in the main text is shown by the blue-gold line in Fig. (4).Analogously to I L, max and I H, max in the main text, we denote the largest value of the mean magnetic-helicity density h that is consistent with EWPT magnetogenesis by h max [this corresponds to ρB (t * ) = ργ (t * ) and λ B (t * ) = r H (t * )].For h 10 −15 h max , the decay of PMFs terminates on line (i) in Fig. 4, which corresponds to Eq. ( 14) of the main text with Pm ∼ Pm Sp [Eq.(37)].For 10 −15 h max h 10 −11 h max , the decay of PMFs terminates on line (ii) [Eq.(45) in Methods], which corresponds to Eq. ( 14) of the main text with Pm = (r L /λ mfp ) 2 Pm Sp .For 10 −7 h max h 10 −5 h max , decays terminate at B ∼ 10 −11 , G, which corresponds to δ c ∼ λ mfp , as explained in the main text.Finally, decays are radiation-drag limited for h 10 −5 h max , and therefore terminate on line (iv) [Eq.(55) in Methods].We note that, for h 10 −5 h max , the role of magnetic reconnection in determining the decay timescale implies significantly stronger relic fields than would be expected under the decay physics envisaged by [24], i.e., Alfvénic [Eq.(9); line (v)] or radiation-drag-limited [line (iv)] decay.
As in the main text, we also indicate by a red-gold line the locus of present-day PMF states if Pm 1 (due to plasma microinstabilities) for B Biso .
B. Decay of non-helical magnetic fields with IH = 0 As explained in the main text, the invariance of I H follows from the conservation of the fluctuation level of magnetic helicity.While we view fluctuations in magnetic helicity to be a generic feature of real MHD turbulence1 , it is nonetheless possible to construct artificial field configurations for which the helicity of each magnetic structure vanishes -this will be the case if they have no twists and do not interlink.Strictly, therefore, the possibility that PMFs might have been generated without helicity fluctuations cannot be ruled out.
A priori, it appears that this kind of field might relax in a fundamentally different manner to the one described in the main text.This was the view that we expressed in Ref. [37]: there, we suggested that fields with I H = 0 might decay subject to the conservation of invariants associated with the velocity, rather than the magnetic, field.This is because individually non-helical structures (unlike helical ones) can relax under entirely flux-frozen dynamics, by driving flows with ũ ∼ B (a process sometimes called kinetic diffusion [13]).Plausibly, the decay of those flows would respect the invariance of the hydrodynamic Loitsyansky integral, which encodes the conservation of angular momentum L = x × u [48] (in the same fluctuating manner as I H encodes helicity conservation [37]). 2 Denoting the characteristic size and scale of the velocity field by ũ and λ u respectively, I L ∼ ũ2 λ 5 u .Conservation of I L therefore implies ũ2 λ 5 u ∼ const.This suggests that B2 λ 5 B ∼ const also, if ũ ∼ B and λ B ∼ λ u , which seems reasonable for, e.g., a magnetic field maintained by the dynamo effect3 .This returns us to Eq. ( 6), i.e., to the same prediction that was shown to be inconsistent with the observational constraints by Ref. [26].
On the other hand, if the magnetic field were maintained by dynamo, then it seems unlikely that I H = 0 would be maintained.This is because random helicity fluctuations could be generated freely at resistive scales (as the Lundquist number is order unity there), where dynamo primarily generates magnetic field (at least in its kinematic stage) [79].Thus, I H could become non-zero, although it would not need to be conserved if magnetic energy remained concentrated at the resistive scales.If, however, the dynamo-replenished magnetic fields later transferred to larger scales and saturated with λ B ∼ λ u , as supposed above, while still having helicity fluctuations, then I H would become invariant, because the integral scale of the magnetic field would be much larger than the resistive scale.This would push us back to the scaling I H ∼ B4 λ 5 B ∼ const [Eq.(13)].Moreover, we conjecture that the size of the conserved product B4 λ 5 B would be of the same order as its value for the initial field with I H = 0, because memory of its B and λ B would be retained by the velocity field.In Fig. 5, we present results from a numerical simulation (Simulation A) designed to assess these arguments.We initialise a large number of untwisted, non-interlinking magnetic-flux loops (otherwise distributed in a random, statistically isotropic way) in a periodic simulation domain [see Fig. 5(a)].At t = 0, I H = 0 because the loops each have zero magnetic helicity.For the purpose of comparison, we also present a second simulation (Simulation B) with the same setup but without the non-interlinking condition, instead starting with many loops superimposed on top of each other.This field has complex initial topology, as Fig. 5(b) indicates, but no net helicity, in the sense that H V /V B2 λ B for all V .However, unlike the field in Simulation A, it has I H = 0, because superimposing loops creates linkages in the magnetic field.Figs.5(c-d) show the evolution of energy spectra for the two simulations.In Simulation A, unlike in Simulation B, there is an immediate and rapid decay of the magnetic energy as the loops contract and drive flows.There is a corresponding decay of the large-scale (low-k) spectral tail, which demonstrates the non-invariance of I L M .The newly generated kinetic energy is comparable in magnitude to the initial magnetic energy [see Fig. 5(i)], and its spectrum peaks close to the initial peak of the magnetic-energy spectrum. 4On the other hand, the contraction of the loops leaves magnetic energy concentrated at small (resistive) scales, where it can be refuelled by the dynamo effect associated with the newly generated flows.As we anticipated above, this resistive-scale magnetic field has random fluctuations in magnetic helicity: this is shown explicitly in Fig. 5(e) where, for volumes V taken to be spheres of radius R, we plot H 2 V /V vs. R at regular intervals in time (the average is taken over a large sample of spheres with centres throughout the simulation box).While H 2 V ∝ V 2/3 at t = 0 (not shown) because H V is dominated by random surface contributions at this time, this scaling is replaced by H 2 V ∝ V as soon as turbulence develops, indicating I H = 0 [see Eq. ( 12) of the main text].Though I H , which is the value of H 2 V /V in the flat part of the curves in Figs.5(e-f), decays by around an order of magnitude during the first few eddy-turnover times, Fig. 5(e) shows that its decay ceases after that.This is consistent with our suggestion above that I H should become constant when the dynamo saturates, due to migration of the helicity-containing scale towards the flow scale λ u . 5This interpretation is supported by the evolution of the helicity-variance spectrum Θ(k) [see Fig. 5(g)], which encodes the characteristic size of helicity fluctuations at each scale. 6In Simulation A, Θ(k) is concentrated around the dissipation scale at early times (though after the decay of the magnetic loops), but later moves to larger scales.I H [which is proportional to the coefficient of k 2 in the Θ(k → 0) ∝ k 2 asymptotic [37]] ceases to decay once the peak of Θ(k) is moderately separated from the dissipation scales.
The value of I H ultimately attained by the magnetic field in Simulation A is smaller than the one in Simulation B by a factor of around 10 4 .This appears to contradict our conjecture that dynamo should generate I H of the same size as B4 λ 5 B at the initial time.On the other hand, we note that (i) this factor may well be smaller for a simulation at larger resolution and larger Prandtl number (recent work has shown that extremely large resolutions are required to probe the asymptotic nature of the large-Pm dynamo [81]), and that (ii) the strong scaling of I H with B and λ B means that even a factor-10 4 reduction in I H corresponds only to a factor-10 reduction in B or λ B .This means that a PMF generated with I H = 0 would migrate only a relatively short distance on the ( B, λ B ) plane (Fig. 3) before settling to decay with B4 λ 5 B ∼ const.To summarise, there appear to be both theoretical and numerical reasons to believe that a PMF generated with I H (t * ) = 0 at the initial time t * would, via an initial period of rapid decay and subsequent regeneration via dynamo, develop I H ∼ B(t * ) 4 λ B (t * ) 5 ∼ const.At later times, a magnetically dominated state would likely be re-established because the flows will drive Alfvénic turbulence, which cascades to small scales and is dissipated by viscosity (which may be large, if associated with neutrinos or photons), while background "quasi-force-free" magnetic fields persist, decaying only on the magnetic-reconnection timescale, as described in the main text.There is some evidence of this in Fig. 5(i), which shows that magnetic energy becomes larger than kinetic in Simulation A at late times.
C. The effect of the large-scale spectral slope: coexistence of flux and helicity invariants In this paper, we have contrasted our theory of I H -conserving PMF decay with the previously accepted theory based on "selective decay of small-scale structure", i.e., the invariance of the large-scale asymptotic of the magnetic-energy spectrum.One success of our theory is that it explains the inverse-transfer effect observed in simulations of magnetic fields initialised with E M (k → 0) ∝ k 4 ( [32,33]; see main text); this effect is manifestly not compatible with selective decay.On the other hand, Ref. [49] observe that inverse transfer is not present in simulations that are initialised with sufficiently shallow large-scale spectra [namely, with E M (k → 0, t = 0) ∝ k n , where n < 3].Instead, they find that the k → 0 asymptotic of E M (k) is preserved.This result raises questions of whether a "selective-decay-like" principle might be at work in such decays, and what its effect might be on the laws for the decay of energy and growth of the integral scale.In this Section, we explain the invariance of this k n asymptotic as a consequence of the conservation of magnetic flux, but also argue that, beyond an initial transient, flux conservation does not affect the decay laws if n > 3/2 (as is the case in all models of EWPT magnetogenesis of which we are aware).It is therefore not necessary to know the precise value of n to compute the present-day properties of EGMFs under the relic-field hypothesisthe theory presented in the main text is valid independently of it.
1. Invariance of the large-scale spectral asymptotic for n ≤ 3 In general, the large-scale spectral asymptotic is frozen in time when the coefficient of k n in E M (k → 0) is proportional to some statistical invariant.As explained in the main text, this is not the case when correlations in B decay rapidly with distance, because then E M (k → 0) ∝ I L M k 4 [Eq.( 4)] where I L M = const.However, for n ≤ 3, it turns out that the coefficient of k n is proportional to an invariant that is related to the conservation of magnetic flux.Physically, this invariant encodes the fact that, over sufficiently large volumes, local fluctuations in magnetic flux may sum to a non-zero net fluctuation level, which must be conserved as the field decays.Spatial correlations must be long (and hence spectra must be shallow) for the fluctuation level to be non-zero, because ∇ • B = 0 means that magnetic structures without sufficiently strong far-field components have net zero flux.The relevant measure of correlation strength is the large-r asymptotic of the magnetic field's longitudinal correlation function, χ B (r) ≡ B r (x)B r (x + r) / B 2 r , where B r = B • r/r.The argument is particularly transparent if χ B (r → ∞) ∝ r −3 (as, for example, would be the case for a superposition of many randomly positioned and oriented magnetic dipoles), as then it can be shown that where is the Saffman flux invariant [37].The invariance of I B encodes conservation of the fluctuation level of magnetic flux in the same manner as the invariance of I H does for magnetic helicity.More generally, if E M (k → 0) = Ck n with spectrum, which on dimensional grounds is of size ∼ B2 λ 5 B k 4 : From B4 λ 5 B ∼ const, we have a decreasing function of time.
The evolution of the magnetic-energy spectrum depicted in Fig. 6 is manifestly non-self-similar.As a result, transient order-unity changes in B4 λ 5 B and the decay timescale as a function of B and λ B should be expected at early times.On the other hand, the decay does become approximately self-similar at late times, when k c 1/λ B , so any deviation from the theory proposed in the main text becomes small as t becomes large.Finally, we acknowledge that, while we expect the evolution depicted in Fig. 6 to be valid for any initial spectrum with n > 3/2, Ref. [49] do not observe the formation of a "k 4 bulge" in their simulation with n = 2.We believe this to be a result of insufficient scale separation in that simulation: with n = 2, Eq. ( 77) implies k c λ B ∼ [λ B /λ B (0)] −1/4 , so with λ B /λ B (0) 10 (see Fig. 16 of Ref. [49]), k c λ B 0.6.It is therefore not surprising that these scales cannot be distinguished.

Figure 2 .
Figure 2. Slice of magnetic-helicity density from a simulation of decaying non-helical MHD turbulence.The turbulence breaks up into patches of positive and negative helicity h (computed in the Coulomb gauge; ∇ • Ã = 0), shown in red and blue, respectively.The invariance of IH (see main text) is a manifestation of the conservation of the net magnetic-helicity fluctuation level arising in large volumes.Because of the complex magnetic-field topology, the rate-setting process for the decay is magnetic reconnection: reconnection sites, indicated in the figure by patches of large current density |∇ × B| (black; variable opacity scale), typically form between the helical structures.See Methods for details of the numerical setup.

9 . ( 59 )
Using this lower bound on B, we can obtain an upper limit on |βi∆|

Figure 4 .
Figure 4. Evolution of a maximally helical PMF.Analogue of Fig. 3, showing the decay of a maximally helical magnetic field generated at the EWPT.
(29) Eq.(29)in Methods].An approximate upper bound, I L M , max , on I L M follows from assuming that the magnetic-energy density ρB ≡ B2 /8π and the electromagnetic-radiation density ργ were equal at the time t * of the EWPT while λ B (t * ) was equal to the Hubble radius r H (t * ).This corresponds to B(t * ) ∼ 10 −5.5 G and λ B (t * ) ∼ r H (t * ) ∼ 10 −10 Mpc (37)d(37), we find the Lundquist number