Evolution of primordial magnetic fields in mean-field approximation

We study the evolution of phase-transition-generated cosmic magnetic fields coupled to the primeval cosmic plasma in turbulent and viscous free-streaming regimes. The evolution laws for the magnetic energy density and correlation length, both in helical and non-helical cases, are found by solving the autoinduction and Navier-Stokes equations in mean-field approximation. Analytical results are derived in Minkowski spacetime and then extended to the case of a Friedmann universe with zero spatial curvature, both in radiation and matter dominated eras. The three possible viscous free-streaming phases are characterized by a drag term in the Navier-Stokes equation which depends on the free-streaming properties of neutrinos, photons, or hydrogen atoms, respectively. In the case of non-helical magnetic fields, the magnetic intensity $B$ and the magnetic correlation length $\xi_B$ evolve asymptotically with the temperature $T$ as $B(T) \simeq \kappa_B (N_i v_i)^{\varrho_1} (T/T_i)^{\varrho_2}$ and $\xi_B(T) \simeq \kappa_\xi (N_i v_i)^{\varrho_3} (T/T_i)^{\varrho_4}$. Here, $T_i$, $N_i$, and $v_i$ are, respectively, the temperature, the number of magnetic domains per horizon length, and the bulk velocity at the onset of the particular regime. The coefficients $\kappa_B$, $\kappa_\xi$, $\varrho_1$, $\varrho_2$, $\varrho_3$, and $\varrho_4$, depend on the index of the assumed initial power-law magnetic spectrum, $p$, and on the particular regime, with the order-one constants $\kappa_B$ and $\kappa_\xi$ depending also on the cut-off adopted for the initial magnetic spectrum. In the helical case, the quasi-conservation of the magnetic helicity implies, apart from logarithmic corrections and a factor proportional to the initial fractional helicity, power-like evolution laws equal to those in the non-helical case, but with $p$ equal to zero.


I. INTRODUCTION
All observed galaxies and cluster of galaxies present microgauss, large-scale magnetic fields. The origin of these cosmic magnetic fields has not yet been fully understood, although several proposals, especially in the last years, have been put forward to explain why our universe is magnetized (for reviews on cosmic magnetic fields, see [1][2][3][4][5]). These proposals belong to two distinct classes of generating mechanisms characterized, roughly speaking, by the time when they operate. Astrophysical mechanisms work on at the epoch of large-scale structure formation or later, and can be based on the Biermann battery effect [6] in the first supernova remnants [7], or on the Weibel instability [8] in intergalactic plasmas [9]. Mechanisms performing in the early universe can, on they part, classified in mechanisms operating in the very early universe, namely during an inflationary epoch of the universe [10], and mechanisms operating after inflation [11][12][13][14], for example during electroweak (QED) or quark-hadron (QCD) cosmological phase transitions [15].
Magnetic fields generated during large-scale structure formation through astrophysical processes usually have large correlation scales but low intensities. Then, the precondition for explaining the observed fields is an amplification of these seeds. In principle, dynamo mechanisms operating in gravitationally bound large-scale structures (such as galaxies and clusters of galaxies) could substantially increase the intensity of such fields up to observed values [16]. However, astrophysical mechanisms cannot explain the presence of large-scale magnetic fields recently detected in cosmic voids [17].
Inflation-produced magnetic fields have the potentiality to explain the origin of cosmic magnetic fields since they possess a large correlation. The unpleasant feature of this class of mechanisms is that, in order to obtain astrophysically interesting magnetic intensities, one has to introduce some nonstandard (to wit, speculative) interaction term in the photon field Lagrangian. A possible exception is represented by a generating mechanism first noticed in [18] and successively analyzed in the series of papers [19]. Here it is shown that, starting from the standard Maxwell Lagrangian and assuming a small negative curvature term (not present in all the other inflationary mechanisms), a very strong magnetic field can emerge as the result of photons creation from a varying gravitational field. The validity of the results in theses works, however, has been put into question in [20] and [21]. According to the authors of [21] indeed, the intensity of the produced field would be today, in the best case scenario, as low as B 0 ∼ 10 −59 G and thus astrophysically unimportant. Another exception (the last one to our knowledge) of the use of nonstandard physics for the case of inflationary magnetic fields is represented by the recent work [22] (see also [23]). Here, criticisms to all previous inflationary generating mechanisms for cosmic magnetic fields have been drawn. The main criticism concerns the use of unrenormalized, and then unphysical, vacuum magnetic fluctuations instead of the renormalized ones (the only ones which are physically acceptable). In [22], it is shown that in standard Maxwell theory, the actual magnetic field resulting from inflationary renormalized vacuum magnetic fluctuations is scale independent and has an intensity which depends only on the scale of inflation. If this scale is around 10 16 GeV (a plausible value for such a scale), that field directly accounts for galactic and galaxy cluster magnetic fields.
Magnetic fields created during cosmological phase transitions have small length scales. This is due to the fact that microphysical processes which participate in the generation of the magnetic field are necessarily uncorrelated on scales greater than the Hubble scale which represents an (event) horizon for all causal physical processes. Nevertheless, the primordial universe could have been very turbulent at the time of QED and QCD phase transitions. This means that any magnetic field coupled to the primeval plasma could have been processed by magnetohydrodynamic (MHD) effects. Indeed, it is well-known in the literature that turbulence effects increase the typical magnetic length scale [24]. Moreover, if the magnetic field is helical the growth of magnetic correlation could be much more significant. It is interesting to observe that if the large-scale magnetic fields we observe today were relic helical fields from cosmological phase transitions (or inflation), then this would be a manifestation of a macroscopic and primeval P and CP violation, since magnetic helicity is odd under discrete P and CP transformations. 1 Moreover, models for generating helical magnetic fields in the early universe there exist in the literature [25,26].
It is worth noting that magnetic fields produced in the early universe are subject to a variety of constraints coming from the fact that the presence of a primordial magnetic field could spoil the predictions of the standard cosmological model. The most important limits came from the study of Big Bang Nucleosynthesis [27] and Cosmic Microwave Background (CMB) anisotropies [28][29][30]. In particular, large-scale, primordial helical magnetic fields could leave peculiar signatures in CMB radiation [31], since they would introduce a parity-odd cross-correlations of the CMB anisotropies, not present in the case of a non-helical magnetic field.
Also, primordial helical magnetic fields could affect the phenomenology of axions since the latter are coupled, through the axion-photon interaction term, to external magnetic fields. Indeed, the generation of helical magnetic fields before the axion coherent oscillations start at a temperature around f ew GeV, could be in contradiction with the existence of the axion [32], since a strong helicity production would in turn produce too much of an axion relic abundance, in disagreement with astrophysical observations. Finally, it is interesting to observe that the presence of large-scale helical magnetic fields in the actual universe could be directly detected by analyzing the propagation properties of charged cosmic rays, as shown in [33].
Magnetic fields generated during inflation or cosmological phase transitions can evolve as the universe expands, since their properties can be modified by magnetohydrodynamic turbulent effects due to the coupling with the primeval cosmic plasma. As recently shown in [34], typical inflation-produced magnetic fields remain almost unchanged on scales of cosmological interest even in the presence of a turbulent plasma (as in the case, for example, of a cosmic phase transition) although on small scales, which are however astrophysically uninteresting, their power gets gradually suppressed.
The fact that the universe during a cosmological phase transition could be very turbulent was first realized in [35], where the importance of MHD turbulent effects on the evolution of phase-transition-produced magnetic fields was stressed. Here, in order to study the correlation properties of an evolving cosmic magnetic field, however, a simple toy model was employed, which replaced the full MHD equations. Since then, there are been many efforts to study, more accurately, the evolution of turbulent magnetic fields [36][37][38][39][40][41][42][43][44][45][46][47][48][49][50][51][52].
However, only in the work [53], it was realized that, other than turbulent phases, phase-transition-generated magnetic fields can experience dissipation processes induced by free-streaming neutrinos, photons, or hydrogen atoms. Indeed, such cosmic components can have mean free paths which exceed the typical correlation scale of a primordial magnetic field, so they can be considered as free-streaming species with respect to the magnetic field itself. In this case, the dissipation of magnetic energy is ruled by a "drag" term in the Navier-Stokes equation, instead of the usual diffusion term proportional to the kinematic viscosity.
Such a possible phase in MHD, called "dragged magnetohydrodynamics" [54,55], has been fully analyzed in [54,55], where a systematic study of the evolution of cosmic magnetic fields has been performed by direct numerical integration of MHD equations. To our knowledge, numerical simulations of dragged MHD have so far been performed only in [54,55] where, moreover, a justification of the results has been given by using simple scaling arguments.
The aim of this paper is to analyze the evolution of phase-transition-generated magnetic fields, in turbulent and viscous free-streaming regimes, by using suitable mean-field-approximated magnetohydrodynamic equations.
In particular, we will derive the evolution laws for the magnetic energy and correlation length both in non-helical and helical cases. We anticipate that our results agree with the numerical results in [54,55].
The paper is organized as follows. In the next section, we introduce the relevant equations and the physical quantities of dragged MHD. Also, we find some new exact results in dragged MHD, and we introduce the transformations connecting the evolution of a primordial magnetic field in a (static) Minkowski spacetime to that in a flat Friedmann (expanding) universe. In section III, we apply some rigorous scaling arguments, first introduced by Olesen for the study of turbulent MHD, to dragged MHD. In section IV, we find the equations governing the evolution of the magnetic energy and helicity by using three different mean-field-theory approaches. In section V, we solve these equations for some specific initial conditions, and find the evolution laws of the magnetic intensity and correlation length in radiation dominated universe for the case of free-streaming neutrinos and photons. In section VI, the case of free-streaming photons and hydrogen atoms in matter dominated universe is analyzed. In section VII, we justify some scaling arguments of [54,55] in the light of our exact analytical results obtained in section II. In section VIII, we discuss the evolution of freely-decaying turbulent magnetic fields, comparing the various and often different results present in the literature. Also, we complete our work started in [49] by extending our previous results, obtained in Minkowski spacetime, to the case of a flat Friedmann universe. In section IX, we study the dependence of our results on the choice of the initial cut-off of the assumed, initial magnetic power spectrum. In section X, we discuss our results and compare them to the results of [54,55]. In section XI, we draw our conclusions and give possible suggestions for further investigations. Finally, appendixes A,B, and C contain some technical details of our analytical computations.

II. DRAGGED MAGNETOHYDRODYNAMICS
In this section, we introduce the equations and physical quantities of interest for dragged magnetohydrodynamics. Some new exact results are derived, some scaling arguments are discussed, and the transformations relating the evolution of a magnetic field in Minkowski and Friedmann spacetimes are introduced.

IIa. Preliminaries
Full magnetohydrodynamic equations.-The full magnetohydrodynamic equations in an expanding Friedmann universe for a relativistic imperfect fluid were first derived by Jedamzik, Katalinić, and Olinto in [53]. A full analysis of these equations would require the inclusion of effects deriving from inhomogeneities in the matter density, pressure field, particle number, viscosity, temperature, etc., those coming from compressibility of the plasma flow, as well as the study of relativistic effects in the case of ultra-relativistic fluid motion.
Inhomogeneities in the magnetohydrodynamic variables generally leads to magnetohydrodynamic modes (e.g., Alfvén waves, slow and fast magnetosonic waves [53]) which propagates through the magnetized plasma, and their decay can result in a dissipation of magnetic energy stored in the small-scale fluctuating part of a cosmic magnetic field.
Compressibility of plasma flow is generally expected to be realized in a cosmological context. For example, a recent analysis by Giovannini [56] has argued that, due to adiabatic inhomogeneities of the scalar curvature, the cosmic plasma flow at large scales could be compressible in the period of time between electron-positron annihilation around T e + e − ∼ 1MeV [57] and last scattering at T ∼ 0.3eV [57].
Also, ultra-relativistic fluid motion can develop during radiation era and could play, in principle, in important role in the evolution of cosmic magnetic fields during radiation era.
Nevertheless, in the following and due to the complexity of the full MHD equations, we limit ourselves to the study of a simplified version of magnetohydrodynamic equations. We will work in the hypothesis where the small-fluctuating parts of MHD variables do not propagate (thus neglecting the development and the damping of the previously indicated MHD modes), as well as in the limit of incompressibility of the cosmic plasma and small velocity fields (v/c ≪ 1, where v is the typical velocity of bulk fluid motion and c is the speed of light). The first approximation implies that the decay laws of a cosmic magnetic field we will find below must be considered as conservative, in the sense that a (slightly) faster magnetic energy decay could happen as a result of the inclusion in the analysis of propagating MHD modes. The validity of the second approximation, as discussed in [55], depends on the initial magnetic field intensity as well as on the epoch considered. A general expectation is that very strong magnetic fields render compressible the fluid motion after the decoupling of photons from the cosmic flow. In the analysis of [55], however, it is shown that incompressibility is justified for comoving magnetic fields intensities below the value 6 × 10 −11 G.
Having pointing out the possible limitations of our simplified analysis, it is worth noticing that we will use the same set of approximate MHD equations used in [54,55] to derive, numerically, the evolution properties of a cosmic magnetic field.
Reduced magnetohydrodynamic equations.-The "reduced", incompressible, 2 Newtonian magnetohydrodynamic equations in the presence of a drag term and in Minkowski spacetime are [54,55] Dv where D/dt = ∂ t + v· ∇ is the so-called hydrodynamical derivative, v is the velocity of bulk fluid motion and ∇ · v = 0, v A = B/ √ ρ + P is the Alfvén velocity, B is the magnetic field, and ρ and P are the energy density and the thermal pressure of the fluid, respectively. The frictional coefficient α and the resistivity η are dissipative parameters, and are determined by microscopic physics.
We will consider throughout the paper the case of constant ρ and P . In the case of a Friedmann universe, where ρ and P evolve in time, the quantities of interest, as we will see in section IIb, areρ andP , instead of ρ and P , with the "tilded" quantities being time independent. For this reason, and for the sake of simplicity, we can take ρ + P = 1 in Eq. (1) The Navier-Stokes equation (1) and the autoinduction equation (2) can then be written, respectively, as where we have introduced the Lorentz force F L = J × B, and the damping force F d = αv, with J = ∇ × B being the magnetic current. Let us now introduce the relevant physical quantities in (dragged) MHD. The magnetic and kinetic energy densities are defined as where E B = 2πk 2 |B(k)| 2 and E v = 2πk 2 |v(k)| 2 are the magnetic and kinetic energy density spectra. Here, B(k) and v(k) are the magnetic and kinetic fields in Fourier space, with k being the wavenumber and k = |k|. The magnetic and cross helicity densities are respectively, where H B = 4πk 2 A(k)·B * (k) and H c = 4πk 2 v(k)·B * (k) are the magnetic and cross helicity density spectra, A is the vector potential, and an asterisk denotes complex conjugation. It is worth noticing that, for all magnetic and kinetic field configurations, the magnetic and cross helicity spectra satisfy the following realizability conditions: The total energy E tot = E B + E v and the cross-helicity are conserved quantities when α = η = 0, while the magnetic helicity is conserved when η = 0. This follows directly from their evolution laws, which can be straightforwardly derived from the MHD equations (3) and (4). They are 2 The kinetic energy associated to a turbulent fluid is Ev = (1/2)ρ v 2 , where ρ is the energy density of the universe, and v is the typical velocity of bulk fluid motion associated to turbulence. Incompressibility of the primordial plasma means that the acoustic Mach number, Ms = v/vs, where vs is the speed of sound, is much less then unity [24]. For example, in the radiation era, when the cosmic fluid is ultra-relativistic and then vs = 1/ √ 3, the incompressibility condition reads v ≪ 1/ √ 3.
where ω = ∇ × v is the so-called vorticity. We also introduce the relevant length scales in (dragged) MHD, i.e., the so-called correlation lengths, which are the characteristic lengths associated with the large magnetic and kinetic energy eddies of turbulence. These are defined by With the aid of the magnetic and kinetic correlation lengths, it is straightforward to transform the above "local" realizability conditions on the magnetic and cross helicity spectra, into the following "global" realizability conditions: Finally, it is useful to define the kinetic Reynolds number, the magnetic Reynolds number, and the Prandtl number as Re = v/(lα), Re B = vl/η, and Pr = Re B /Re, respectively, where v and l are the typical velocity and length scale of the fluid motion.

IIb. The equilibrium state
A dragged phase develops when the kinetic Reynolds number is small [54,55]. For this reason, we work in the case where Re ≪ 1. Defining the "drag time" τ d = α −1 , we have |D t v|/|F d | ∼ τ d /τ eddy = Re, where we used ∂ t ∼ 1/τ eddy and introduced the so-called (kinetic) eddy turnover time τ eddy = l/v. Accordingly, for very low kinetic Reynolds numbers, we can neglect the left-hand side of Navier-Stokes equation (3). This means that the system will approach asymptotically a state where the Lorentz force equilibrates the damping force [54,55], In this "equilibrium state", the autoinduction equation (4) reduces to an equation which depends only on the magnetic field, and will be basis for the study of the evolution of the magnetic field itself. Multiplying both sides of Eq. (14) by B, we find that the cross helicity is zero in the equilibrium state, Also, there is no equipartition between kinetic and magnetic energy in this state. Indeed, defining the Γ-ratio we get where B is the typical intensity of the magnetic field. Moreover, τ d and τ eddy are related by , where a dot indicates a time derivative. Consequently, in the limit of low kinetic Reynolds number and large magnetic Reynolds number, Eq. (9) can be approximated with The above equation connects the evolutions of magnetic and kinetic energies, and will be useful in the following. Finally, since we are considering a very large Prandtl number, we can neglect the second term in the right-hand side of Eq. (10). In fact, we have |η d 3 x J · ω|/|αH c | ∼ 1/Pr. Accordingly, we getḢ c ≃ −αH c , whose solution gives the evolution law of the cross helicity in the equilibrium state. In the next subsection, we show that the above equation is indeed (asymptotically) compatible with Eq. (16).

IIc. Scaling arguments
Let us suppose, and this is indeed the case in a cosmological context [54,55], that the drag coefficient α evolves in time as a simple power law, where t i is a reference time.
In numerical simulation of dragged magnetohydrodynamics [54,55], it is observed that the magnetic energy scales in time as a power law if the initial magnetic spectrum is assumed to be a simple power of the wavenumber. Let us then parameterize the magnetic energy as Inserting Eqs. (22) and (23) into Eq. (20), we find the evolution laws for the kinetic energy and Γ-ratio, respectively, where Γ(t i ) = (β/2)[τ d (t i )/t i ]. From this last relation we get β > 0, which means that the magnetic energy decreases in time during a dragged MHD phase (as in the case of turbulent MHD). Moreover, recalling that Γ ∼ Re = v/αl and taking v ∼ E 1/2 v , we find that the typical length scale evolves in time as giving, in turn, the evolution law for the eddy turnover time. Finally, inserting Eq. (22) in Eq. (21), we get if a = −1, and if a = −1. In both cases, since t i ∼ τ eddy (t i ) ≫ τ d (t i ) [see Eqs. (18) and (19)], we get that the cross helicity goes to zero asymptotically (t ≫ t i ), in agreement with Eq. (16). It is amazing to observe that, even without knowing the value of the exponent β, simple scaling arguments, when combined with the exact laws (20) and (21), give the evolution laws for the Γ-ratio, the eddy turnover time, and the cross helicity.

IId. Expanding universe
In the case of an expanding universe (with zero spatial curvature) described by a Friedmann-Robertson-Walker line element, it has been shown that the dragged MHD equations are the same as the Eqs. (1) and (2) provided that time, coordinates and dynamical variables are replaced by suitable quantities. The form of these new "tilded" quantities is different in radiation an matter eras.
Radiation era: comoving variables.-In a radiation-dominated universe the tilded quantities are [54,55]: where R(t) is the expansion parameter,t is the so-called conformal time, and we note that v is not scaled.
Matter era: supercomoving variables.-In a matter-dominated universe the tilded quantities are [54,55]: Due to the formal coincidence of the dragged MHD equations in Minkowski and Friedmann spacetimes, we can analyze the evolution properties of a primordial magnetic field in both cases in a similar way. For definiteness, we consider first the case of a non-expanding universe and then translate the results to the case of interest of an expanding universe. We observe that, since ρ ∝ P ∝ R −4 in the radiation-dominated era and ρ ∝ P ∝ R −3 in the matter-dominated era, the quantityρ +P is a constant during the evolution of the universe. This result was used in section IIa.

III. OLESEN'S APPROACH
In this section, using a scaling approach first introduced by Olesen in [36] for the case of freely-decaying MHD turbulence, we express the unknown exponent β introduced in section IIc as a function of the exponents of the assumed initial power-law magnetic spectrum, E B (k, t = 0) ∝ k p , and of the power law for the drag coefficient, Non-helical case.-We are interested in the case of large magnetic Reynolds numbers, and then we neglect the dissipative term in the autoinduction equation (4). We can then write the dragged MHD equations as A direct inspection shows that, under the scaling transformations x → ℓ x and t → ℓ 1−u t, equations (38) and (39) admit solutions of the type where ℓ > 0 is the "scaling factor", and u and m real parameters. Differentiating Eq. (42) with respect to ℓ, and putting ℓ = 1 afterwards, we get α(t) ∝ t m/ (1−u) , and then m = a(1 − u), since we are assuming that α ∝ t a . Defining p = −(2 + u + m) and using Eqs. (40) and (41), we straightforwardly obtain where we have defined r = [3(1 + a) + 2p ]/(1 − a) for notational convenience. Here, λ B and λ v are constants, while ψ B and ψ v are arbitrary scaling-invariant functions of their arguments. What is observed in dragged MHD [54,55] (as well as in turbulent MHD [54,55]) is that the evolution of the magnetic spectrum in the non-helical case proceeds through the so-called "selective decay" (see section V), so that the initial magnetic spectrum retains its form (for lengths below the characteristic dissipation scale) for all times. This means that if we assume that the initial magnetic spectrum is a simple power of the wavenumber k, then the coefficient p is the index of that power. In fact, we have E B (k, 0) = λ B k p ψ B (0), a result which is true also in turbulent MHD (see [36]). Inserting Eq. (43) into Eq. (5) and Eq. (44) into Eq. (6), we obtain Eqs. (23) and (24), respectively, where now β = (1 − a)(1 + p)/(3 + p) is expressed as a function of p and a, Equations (25) and (27) are accordingly confirmed. Also, we find that confirming the scaling relation (26). Helical case.-In the case when the magnetic helicity is different from zero, we can still apply the Olesen's approach to find the scaling behavior of the magnetic helicity spectrum. We find, straightforwardly, where µ B is a constant and φ B an arbitrary scaling-invariant function. In dragged MHD [54,55] (as well as in turbulent MHD [24,54,55]), the magnetic spectrum evolves through a phase of so-called "inverse cascade" (see section V), during which magnetic modes with small wavelengths get modified by the longer ones and viceversa. During this process, the information enclosed in the initial magnetic spectrum is lost, so that the coefficient p in Eq. (43) does not represent, necessarily, the index of the initial power-law spectrum. Indeed, since magnetic helicity turns out to be (quasi-)conserved in dragged MHD [54,55] (as well as in turbulent MHD [24,54,55]), we must have p = 0. In fact, equations (7) and (48) imply . Accordingly, the evolution laws for the magnetic and kinetic energies and correlation lengths are obtained in the helical case taking p = 0 in the evolution laws previously found for the non-helical case: Equations (45), (46), (47), and equations (49), (50), (51), are in agreement with the results of [54,55].

IV. MEAN-FIELD APPROXIMATION
In this section, working in mean-field-theory approximation, we derive the evolution integro-differential equations for the magnetic energy and helicity spectra, which will be solved later on in section V. We start from the "simplest" mean-field approximation first introduced in [58], we then study the dragged MHD equations in the so-called "onepoint-closure" [24] and "quasi-normal" [49] approximations.

IVa. Cornwall's approximation
It is useful to re-write the autoinduction equation (15) as We now proceed as in [58] by replacing the quadratic terms J · B and B 2 in the above equation with J · B → 1 3 J · B and B 2 → 1 3 B 2 , where the brackets ... indicate a suitable average to be defined later. This operation allow us to linearize the autoinduction equation, Whatever is the averaging operation, the averaged quantities in the above equation must satisfy a consistency relation coming from conservation of magnetic helicity in the case of null dissipation. Indeed, observing that Eq. (53) gives we get Assuming now, as in [58], that the operation of averaging is just a volume average, (55) is satisfied for all magnetic field configurations and, taking into account Eqs. (11) and (5), that J · B = −Ḣ B /(2η) and B 2 = 2E B . Finally, we can re-write Eq. (53) as where we have introduced the quantities α (56) describes the well-known α−dynamo effect [24], in which the dynamo, the turbulent diffusion, and the total effective diffusion coefficients, are given by α B , η T , and η eff , respectively.

IVb. One-point-closure approximation
Cornwall's arguments find a full justification in the framework of one-point-closure approximation, originally introduced in turbulent MHD [24], and that we apply now to dragged MHD. Equation.
-Let us assume that v and B can be decomposed into an average part varying only on large scales and a weak, small-scale fluctuating part, with v = B = 0, and | v| ≪ |v 0 |, | B| ≪ |B 0 |. In a moving coordinate system such that v 0 = 0, the equation governing the evolution of the mean magnetic field is a dynamo equation with dynamo and turbulent diffusion coefficients given by [24] α ′ In Eq. (59), we have introduced the the so-called kinetic helicity density [24], while in the last equality of Eq. (60), we used Eq. (17). Starting from Eq. (58) and imposing the conservation of the mean magnetic helicity T and in turn, using Eqs. (59) and (60), we find H v = (Γ/2η)Ḣ B . Therefore, we can re-write the dynamo coefficient as α ′ B (t) ≃ −Ḣ B Γτ eddy /(6η). Both the dynamo and turbulent diffusion coefficients are then consistent with the expressions for α B and η T because of Eq. (19). Solution.
-In order to solve the dynamo equation (56) [or, which is the same, Eq. (58)] it is useful to introduce the orthonormal helicity base {e + , e − , e 3 }, with e ± = (e 1 ± ie 2 )/ √ 2 and e 3 = k/k, where {e 1 , e 2 , e 3 } is a righthanded orthonormal base. In the helicity base, the magnetic field in Fourier space can be decomposed as B(k, t) = B + (k, t) e + + B − (k, t) e − , where B + (k, t) and B − (k, t) represent the positive and negative helicity components of B(k, t), respectively. Equation (56) becomesḂ ± = ±α B kB ± − η eff k 2 B ± in Fourier space. The solution of the above equation is easily found, where we have defined the "dynamo" and "dissipation" lengths, ℓ α (t) = t 0 dt α B and ℓ 2 diss (t) = t 0 dt η eff , respectively. In the helicity base, the energy and helicity spectra read To proceed further, let us suppose that . Indeed, we have made the simplifying assumption that the initial helicity is just a fraction of the initial maximal helicity. In other words, we restrict our analysis only to so-called magnetic fields with (initial) "fractional helicity". Inserting Eq. (61) in Eqs. (62) and (63) we get, respectively, The above equations are integro-differential equations for the magnetic energy and helicity, and will be solved in section V.

IVc. Quasi-normal approximation
We now show that the results in the above two subsections, can be also derived using a two-point-closure approximation, also known as "quasi-normal" approximation.
We start by observing that the autoinduction equation (52), in Fourier space, reads: where ε ijk is the totally antisymmetric tensor and summation over repeated indexes is understood. In quasi-normal approximation, the four-point magnetic correlator is decomposed, in terms of two-point correlator, as [24]: where ... denotes an ensemble average. Multiplying Eq. (66) respectively by B * i (k) and A * i (k), and then averaging out, we get As it is easy to check, the solution of the above equations in the case of magnetic fields with initial fractional helicity is given by Eqs. (64)- (65), so that the one-point-and two-point-closure approximations give exactly, and perhaps surprisingly, the same identical results.

V. RESULTS
In this section, we find the evolution laws of the magnetic intensity and correlation length in a radiation-dominated universe for the case of free-steaming neutrinos and photons, by solving Eqs. (64)-(65) for some specific initial conditions. When working in a Minkowski spacetime, we take t = 0 as the initial time for the sake of simplicity.

Va. Initial conditions
Magnetic spectrum.-For the sake of simplicity, we assume that the initial magnetic energy spectrum is represented by the simple function where λ B and ℓ B are constants. For k ≪ ℓ −1 B , the magnetic energy spectrum possesses a power-law behavior with exponent p which, to avoid infinities in the magnetic energy, we assume to be greater than 1. For the same reason (finiteness of energy), an exponential cut-off has been introduced. The quantity ℓ B can be related to the initial correlation length, is the Euler gamma function [59]. It is worth noting that the power-law behavior for the initial magnetic spectrum is predicted by many models of generation of cosmic magnetic fields in the early universe. In particular, inflation-produced magnetic fields [10] usually exhibit a power-law spectrum when they re-enter the horizon. Also, most of the generating mechanisms after inflation [26] repose on microphysical processes acting only on small scales. Therefore, the structure of the resulting magnetic field appears as a set of uncorrelated (Gaussian-distributed) eddies, which implies a k 2 -like spectrum. To see this, suppose that the magnetic field at the initial time t = 0 is a stochastic variable with gaussian distribution, is then proportional to k 2 . In fact, observing that where in the last equality we used the properties of Gaussian integral [60], we get E B (k, 0) = λ B k 2 .
Recently enough [61], however, it has been shown that the analyticity of the magnetic field correlator B i (x)B j (y) defined on a compact support, together with the divergenceless of B, implies that the spectral index p has to be even and equal or larger than 4. Moreover, in [51], it has been shown that, starting from suitable initial conditions which should reflect the production mechanism of a magnetic field from bubble collision in a first-order phase transition, a Batchelor spectrum E B ∝ k 4 is established after a short time interval and at small wavenumbers. Nevertheless, we can always assume that on very large scales (very small wavenumbers k) the initial spectrum satisfies the above requirements, but that on smaller scales, where the gross of the magnetic energy is stored, it behaves like k 2 . In any case, in the following we leave the index p as a free parameter, privileging the cases p = 2, 3, 4 when we show graphically our results (for the choice p = 3, see below).
Drag coefficient.-We assume that the drag coefficient α scales in time following the simple law where τ is the normalized time and γ(0) and a are constants. The above parametrization is useful since α(t) ≃ α(0) for t ≪ τ eddy (0), when the system is, as we will see below, in a quiescent phase [i.e. the integral quantities like E B (t) and ξ B (t) remain almost constant in time], and α(t) follows asymptotically a power-law behavior α(0)[τ /γ(0)] a for large times, when magnetohydrodynamic effects operate in changing the state of the system. Using the results in [54,55], we find that in comoving variables a = −3 for the case of free-steaming photons, while a = −4 for the case of free-steaming neutrinos. The constant γ(0) is generally different from unity, and its explicit expression will be derived in section Vh.
Setup.-Numerical integration of dragged MHD equations has been performed only in [54,55]. There, the simple case of constant dragged coefficient was analyzed. Also, one of the cases discussed in [54,55], was that with p = 3. Hence, in the light of the above discussions, we consider the three cases for the non-helical case, and the three cases Gaussian, photon case, (3, 0, 1), Banerjee-Jedamzik case, (4, −4, 10 −15 ), causal, neutrino case, for the helical case. Note that in the helical, Banerjee-Jedamzik case, the magnetic field is taken to be maximally helical, as it is in the numerical analysis of [54,55]. Also, the very small value h = 10 −15 is taken to clearly show the transition from the selective decay phase to the inverse cascade phase to be discussed in section Vd.

Vb. Master equations
Inserting Eqs. (64) and (65) in Eqs. (5) and (7) we find, respectively, with is the Kummer confluent hypergeometric function [59], and we have introduced the quantities It is useful, for the following discussion, to define accurately the magnetic Reynolds number and the eddy turnover time as where we used the magnetic correlation length and the root-mean-square value of the velocity field, v 2 rms = d 3 x v 2 (x, t) = 2E v , as typical length scale and velocity, respectively. With the aid of the above definitions, the integro-differential equations (77) and (78) for E B and H B can be transformed into ordinary differential equations for the quantities ζ diss and ζ α . We have with E B and H B as a function of ζ diss and ζ α given by Eqs. (77) and (78), respectively. The constant is of order unity [see Eq. (19)]. In section Vf, we will found that this order-unity constant is indeed a parameter which weakly depends on p.
Equations (77)-(78) and (83)-(84) are our master equations in the study of the evolution of the magnetic field.

Vc. Solutions: non-helical case
In the non-helical case, h B = 0, the solution of Eqs. (83)-(84) is such that ζ α (t) = 0, that is H B (t) = 0 for all times. Then, from Eqs. (77), (81), and (64) we get For large magnetic Reynolds numbers, the first term in the right-hand-side of Eq. (83) can be neglected with respect to the second one. In this case, the solution of Eq. (83) is where c 0 is given in Appendix A. For τ ≫ 1, we have where c 1 is given in Appendix A, so that where c 2 , c 3 , c 4 are given in Appendix A. Equations (91), (92), and (93), are in agreement with Eqs. (45), (47), and (43), respectively. We deduce that the mean-field-theory approach and the simple scaling arguments applied to MHD equations are compatible. The mean-field-theory approach gives us two more pieces of information with respect to scaling arguments. Firstly, the time when the system enters into the scaling regime is, looking at Eqs. (86), (87), (88), and (89), approximatively given by the initial eddy turnover time times γ(0). Secondly, comparing Eq. (93) with Eq. (43) we get the expression for the scaling function ψ B (x): In Fig. 1, we plot the spectrum of the magnetic energy at different times for the three cases in Eq. (75) [in all figures of this paper, we take γ(0) = 1]. It is clear that the decay of the magnetic energy and the growth of the magnetic correlation length (which is, roughly speaking, the inverse of the wavenumber corresponding to the peak in the magnetic spectrum) advance through selective decay, a phenomenon well described in [39] in the case of turbulent MHD (see also [24]). There is no direct transfer of magnetic energy from small scales to large scales but modes with larger wavenumbers just decay faster than those whose wavenumbers are small. Accordingly, as the time goes on, the magnetic field survives only on larger and larger scales, as it is evident in Fig. 1. The dotted lines in Fig 1 represent the magnetic energy spectra at a certain time t/τ eddy (0) = τ E (to be defined in section Vd) when, if the field were helical, the magnetic field would enter in a phase of inverse cascade (see below). In the case at hand, the magnetic field has no helicity and, indeed, nothing special happens when τ reaches τ E .

Vd. Solutions: helical case
In the helical case, the dynamo and dissipation lengths, as well as the ratio ζ α /ζ diss , are increasing functions of time. This is inferred by the results of numerical integration of Eqs. (83)-(84) as shown, for example, in the left panel of Fig. 2. [Here, and in the following, we consider the case of a very high magnetic Reynolds number, Re B = 10 15 ]. In the Appendix B, we show that there are two different regimes depending on the value of ζ α /ζ diss . For ζ α /ζ diss ≪ 1 the system behaves as if the magnetic helicity were zero, so that the system evolves by selective decay. The asymptotic solutions are then the same as obtained in section Vc. For ζ α /ζ diss ≫ 1, instead, the mechanism operating in the evolution of the system is such that energy is transferred from small to large scales. In turbulent MHD, this mechanism is known as inverse cascade and will be discussed, in more detail, later on in this subsection.
The following asymptotic (τ → ∞) solutions of Eqs. (83) and (84) are derived in the Appendix B in the case of very large magnetic Reynolds numbers: and where c 5 , c 6 , c 7 , c 8 are given in Appendix A. The above relations have been obtained supposing that the magnetic helicity is almost constant during the evolution of the system. This assumption is well justified by the results of the numerical integration of Eqs. (83) and (84), as shown, for example, in the right panel of Fig. 2.
In the upper panels of Fig. 3, we show the magnetic energy and correlation length as a function of time for the three cases in Eq. (76). Dotted lines are the asymptotic expansions (97) and (98) for τ ≥ τ E and τ ≥ τ ξ , respectively, and the analytical results (86) and (87) [with ζ diss given by Eq. (89)] for 0 ≤ τ ≤ τ E and 0 ≤ τ ≤ τ ξ , respectively. Here, τ E and τ ξ are the times when, approximatively, the integral quantities E B and ξ B enter into the inverse cascade regime. We can find them by simply matching the asymptotic solutions in the two different regimes, namely equating Eqs. (91) and (92) to Eqs. (97) and (98), respectively. We find, a)]. In the case h B ≪ 1, the above results can be approximated as where c 9 , c 10 are given in Appendix A.
We also note that, in the inverse cascade regime, the product of magnetic energy and correlation length is a p, a, h B t Τ eddy 0 quasi-conserved quantity. In fact, as we show in the Appendix B, we have From the above equation, we deduce that a magnetic field with fractional helicity becomes maximally helical approximatively after the system enters into the inverse cascade phase. The time when this happens, τ H , can be found by matching the product of asymptotic solutions Eqs. (91)- (92) and (101). We get where c 11 is given in Appendix A. In the lower panel of Fig. 3, we show the quantity 2E B ξ B /H B as a function of time for the tree cases in Eq. (76). The asymptotic expansions for τ ≤ τ H and τ ≥ τ H are practically indistinguishable from the numerical solutions.
In the helical case, the evolution laws of magnetic energy and correlation length do not depend on the index of the initial magnetic energy spectrum. Moreover, apart logarithmic corrections, and for the case of constant drag coefficient (a = 0), we find the t −1/3 and t 1/3 laws which are indeed observed in numerical simulations of dragged MHD [54,55].
It is evident from all the above figures that the (approximate) analytical expansions fit very well the numerical solutions. Moreover, the magnetic helicity is, as expected, an almost constant function of time. Because of the conservation of the magnetic helicity, small-scale modes are not dissipated during the decay but their energy is transferred to larger scales. This process of inverse cascade is manifest in Fig. 4 (see, in particular, the magnetic spectra). It is worth noting that, as explained above, this mechanism begins to operates only after a certain time, that is when the selective decay mechanism ends. We can suitably define the time when the inverse cascade begins, τ E , as the time when the maximum of the magnetic spectrum meets the initial spectrum (see the magnetic spectra in Fig. 4). Assuming for simplicity h B ≪ 1, we get (see Appendix C) where c 12 is given in Appendix A. Looking at Eqs. (74), (99), (100), (102), (103), and roughly speaking, we can say that the time when system enters into the inverse cascade regime is t inv cas ∼ h τ eddy (0). Finally, we note that magnetohydrodynamic effects cannot change the characteristics of the magnetic field on scales well above the integral scale (the magnetic correlation length); in particular, in the limit k → 0, the initial spectrum must retain its form for all times. This is due to the fact that the transfer of magnetic energy from small to larger scales (inverse cascade) is not instantaneous and needs some time to develop. As a result, there will be always, at a given time, a tail in the small wavenumber part of the magnetic spectrum which is leaved unprocessed by the inverse cascade. This phenomenon of inefficiency of inverse cascade for k → 0 is visible in the case shown in the upper left panel of Fig. 4, where it is evident that the magnetic energy spectrum retains its initial power law form at very large scales. Such a phenomenon is however present in all the cases we have analyzed, although it is not evident in the other panels of Fig. 4 since the magnetic energy and helicity spectra are not shown for very small wavenumbers.

Ve. Solutions: maximally helical case
In the case where the initial magnetic field is maximally helical, h = 1, we get two exact analytical results. Firstly, from Eqs. (64) and (65) we find that the magnetic field remains maximally helical for all times, H B (k, t) = 2k −1 E B (k, t). Secondly, integrating the above expression in k we have E B (t) ξ B (t) = H B (t)/2, a result usually used in the literature for the evolution of a (maximally) helical magnetic field, but never fully justified.
where we used Eq. (91) for the non-helical case, and Eq. (97) for the helical case. Equation (104) Finally, starting from the definition of δ(0) in Eq. (85), and using Eqs. (17) and (82) evaluated at τ = 0, we find which relates the kinetic energy at the onset of the dragged phase to the initial values of the drag coefficient, the magnetic field energy, and the correlation length.

Vg. Exiting the dragged phase
Neutrino and photon mean free paths.-The dragged phase is characterized by the condition that the comoving mean free path of neutrinos, ℓ (ν) mfp , and/or photons, ℓ (γ) mfp , is much greater than the comoving magnetic correlation length. Therefore, it is not obvious that this condition in maintained during the evolution of the system, since ξ B and both ℓ mfp ∝ T −2 ∝ R 2 ∝t 2 in radiation-dominated era [54,55], we have ℓ (ν) mfp ∝ τ 4 and ℓ (γ) mfp ∝ τ 2 . Here, we used the results in section IId, the fact that R ∝ t 1/2 in radiation era (with t being the cosmic time), and the fact that the temperature T and the scale factor R are related by R ∝ g −1/3 * ,S T −1 [57], with g * ,S being the effective number of entropy degrees of freedom at the temperature T . On the other hand, the maximum growth for the correlation length happens for the helical case and for the neutrino free-streaming case (a = −4) [see Eqs. (92) and (98)]. Namely, ξ B increases at most as ξ B ∝ (ln τ ) −1/3 τ 5/3 . As a consequence, if the system is in a dragged phase such that ξ B ≪ ℓ (ν) mfp , ℓ (γ) mfp , it will remain in it for all times.
However, the evolution laws for the magnetic energy and correlation length has been obtained in the equilibrium state, where the kinetic Reynolds number is much smaller than unity, Re ≪ 1. This quantity generally evolves in time as the typical velocity and correlation length of the fluid motion, and the drag parameter change in time. Consequently, there could exist a time when the condition Re ≪ 1 ceases to be satisfied and the system exits the dragged phase.
Kinetic Reynolds number.-Defining the kinetic Reynolds number accurately as we find Re = (−Ė B ) 1/2 ξ −1 v α −3/2 , where we used Eq. (20). The asymptotic expansions (for τ ≫ 1) of the kinetic Reynolds number in the non-helical and helical cases, and for large magnetic Reynolds numbers, are where we used the fact that ξ v (t) ∝ ξ B (t) [see the first relation in Eq. (47)

Vh. Evolution laws in expanding universe
Returning to the "tilde" notation and indicating the time when the dragged phase begins ast i , equation (22) reads Introducing the new variableτ asτ =t τ eddy (t i ) , we can write Eq. (109) as α(t) = α(t i )[τ /γ(t i )] a with γ(t i ) =t i /τ eddy (t i ). In radiation-dominated era, it is easy to see that γ(t i ) can be expressed as γ(t i ) = [d H (t i )/ξ B,phys (t i )] v rms (t i ), where d H = 2t is the length of the Hubble horizon at the time t, with t being now the cosmic time. Therefore, in the notation of section Va, we have γ(0) = N i v i , where we have introduced the initial number of magnetic domains per horizon length and we have defined v i = v rms (0) for the sake of simplicity. Moreover, the normalized timeτ in the radiation-dominated era is easily found to beτ Finally, using Eq. (106), we can relate v i to the initial values of the magnetic field strength, B 2 = B 2 rms = d 3 x B 2 (x, t) = 2E B , the physical magnetic correlation length, ξ B,phys = Rξ B , and drag coefficient, as where κ 0 is given in Appendix A.
[We note that in the notation in which the quantity ρ + P is not taken to be 1, the magnetic field strength B(R i ) in Eq. (113) should be replaced by Non-helical case.-Translating the results previously found in the case of a Minkowski spacetime to the case of a flat Friedmann universe, we get for the comoving magnetic field strength, and for the comoving magnetic correlation length. Here, and κ 1 , κ 2 , κ 3 are given in Appendix A. In Eq. (116), we have assumed that κ 3 (N i v i ) −1/(1−a) > 1. If this is not the case, R 1 is (approximatively) equal to R i . Helical case.-In the case of helical magnetic fields, we find for the comoving magnetic field strength, and for the comoving magnetic correlation length. Here, , and κ 4 , κ 5 , κ 6 , κ 7 are given in Appendix A. In the above equations, we have assumed that ω is much greater than unity.

VI. DRAGGED MHD IN MATTER-DOMINATED UNIVERSE: FREE-STREAMING PHOTONS AND HYDROGEN ATOMS
We now discuss the evolution of a magnetic field in matter-dominated universe, when the drag term in the MHD equations is ruled by the free-streaming properties of photons or hydrogen atoms.

VIa. Evolution laws in supercomoving variables
The drag coefficients for photons and hydrogen atoms evolve as a function of the temperature as [54,55] α γ (T ) ∝ T 4 and α H (T ) ∝ T 3 , respectively. In supercomoving variables (namely in Minkowski spacetime), their evolution with time is then straightforwardly given by where, taking into account the results of section IId, we used the fact that the normalized timeτ in Eq. (110) is in matter-dominated era, when the scale factor evolves as R ∝ t 2/3 . In Eq. (120), and γ(0) = N i v i /3. The initial number of magnetic domains per horizon length, N i , is the same as in Eq. (111), the length of the Hubble horizon in matter-dominated era being d H = 3t (with t the cosmic time), and v i = v rms (0) is given, as in the case of radiation-dominated universe, by Eq. (113). In supercomoving variables, the master equation (83) reads, for the case at hand, Defining the new variable γ(0)[e −aτ /γ(0) − 1] = N i v i (R/R i ) 3/2 /3, the above equation becomes Comparing Eq. (124) with Eq. (83), we see that the evolution of the magnetic energy and correlation length can be derived by applying the substitutions τ →τ and a → a+ 1 in the solutions found in the case of a Minkowski spacetime in sections Vc, Vd, and Ve.

VIb. Evolution laws in expanding universe
Non-helical case.-Using the results in the above subsection, we straightforwardly find the evolution laws for the the comoving magnetic field strength, and for the comoving magnetic correlation length, where and κ 8 , κ 9 , κ 10 are given in Appendix A. In Eq. (127), we have assumed that κ 10 (N i v i ) 2/(3a) > 1. If this is not the case, then R 1 ≃ R i . Helical case.-In the helical case, we get for the comoving magnetic field strength, and for the comoving magnetic correlation length, where , and κ 11 , κ 12 , κ 13 , κ 14 are given in Appendix A. In the above equations, we have assumed thatω ≫ 1.

VIc. Exiting the dragged phase
Photon mean free path.-In the matter-dominated era, the comoving mean free path of photons scales in time as ℓ (γ) mfp ∝ T −2 ∝ R 2 . On the other hand, the maximum growth for the correlation length happens for helical fields in the case of free-streaming photons. Accordingly, ξ B increases at most as ξ B ∝ (ln R 3/2 ) −1/3 R 2 . As a consequence, if the system is in a dragged phase such that ξ B ≪ ℓ (γ) mfp , it will remain in it for all times. Hydrogen atoms mean free path.-The comoving mean free path of hydrogen atoms evolves as ℓ The correlation length in the hydrogen atoms free-streaming case, instead, increases as ξ B ∝ R 15/[2(3+p)] in the non-helical case, and ξ B ∝ (ln R 3/2 ) −1/3 R 5/2 in the helical case. As a consequence, the system will eventually exit the dragged phase in the helical case, and in the non-helical case for the case p ≤ 4/3 (we remember that we are assuming p > 1). Let us call t 1,exit the time when this happens, which can be defined by the condition ξ B (t 1,exit ) = ℓ From the above equation, we see that the kinetic Reynolds number is an increasing function of time in the physical cases of interest, a = −3, −5. Therefore, there will be a time t 2,exit when the system will leave the dragged phase characterized by Re ≪ 1. This time can be operatively defined by the condition Re(t 2,exit ) = 1. Accordingly, and in the cases of physical interest p = 2, 4, we have that the dragged phase terminates at the time t exit given by hydrogen, non-helical case, min{t 1,exit , t 2,exit }, hydrogen, helical case. (132) Finally, we observe that, as in the radiation-dominated era, we find that Re ∼ Γ, in agreement with Eq. (25).

VII. JUSTIFYING THE BANERJEE AND JEDAMZIK SCALING ARGUMENTS
In this section, we reconsider some scaling arguments used by Banerjee and Jedamzik in [54,55] to explain their results obtained by direct numerical integration of dragged MHD equations. Our goal here, it to show that some working hypotheses used in [54,55] find indeed a justification in the light of our analytical results obtained in sections II, V and VI.
Master equation.-In order to discuss the Banerjee and Jedamzik's scaling arguments, let us introduce the notations used in [54,55]. We define the magnetic energy in the mode k, . This is related to the magnetic energy spectrum by E B (k, t) = kE B (k, t). Accordingly, we can define the magnetic energy on the scale l ∼ 1/k as E B (l, t) = (1/l) E B (1/l, t). At the initial time t = 0, the magnetic energy on the scale l is proportional to a power of l, E B (l, 0) ∝ l −n , where n = 1 + p in the notation of [54,55].
In [54,55], it is assumed (without a full justification) that the evolution of the magnetic energy on the integral scale l B (t) (namely the scale under which the magnetic energy is dissipated) is ruled by the simple scaling equation where The above equation can be justified in the light of our results in section IIb. In fact, from Eq. (20), it follows that where E B (t) is the magnetic energy. In order to obtain the above equation, we eliminated E v in Eq. (20) in favor of Γ and then used the fact that αΓ ∼ 1/τ eddy [see Eq. (19)]. We now use the fact that the magnetic spectrum, in both the helical and non-helical cases, is peaked at a wavenumber k max (t) ∼ 1/ξ B (t) ∼ 1/l B (t), as it is not hard to see using the expressions for the magnetic energy spectrum and correlation length in section V, so that . This last observation, when combined with Eq. (134), gives Eq. (133). Now, we can transform Eq. (133) into an equation which does not involve the bulk field v (which enters in the equation through the quantity τ eddy ). To this end, we use Eq. (38) and approximate the nabla operator by ∇ ∼ 1/l B . We find that the typical velocity v and the typical magnetic field strength B are related through v ∼ B 2 /(αl B ) ∼ E B /(αl B ). In turn, this gives for the eddy turnover time τ eddy ∼ l 2 B /(αB). Inserting the above equation in Eq. (133) we get This is the master equation used in [54,55] to derive the evolution laws for the magnetic energy and correlation length. As we have just shown, this equation is a direct consequence of Eqs. (14) and (20). Non-helical case.-In [54,55], it is assumed that the magnetic energy on the integral scale l B is related to the integral scale itself by a simple power law This is indeed legitimated by the fact that, as we have shown in sections V and VI, the evolution of the magnetic spectrum in the non-helical case proceeds through selective decay, so that the initial magnetic spectrum retains its form (for length greater than the integral scale) for all times.
The above evolution laws are in agreement with those derived in mean-field approximation (see sections V and VI). Helical case.-In the helical case, the magnetic helicity is related to the magnetic energy on the integral scale l B through E B (l B ) l B ≃ H B (l B ). Since the magnetic helicity on the integral scale l B is a quasi-conserved quantity, and ξ B (t) ∝ t (1−a)/3 for free-streaming neutrinos and photons in radiation-dominated era, and E B (t) ∝ e (a/3) t and ξ B (t) ∝ e −(a/3) t for free-streaming photons and hydrogen atoms in matter-dominated era. Also in the helical case, the above evolution laws are in agreement with those derived in mean-field approximation (see sections V and VI).
The invariance of turbulent MHD equations under scaling transformations has been firstly used by Olesen [36] to study the evolution properties of the magnetic spectrum and correlation length. Although this is a very simple and elegant way to study MHD turbulence, this approach lives some unresolved questions, such as when the system enters into the scaling regime, what is exactly the shape of magnetic spectrum, etc.
The most famous toy model for hydrodynamic turbulence is the so-called shell model based on the original idea by Desnyanski and Novikov [62], developed later by Gledzer [63], Ohkitani and Yamada [64], and now know as GOYmodel. Subsequently, it was generalized by Gloaguen et al. to include the case of magnetic turbulence [65]. The basic idea on which repose the GOY-model is that the interaction due to non-linear terms in the MHD equations are local in k-space. Then, convolutions in MHD equations are approximated by sums over the nearest and next nearest neighbors in a discretized wavenumber space (the shell space). Moreover, the shell-model equations mimic the full MHD equations retaining their main characteristic, such as energy conservation in the case of null dissipation, phase space conservation, etc. Also, there exists a continuous version of discretized shell model in which the distance between shells goes to zero. It was firstly introduced by Parisi [66] in a hydrodynamical context and later extended to the case of MHD turbulence by Brandenburg,Enqvist,and Olesen [37]. Although many properties of turbulence, such as energy transfer, general spectral properties, etc., have been studied successfully using shell models, the physical validity of such toy models remains an open question (for a review on shell models of magnetohydrodynamic turbulence, see [67]).
The EDQNM approximation is the most widely two-point-closure approximation used for studying hydrodynamical turbulence. It was first proposed by Millionshtchikov [68], developed by Proudman and Reid [69], and extended to MHD turbulence by Pouquet et al. [70]. It consists in replacing correlation functions by suitable statistical averages [45], neglecting (conveniently) high-order correlators, and using a phenomenological expressions for the so-called "eddy damping rate". Even if the resulting equations are more simple than the full MHD equations, they can be solved only numerically and, moreover, some unresolved questions concerning EDQNM approximation still exist [24].
In [49], analytical results for the evolution of the magnetic energy density and correlation length were obtained in mean-field-theory approximation. Below, we briefly review such an approach to the study of turbulent MHD equations, and extend the results obtained in non-expanding universe to the case of a flat Friedmann spacetime.
A direct integration of the full set of turbulent MHD equations is certainly the best way to study the dynamics of freely-decaying MHD turbulence. However, due to the complexity of such highly non-linear equations, the numerical results in the literature are often contrasting I: Evolution laws for the magnetic energy, EB(t) ∝ t −σ , and magnetic correlation length, ξB(t) ∝ t ς , in helical and non-helical freely-decaying magnetohydrodynamic turbulence. All results were obtained assuming an initial power-law for the magnetic energy spectrum, EB(k, 0) ∝ k p , and constant dissipation parameters. The asterisk refers to analytical solutions of continuous shell model. In Tab. 1, indeed, we show some results found in the literature concerning the evolution laws of magnetic energy and correlation length in freely-decaying MHD turbulence in Minkowski spacetime. In all cases, the assumed initial magnetic energy spectrum has a simple power-law behavior, and the dissipation parameters are constants. Obviously, one should keep in mind that these results were obtained by using different initial conditions, such as initial magnetic Reynolds number, initial Γ-ratio, an so on. As it is evident from the Table, the situation is still unclear. Nevertheless, we see that convergence in the results, at least for the most recent works and for the case of helical magnetic fields seems achieved (see, in particular, Refs. [49,51,52,54,55]).

VIIIa. Evolution laws in non-expanding universe
Prior to recombination, namely for temperatures T T rec ≃ 0.3 eV [57], the equations of turbulent MHD are given by Eqs. (1) and (2) by the replacement where ν is the kinematic viscosity. After recombination, instead, and in the tight-coupling regime between ions and neutral particle species, the turbulent MHD equations are given by Eqs. (1) and (2) by the replacement [54,55] − αv → 1 where ρ i is the matter density of ions and α in is the momentum transfer rate due to ion-neutral collisions [54,55].
Although the form of the dissipation term is different before and after recombination, it is inessential to our discussion, since we work in the hypothesis of large kinetic Reynolds numbers, namely assuming that dissipation is negligible. [After recombination, the condition of tight-coupling is equivalent to assuming that the dissipative term in Eq. (138) is negligible in the Navier-Stokes equation [54,55].] Proceeding as in section IVb, we can decompose the magnetic field according to Eq. (57). In this case, the quadratic term v · ∇v can be approximated by v · ∇v ≃ v 0 · ∇ v. This term in the Navier-Stokes equation can be neglected if the condition Γ| v|/|v 0 | ≪ 1 is satisfied. In fact, comparing it with the Lorentz force, we have |v · ∇v|/|F L | ∼ Γ| v|/|v 0 |. Since we are assuming | v| ≪ |v 0 |, the above condition is certainly valid if the Γ ratio does not assume very large values (and this is the case studied in [49] and in this paper). The approximation just discussed corresponds to assuming that the small-scale components of the velocity field do not play any role in the development of MHD turbulence, at least on large scales. This means, in turn, that the Lorentz force acting on the charged particles of the fluid is responsible for the development of turbulence on large scales. In the limit of vanishing dissipation or high Reynolds numbers (see below), we have from Eq. (3) To further simplify the Navier-Stokes equation, we write ∂ t ∼ 1/τ d , giving v ≃ τ d F L , where τ d is the fluid-response time to the Lorentz force. This sort of "drag time' was first introduced in [43], and its explicit expression as a function of time was derived in [49], .
The quantity γ ∞ is a constant since the results of numerical integration of turbulent MHD equations give, asymptotically, Γ(t) ∼ 1 in the non-helical case and Γ(t) ∼ constant (eventually different from unity) in the helical case [54,55], while scaling arguments indicate that τ eddy (t) ∝ t [49]. Consequently, the turbulent MHD equations are, formally, the same as the dragged MHD equations and then the master equation for the quantity ζ diss , which rules then the evolution of the magnetic energy and correlation length is Looking at Eqs. (140) and (83), we see that the substitutions δ(0) = 1, γ(0) → γ ∞ , and a = −1 in the solutions found for the dragged case (see sections Vc, Vd, and Ve), give directly the solutions for the magnetic energy and correlation length in freely-decaying MHD turbulence.

VIIIb. Evolution laws in radiation-dominated universe
Non-helical case.-Using the results in the above subsection and taking into account Eq. (112), we easily find the evolution laws of the magnetic field in radiation-dominated universe: for the comoving magnetic field strength, and for the comoving magnetic correlation length, where v i = Γ should be replaced by B(R i )/ ρ(R i ) + P (R i ) in the notation in which the quantity ρ + P is not taken to be 1, as already noted in section Vh]. In the above equations, and κ 15 , κ 16 , κ 17 are given in Appendix A. In Eq. (143), we have assumed that κ 17 (N i v i ) −1 > 1. In the opposite case, we have approximatively R 1 ≃ R i .
Helical case.-In the helical case, we find for the comoving magnetic field strength, and for the comoving magnetic correlation length, where with q = 3(3 + p)/(4p), and κ 18 , κ 19 , κ 20 , κ 21 are given in Appendix A. In the above equations, we have assumed that is a quantity much greater than unity.

VIIIc. Evolution laws in matter-dominated universe
Non-helical case.-In the case of a matter-dominated universe, we find for the comoving magnetic field strength, and for the comoving magnetic correlation length, where we have taken into account Eq. (121). In the above equations, and κ 22 , κ 23 , κ 24 are given in Appendix A. Helical case.-In the helical case, and neglecting loglog factors, we find for the comoving magnetic field strength, and for the comoving magnetic correlation length. Here, and κ 25 , κ 26 , κ 27 , κ 28 are given in Appendix A. The above equations are valid in the limit h B ≪ 1.

VIIId. Exiting the turbulence phase
The turbulence phase prior to recombination is characterized by high values of the kinetic Reynolds number, Re = vl/ν ≫ 1 where, as usual, v and l are the typical velocity and length scale of the fluid motion. Since v, l, and ν are functions of time in an expanding universe, there could be a time when the system, initially in a turbulent phase, undergoes a transition to a viscous phase characterized by Re < 1.
After recombination, the kinetic Reynolds number, known in this case as "ambipolar" Reynolds number, has the form Re amb = vlα in X e /B 2 , where X e is the (constant) ionization fraction after recombination, and B is the typical magnetic field strength. Turbulence operates if Re amb ≫ 1, a condition that, although initially satisfied, can be (and indeed is) successively violated.
Kinetic and ambipolar Reynolds numbers.-In order to study the evolution of the kinetic and ambipolar Reynolds numbers, let as define them accurately as In a flat Friedmann universe, the kinematic viscosity re-scales as [35,54,55] ν →ν = R −1 ν, radiation-dominated era, where ν as a function of the temperature is given by [39] ν(T ) ≃ [α em log(1/α em ) T ] −1 for T ≫ m e and ν(T ) ≃ (n e σ T ) −1 for T ≪ m e , with α em being the fine structure constant, n e ∝ T 3 the electron density, and σ T the Thompson cross section. The momentum transfer rate α in evolves with the expansion parameter as α in ∝ R −3 [54,55], so that in "tilde" variables we haveα in = R 3/2 α in ∝ R −3/2 . Radiation-dominated universe.-In radiation-dominated universe, the asymptotic expansion of the kinetic Reynolds number is in the non-helical case, and in the helical case.
Matter-dominated universe.-In a matter-dominated universe, instead, we find and in the non-helical and helical cases, respectively. Regarding the ambipolar Reynolds number, we have Re amb =ṽ rms ξ vαin X e /B 2 rms in expanding universe. Using the fact that in a turbulent phaseṽ rms ≃τB 2 rms /ξ B (see section VIIIa), we get Re amb ≃τα in X e . Sinceτ ∝ ln R in a matter-dominated universe, we finally get in both non-helical and helical cases.
Excepted the case of helical magnetic fields in radiation-dominated universe, we see that the kinetic and ambipolar Reynolds numbers are decreasing functions of time (we remember that we are assuming p > 1). Hence, there will exist two times t exit,1 and t exit,2 when Re(t exit,1 ) = 1 and Re amb (t exit,2 ) = 1, thus defining the end of turbulence. 3

IX. ON THE CUT-OFF OF THE INITIAL POWER SPECTRUM
In section Va, in order to get finite results for the magnetic energy, we introduced a Gaussian cut-off, ∝ e −k 2 , for the initial magnetic energy spectrum. In the literature (see, for example, [40,41,48]), it is sometimes used a cut-off of the form e −k 4 , sharper than the Gaussian one. Therefore, it is of some interest to consider the initial magnetic spectrum instead of Eq. (70). e −k 4 cut-off: Non-helical case.-Following the same arguments in section Vb, we find for the magnetic energy and correlation length instead of Eq. (77) and (78), respectively. In the limit of large ζ diss , which corresponds to large values of the time, we have the asymptotic expressions where we used the asymptotic expansion of the confluent hypergeometric function U (a, b, z), namely U (a, b, z) ∼ z −a for z ≫ 1 [59], valid in the case a > 0. The master (differential) equation for ζ diss is given by Eq. (83), with the formal expression of δ(0) given by Eq. (85). The solution of this master equation gives Eqs. (91) and (92), with the constants c 2 and c 3 replaced by respectively. Also, we find straightforwardly that the expression of δ(0) as a function of the index p is instead of Eq. (105).
e −k 4 cut-off: Helical case.-Following again the analysis in section Vb, we find that in the helical case, the expressions for the magnetic energy, helicity, and correlation length, take the form in the limits ζ diss ≫ 1 and ζ α /ζ diss ≫ 1. Note that they differ from the case of a Gaussian cut-off for the introduction of the coefficients c 13 , c 14 , c 15 (given in Appendix A), not present in Eqs. (186), (187), and (188). Therefore, following the analysis in the Appendix B, we find that the expressions for the magnetic energy and correlation length are given, respectively, by Eqs. (97) and (98), with the replacements Other cut-offs.-We conclude this section by observing that very similar results are obtained if the cut-off of the initial spectrum is even sharper than the e −k 4 one, as for example in the case In general, what we observe is that changing the form of the cut-off, just changes the form of the coefficients c 2 , c 3 , c 7 , c 8 , and the expression for δ(0) as a function of p. This, in turn, changes the form of the κ's coefficients which enter in the expression for the magnetic field intensity and correlation length in the case of interest of an expanding universe. However, the power-law behaviors, and logarithmic corrections in the helical case, do not depend on the form of the cut-off. Finally, we note that the authors in [51], as already explained in section Va, found that soon after a first-order phase transition, the magnetic spectrum has the form of a Batchelor spectrum E B ∝ k 4 at small wavenumbers. Also, they found that during a turbulent phase, a Kolmogorov spectrum of the type E B ∝ k −5/3 develops at intermediate wavenumbers, k B < k < k diss , where k B = 2π/ξ B and k diss = 2π/ √ ηt, with √ ηt being the dissipation scale. This means that, if a dragged phase follows a turbulent phase (as it could be in a cosmological context [54,55]), a more realistic form of the initial magnetic spectrum in dragged phase could be given by where λ 1 , λ 2 , and λ 3 are constants, k B = 2π/ξ B (t i ), k diss = 2π/ √ ηt i , and t i the initial time. So, in this case the cut-off is represented by a Kolmogorov spectrum up to the initial dissipative scale, 2π/k diss , followed by a Gaussian cut-off due to resistivity.
In the light of the above discussions, however, we expect that even in this more realistic and complicated case, the main results we found, namely the power laws for the magnetic intensity and correlation length (corrected by logarithmic factors in the helical case), do not change.

X. DISCUSSION
In this section, we discuss our results and write down the final and simplified expressions for the magnetic field intensity and correlation length in a Friedmann universe, which are relevant when studying the evolution of a primordial, phase-transition-generated, cosmic magnetic field.  (175) and (176) of a primordial magnetic field in dragged and turbulent magnetohydrodynamic phases, in radiation (RD) and matter (MD) eras. For the dragged phase, it is indicated the free-streaming species which determines the drag coefficient. The parameter p is the index of the initial magnetic power-law spectrum, EB(k, ti) ∝ k p .

Phase
Era Streaming particle ̺1(p) ̺2(p) ̺3(p) Simplified evolution laws.-We saw, in the previous section, that the κ coefficients entering in the expressions for the magnetic field and correlation length depend on the choice of the cut-off of the initial magnetic spectrum. In the case of magnetic fields generated in primeval phase transitions, the exact form of this cut-off is not precisely known, and this introduces a factor of arbitrariness in the evolution laws of B and ξ B . Fortunately, since the κ coefficients are all of order unity (see section IX and Appendix A), this factor of arbitrariness can be safely neglected when studying the evolution of a phase-transition-generated magnetic field.
If on the one hand the κ coefficients depend on the initial cut-off, on the other hand the logarithmic factors we found in the case of helical magnetic fields are a general consequence of the fact that the magnetic helicity is a quasi-conserved quantity in magnetohydrodynamics (see Appendix B). However, these factors introduce just a small correction in the power-law evolution of B and ξ B and then, in first approximation, can be neglected.
According to the above discussion and given an initial magnetic field with spectrum of the form E B (k, t i ) ∝ k p , the (physical) magnetic field intensity, B, and the comoving magnetic correlation length, ξ B , evolve as a function of the temperature T approximatively as and respectively, where B i = B(T i ) and ξ B,i = ξ B (T i ), and we took R ∼ T −1 [neglecting small corrections due to g * ,S (T )]. We remember that N i is the initial number of magnetic domains per horizon length, namely N i = d H,i /ξ B,phys,i , where d H,i and ξ B,phys,i are the length of the Hubble horizon and the physical magnetic correlation length at the initial time. The value of the bulk velocity at the onset of the particular regime, v i , is not an independent parameter in both dragged phase and turbulent non-helical phase, being related to B i , ξ B,phys,i , and α i by the relations In equations (175) and (176), we quoted the expressions for B and ξ B in the helical case. In the non-helical case, those expressions are still valid, the inverse cascade phase being absent. The transition temperatures T 1 and T 2 can be easily found by matching the expressions of B (or, which is the same, those for ξ B ) in two consecutive phases. We summarize our results on the evolution laws for the magnetic field and correlation length in Table 2.
Comparison with the Banerjee and Jedamzik's results.-Let us now compare the above results with those of [54,55], which are, to our knowledge, the only ones in the literature which consider other than the case of freely-decaying magnetic fields, also the case of magnetic fields in dragged phase. To this end, we introduce the so-called Alfvén eddy turnover time, τ A , as τ A = ξ B,phys /B. We can then write the quantity N i v i as where we remember that in turbulent regime, Γ i is of order unity for the non-helical case, but can be different from one in the helical case. Banerjee and Jedamzik [54,55] assume (without a full justification) that the initial correlation length is not a free parameter but it is determined, at the initial (cosmic) time t i , by the equality of the Hubble rate H ≃ t −1 and the Alfvén eddy turnover rate 1/τ A . Namely, in their model we have In this case, and taking Γ i of order of unity (as in [54,55]), we get With the choice in Eq. (179), it is easy to check that our formulas agree with those of [54,55] [with the exception of Eqs. (62) and (63) of [54,55], in which the exponents 3/(2 + n) and 3n/(2 + n) should be 4/(2 + n) and 4n/(2 + n), respectively.] For example, for non-helical dragged MHD in the case of neutrino free streaming we have, excluding factors of order unity, where G F is the Fermi constant, m Pl is the Planck mass, and we used the facts that t i ∼ m Pl /T 2 i (in radiationdominated era) and α i ∼ G F T 5 i (for the neutrino case). Equations (181) and (182) agree, respectively, with Eqs. (48) and (47) of [55]. [Let us observe that in the notation of [55], p = n − 1 and B(T ) ∝ r(T ) 1/2 T 2 .] Inertial range.-According to the "Kolmogorov hypothesis"(see, e.g., [24]), the cascade of energy in k-space is a quasi-local process independent on the particular scale considered (although this is strictly possible only for scales larger than the dissipation scale and smaller than the integral scale). The k-space interval where the Kolmogorov hypothesis may apply is called "inertial range". Here, the magnetic spectrum is expected to decay as E B (k, t) ∝ k β , with β being negative. The value of β could in principle depend on the particular phase (dragged or turbulent) and be different for the non-helical and helical cases. Indeed, the pure Kolmogorov spectrum, β = −5/3, has been observed in [51] for the turbulent non-helical case, while [54,55] found β ≃ −2 both in helical and non-helical dragged case.
In our case, we do not find any inertial range where the magnetic spectrum decays as a power law. This is probably to be ascribed to the mean-field approximation, which neglects the small-scale fluctuating part of the velocity and magnetic fields (see section IVb). It is the (quasi-linear) interaction of these small-scale fields with the corresponding large-scale ones which "opens" this particular range in k-space. However, if the inertial range is sufficiently narrow, we do not expect significant modifications to the evolution laws we found for B e ξ B . Indeed, and interesting enough, Ref. [54,55] found that the inertial range is almost absent in the non-helical turbulent case, and that on intermediate scales the magnetic spectrum is more consistent with an exponential law than with a power law. This seems to be in agreement with the results we found in mean-field approximation, namely with Eq. (94).

XI. CONCLUSIONS AND OUTLOOK
Conclusions.-The presence of large-scale magnetic fields in the universe is still an open and unsolved mystery of modern cosmology. A plethora of mechanisms able to generate cosmic magnetic fields in the early universe, as for example during inflation or primeval phase transitions, or during more recent eras by astrophysical mechanisms have been proposed since the first attempt put forward by Harrison in 1970 [11].
If on the one hand, astrophysical mechanisms seem to be ruled out by recent observations of large-scale magnetic fields in cosmic voids [17], on the other hand inflationary mechanisms are able, at least in principle, to explain the magnetization of the universe. However, these latter mechanisms repose on the use of nonstandard physics, an exception being the mechanism recently proposed in [22]. Nevertheless, if the scale of inflation is considerably below 10 16 GeV, also this mechanism fails to explain galactic and galaxy cluster magnetic fields.
If, instead, presently-observed magnetic fields originate from cosmological phase transitions (such as QCD or electroweak phase transitions), another unsolved question arise: how have these magnetic fields evolved from the time of they generation until today? The answer to this question can be given only in the framework of magnetohydrodynamics.
In this paper, we have studied the evolution of phase-transition-generated magnetic fields coupled the to primeval plasma by solving, analytically, the magnetohydrodynamic equations in mean-field approximation. The reduction of full MHD equations to simpler equations due to mean-field approximation causes a lost of information about the transfer of magnetic energy at intermediate scales (corresponding to the inertial range). Nevertheless, the main characteristics of an evolving magnetic field are preserved.
In particular, we have analyzed the decay of primordial magnetic fields, both in radiation and matter eras, in the two regimes which are relevant in a cosmological context, namely the turbulent and viscous free-streaming regimes. During a viscous free-streaming phase, dissipation is ruled by a drag term in the Navier-Stokes equation, the drag term depending on the free-steaming particle species, namely neutrinos, photons or hydrogen atoms.
Our results can be summarized as follows.
If the initial magnetic field is not helical, then the magnetic helicity remains null during the evolution of the system. There is an initial phase in which the system is quiescent: magnetohydrodynamic effects do not operate and the characteristic quantities of the system, such as energy and correlation of the magnetic field, remain almost constant. This phase persists for a period which is proportional to the initial eddy turnover time. After that, the evolution of the magnetic field proceeds through selective decay of magnetic modes, that is there is no direct transfer of magnetic energy from small scales to large scales, but simply modes with larger wavenumbers are dissipated faster than those whose wavenumbers are small. During this process, the magnetic correlation length grows, while the magnetic energy decays in time. The evolution laws for the magnetic field depend on the initial magnetic power spectrum and on the particular regime. They are summarized in Eqs. (175)-(177) and Table II. In the helical case, the magnetic helicity is an almost conserved quantity. The system undergoes three different phases: a quiescent phase and a selective-decay phase in which the system evolves irrespective of the presence of magnetic helicity, and an inverse-cascade phase. As in the case non-helical case, the quiescent phase ends approximatively after a period of time equal to the initial eddy turnover time. Then, the system enters into a selective-decay phase characterized by a magnetic field evolution similar to the non-helical one. This last phase ends due to the conservation of magnetic helicity, favoring an inverse cascade of the magnetic field during which the magnetic energy stored on small scales is partially transferred to larger scales. This process of energy redistribution weakly depends on the properties of the initial magnetic field, so that the evolution laws of magnetic energy and correlation length do not depend on the index of the initial magnetic spectrum power law, although they are different in turbulent and free-streaming regimes [see Eqs. (175)-(177) and Table II]. Moreover, the time when the system enters into the inverse-cascade regime depends both on the form of the initial magnetic spectrum and on the amount of initial magnetic helicity.
Outlook.-Our analytical evolution laws in the different regimes and for helical and non-helical magnetic fields, are in substantial agreement with the numerical results obtained in [54,55]. Nevertheless, it is important to stress that (i) there is no yet convergence in the literature on the results of numerical simulations of MHD turbulence in the case of non-helical fields, and that (ii) the only available numerical results for the dragged phase are those in [54,55]. Therefore, it is desirable, in order to have a firm understanding of the evolution of primordial magnetic fields, to solve this disagreement for turbulent non-helical fields and, at the same time, to have an independent (and, possibly, full numerical) confirmation of the Banerjee and Jedamzik results [54,55] for the evolution of magnetic fields in dragged phase. Finally, (iii) we stress that the study in [54,55] on the properties of phase-transition-generated magnetic fields throughout the evolution of the universe assumed a linear dependence between the initial magnetic correlation length and the initial magnetic field intensity. It would be important to see if the conclusions of [54,55], that phase-transition-generated magnetic fields may directly account for galactic and galaxy clusters magnetic fields, are modified or even invalidated by relaxing the above ansatz on the initial magnetic correlation length [71].
We would like to thank M. Giovannini for useful discussions.