{'fraction_non_alphanumeric': 0.035685637704693525, 'fraction_numerical': 0.014458161328820456, 'mean_word_length': 4.345534407027818, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 15, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Primordial Weibel instability

We study the onset of vector instabilities in the post-inflationary epoch of the Universe as a mechanism for primordial magnetic fields amplification. We assume the presence of a charged spectator scalar field arbitrarily coupled to gravity during Inflation in its vacuum de Sitter state. Gravitational particle creation takes place at the transition from Inflation to the subsequent Reheating stage and thus the vacuum field state becomes an excited many particles one. Consequently this state can be described as a real fluid, and we build out the hydrodynamic framework using second order theories for relativistic fluids with a relaxation time prescription for the collision integral. Given the high-temperature regime and the vanishing scalar curvature of the Universe during Reheating (radiation-dominated-type era), the fluid can be regarded as a conformal one. The large quantum fluctuations induced by the rapid transition from inflationary to effectively radiation dominated expansion become statistical fluctuations whereby both a charge excess and anisotropic pressures are produced in any finite domain. The precise magnitude of the effect for each scale is determined by the size of the averaging domain and the coupling to curvature. We look at domains which are larger than the horizon at the beginning of Reheating, but much smaller than our own horizon, and show that in a finite fraction of them the anisotropy and charge excess provide suitable conditions for a Weibel instability. If moreover the duration of reheating is shorter than the relaxation time of the fluid, then this instability can compensate or even overcome the conformal dilution of a primordial magnetic field. We show that the non-trivial topology of the magnetic field encoded in its magnetic helicity is also amplified if present.


Introduction
Except for the Electroweak (EW) and Quantum Chromodynamic (QCD) phase transitions the epoch of the Universe prior to Big-Bang Nucleosynthesis (BBN) is practically unknown [1][2][3]. Processes such as the Reheating of the Universe after Inflation, which lead to the Radiation epoch are still arena of intense research. Since the seminal works of Albrecht, Steinhardt et al. [4], Kofman, Linde et al [5] and others, the understanding of the Reheating of the Universe evolved from an almost instantaneous process to being considered a stage in itself, with several e-folds of duration [6][7][8]. This fact leads to the question of what processes other than the decay of the inflaton might take place during that period, specially those related to the evolution of the gravitational field and primordial magnetic fields.
Among different out-of-equilibrium processes in the Early Universe, the onset of instabilities can play a determinant role during Reheating by accelerating the thermalization of the primordial radiation [6,7], induce or amplify primordial magnetic fields [9][10][11][12][13], and perhaps help to trigger phase transitions as the EW one.
In general the research about the first stages of the Universe after Inflation is carried out using non-equilibrium quantum field theory tools [14], but recent developments in the study of real relativistic fluids [15] began to pave the road toward the treatment of the matter fields in the early Universe within relativistic hydrodynamic formalisms [16,17].
One of the main unsolved cosmological and astrophysical problems is the origin and evolution of magnetic fields that permeate the Universe [72]. Magnetogenesis mechanisms that operate at the end of Inflation, as e.g. the one of Refs. [9,73] consider that the induced field is subsequently diluted by the cosmological expansion. In this manuscript we wonder what processes might act at the same epoch and compensate such dilution. Using relativistic hydrodynamics, in this manuscript we want to investigate one such processes, namely the Weibel instability. As a concrete hydrodynamic model, we take the results of Ref. [74] where using a Grad prescription for the 1pdf, the relativistic Weibel instability was naturally described.
Weibel instability can be important for magnetogenesis during the primordial QCD phase transition as was envisaged long time ago, when it was proved (both numerically and analitically) that if a cold partonic beam passes through a quark-gluon plasma a beam-driven Weibel instability develops and magnetic fields are excited (see section 4.4 of Ref. [75] and references therein). The magnetic fields induced by the instability can be quite intense for relatively small anisotropic velocity distributions, and of coherence length equal to the scale of the patch where the phase transition occurs.
Since Weibel's seminal work on instabilities due to anisotropies in the velocity distribution function of gases [76,77], relativistic generalizations were developed in the framework of kinetic theory [78][79][80], and hydrodynamic descriptions for the non-relativistic case were worked out not long ago [81][82][83]. In addition plasma instabilities related to the quark-gluon plasma were studied in [84,85]. In the model of Ref. [74], the anisotropies in the velocity distribution are encoded in a new tensor fluid variable that represents an anisotropic pressure. This SOT hydrodynamic framework has already been used to describe the primordial plasma in the very Early Universe and, in particular, to capture its interaction with the primordial gravitational waves [69,70].
We consider a cosmological scenario where an inflationary period occurred, and during which a charged spin-0 field coupled to gravity was in its vacuum state. At the transition between Inflation and Reheating, which for simplicity we assume instantaneous, the vacuum state becomes an excited, highly populated state [87,88]. Due to the high symmetry of the considered stages of the Universe, the quantum gas is also highly symmetric and consequently all mean values of vector quantities as e.g., momentum flux, electric four-current, energy flux, will be vanishing. Quantum fluctuations on the other side, build non-vanishing rms deviations. These quantum fluctuations become hydrodynamical fluctuations in the large occupation numbers regime. To obtain the corresponding rms hydrodynamical fluctuations we apply the Landau prescription, whereby the quantum Hadamard two point functions for the electric four-current J µ and spatial components of the energy momentum tensor T µν of the scalar field, averaged on a scale of interest L D , are equated to the corresponding two-point correlation functions of a real relativistic plasma. The temperature of the fluid is estimated from the vacuum expectation value (VEV) of the scalar field energy density T 00 together with a consistent relation between number and charge densities. This estimation and the fact that the scalar curvature of the universe during Reheating is zero, will indicate whether the fluid can be regarded as conformal or not: if the temperature is higher than the effective field mass then we shall consider that it is indeed conformal. This entails that phenomena related to the trace of the energy momentum tensor, such as bulk viscosity, will be neglected in the present work.
To be more precise, the diagonal terms in the spatial part of the EMT of the gas of scalar particles contain information about the isotropic pressure (or hydrodynamic pressure) as well as about the anisotropic (or non-hydrodynamic) contributions. To estimate the later we shall build a Hadamard two-point function of the tensor T ij tls = T ij − (1/3)g ij T k k , i.e., from the spatial components of the EMT without the hydrodynamic components. After filtering over a scale of physical interest we equate the resulting expressions to corresponding fluid variables. Those variables cannot be the usual hydrodynamic ones because they do not sustain anisotropic pressures. In SOTs they are readily identified as higher moments of the one particle distribution function.
We assume a de Sitter inflationary period after which a Reheating processes took place to finally establish the bulk radiation field. According to Refs. [7,8], the duration of those stages depends on the equation of state. The expansion rate is the same as for a radiation dominated universe almost immediatly after the end of Inflation while the radiation field takes several e-folds to establish. This means that during that period the charged scalar field and the bulk fields can be considered as decoupled, thus allowing the use the one-fluid formalism developed in Ref. [74].
Whether a Weibel instability will occur on any given domain depends on several occurrences of essentially a random nature. The most important conditions are of course enough pressure anisotropy and enough charge excess. The statistics of pressure anisotropy is discussed in Section 3.1. It depends on two independent Gaussian parameters whose variance is given in equations (3.30) and (3.32), and on the size of the domain we are looking at in a way given by Eq. (4.14); similarly, charge excess is a third Gaussian parameter whose statistics is discussed in Section 3, see Eqs. (3.26) and (4.20). Fluctuations are larger for small domains and die out for larger ones; of course the Universe as a whole is isotropic. For concreteness, we look at domains which are about ten times the size of the horizon at the beginning of Reheating. Since there are a huge number of them within our own horizon, in a large number of domains the actual anisotropy parameters and charge excess will be about or larger than their variances, which we take as an appropriate order of magnitude estimate.
Given a domain where pressure is sufficiently anisotropic, the anisotropic pressure tensor will pick up a preferred direction, which we may identify with the direction of a particle beam, see the discussion in Section 2.3 and Fig. 1. This direction is essentially equally distributed on the unit sphere. The modes which will be most amplified are those which propagate in an orthogonal direction to the "beam", but a perfect alignment is not necessary, see Fig. 6. This gives enough leeway that an instability will likely occur in a significant fraction of the suitable domains. For the actual strength of the instability see Fig. 5.
The final issue is for how long the instability will be operative. The average pressure anisotropy, which we take as a homogeneous background within a given domain, will decay on time scales of the order of the relaxation time of the fluid. The relaxation time is essentially determined by the fluid temperature, see Eq. (5.6) and Fig. 4, and the fluid temperature depends in a highly nonlinear way on the domain size and the coupling to curvature, see Figs. 2 and 3. On the other hand, as reheating progresses the ordinary radiation background in the Universe will build up, and since our fluid is charged, it will efficiently thermalize with it. So we need two conditions, the fluid must be close enough to free streaming, and moreover the reheating era must last less than the relaxation time. We discuss these conditions in Section 5. If these conditions are met, the Weibel instability may successfully compensate the dilution of the magnetic field due to the expansion of the Universe, see Fig. 7.
In conclusion, although whether a strong enough Weibel instabily actually happens in a given domain is a random phenomenon, the probabilities are significant enough that we may expect it will actually occur in a large number of domains.
To summarize the results, we found that for different values of the coupling of the scalar field to gravity the temperature of the fluid is high enough to consider it as conformal during Reheating. For the scenarios described in the above paragraphs, we find that for several Reheating models [93] the instability compensates de dilution due to the expansion at the end of Reheating. This sets new initial values of primordial magnetic fiels for their subsequent evolution, The manuscript is organized as follows. In Section 2 we briefly review the derivation of the set of tensors and equations for the SOT description of the primordial plasma made in Ref. [74] and analyze its conformal invariance. We linearize the equations around an anisotropic background pressure and for a generic propagation direction, and obtain a dispersion relation for the vector modes, which now posses two sets of solutions associated to two orthogonal polarizations. In Section 3 we build the simplest matching relations between the fluid and the field two point functions, from which we read the VEVs we need to compute.
In section 4 we provide the expressions for T µν and and the four-current J µ . (To make the mansucript as self-contained as possible, in Appendix A we give a short review of charged scalar fields in curved spacetime.) We plot the result of the numerical calculation of the VEV T 00 , which shows that the energy density is weakly dependent on the small mass of the scalar field but sensitive to the value of the coupling to the gravitational field. We provide the definitions for the noise kernels to be used to compute the rms values of the electric charge and anisotropic pressure, and their general expressions in terms of the field modes as well. We compute the filtering on a scale L D of the resulting expressions of those kernels finding that both quantities depend on the averaging scale as L −3 D . Regarding the dependence on the coupling to gravity, we found that the anisotropic pressure depends very weakly on ξ, while for the charge excees the dependence is stronger, going as ξ −1 . These two are the main parameters that determine the intensity of the instabilites.
In Section 5 we apply the results of the previous section to evaluate the fluid temperature and to study the linear unstable modes from the dispersion relation found in Section 2. We confirm that the temperature is high enough to consider the plasma as conformal. For both polarizations the main instability occurs for propagation orthogonal to the direction where the anisotropic pressure is maximal. We continue the analysis taking into account only the largest solution to the dispersion relation. As stated in previous paragraphs, we focus on Reheating scenarios whose duration is shorter than the relaxation time of the fluid. The reasons for this assumption are several: (1) The background radiation field to which the scalar plasma would couple is not settled yet, thus we are able to work with a one-fluid formalism; (2) since the effects of the collisions strongly manifest after the relaxation time τ , we regard an almost collisionless plasma. Under these conditions the Weibel instability is operative during the entire period.
We find that for some models of Reheating, there is a range of values of ξ for which the instability overcomes the dilution by the cosmological expansion at the end of this stage. The associated scales are of the order of the particle horizon at that time. We end this section studying how the topology of the field, encoded in the magnetic helicity, might be affected by the instability. We find that it is also amplified and that the amplification depends on both unstable frequencies (not only on the largest). This result implies that the chances of survival of an helical primordial magnetic field are improved.
Section 6 contains the main conclusions on this work. We work with signature (−, +, +, +) and natural units = k B = c = 1.

Fluid Description of the Primordial Plasma
In this section we shortly review the hydrodynamic description of a conformal plasma developed in Ref. [74], where we wrote down a minimal SOT to describe a real relativistic magnetized plasma, starting from Vlasov-Boltzmann equation and adopting a Grad prescription for the non-equilibrium distribution function. In this section we particularize the results to a fluid in an expanding FRW space-time and re-obtain the linear equations for the vector variables considering propagation along an arbitrary direction. There are two sets of solutions of the dispersion relation that correspond to the two possible polarizations of the modes and hence two sets of possible unstable momenta.
The general form of the tensors that describe the magneto-hydrodynamics is with s = 0, 1 and the general conservation law reads where the notation (µ i ) means that µ i is excluded, the tensor I µ 1 ···µn s represents dissipation and it is given by The fundamental tensors we work with in this paper are the charge current J µ , the energy momentum tensor T µν , two non-conserved tensors, A µν and A µνσ , whose divergencies provide the equations for Ohmic dissipation and viscous stresses respectively, and the Maxwell tensor F µν . The number current N µ is automatically conserved in a conformal theory.
We shall consider the evolution of the real relativistic plasma in a spatially flat Friedman-Robertson-Walker (FRW) Universe. In conformal coordinates x µ = (η, x), the metric tensor is g µν = a 2 (η)g µν . It is straightforward to show that the hydrodynamic theory under study is conformally invariant. The conformal transformations that apply in this case are T =T /a(η), α =α, u µ =ũ µ /a(η), τ = a(η)τ , ζ µ = a −2 (η)ζ µ , ζ µν = a −2 (η)ζ µν and F µν = a −4 (η)F µν . Moreover, to write down dimensionless equations we choose the Hubble constant during Inflation (H) as the fiducial dimensionful energy scale. The dimensions of each variables are For the coordinates we apply the transformation x µ → x µ /H, so each derivative (either time or spatial) will be accompanied by a factor H.
Explicitly we haveÂ where 14) Clearly exprs. (2.13) and (2.14) satisfy the constraint which will be used in Section 5 to evaluate the plasma temperature. The electromagnetic tensor isF whereη µνρσ is the four-dimensional totally antisymmetric tensor. The moments of the collision integral that appear in the r.h.s. of the conservation equations are [74] The conservation equations of the theory areĴ µ ;µ = 0 (2.20) is the transverse projector and is the symmetric, double transverse and traceless projector. The non-equilibrium equations (2.22) and (2.23) must be projected in order to provide the right number of supplementary equations [74].

Linear Analysis
To analyze the evolution of linear perturbations we propose the following decomposition: where η µν is the unperturbed Minkowski metric and U µ 0 is the four-velocity of the fiducial observers from which we define the electric and magnetic fields. Furthermore U µ 0 defines the time direction and the time derivative through the partial derivative along U µ 0 as U µ 0 X β ,µ = dX β /dη. We also consider U µ 0 v µ = 0 and U µ 0 U 0µ = −1, and construct the zeroth order spatial projector ∆ µν = U µ 0 U ν 0 + η µν . This four-velocity for comoving observers is U µ 0 = (1, 0, 0, 0). We use U µ 0 and ∆ µν in order to split the system into its temporal and spatial projections. Given the relations (2.28) and (2.29) . For the Maxwell tensor we write

Solving the Vector Sector
To solve eqs (2.20)-(2.25) we propose a solution that varies as exp ση + ik.r , we extract the vector sector of the theory from the equations for the momentum, conduction current and Ampère's law and the equation for the viscous stress tensor. We only consider the terms which contain vector degrees of freedom, disregarding the first order scalar quantities, in order to capture the simplest dynamical behavior of the vector sector. We first write down the equations in a coordinate system in which the background anisotropic pressure tensor reads Z ij 0 = diag(Z 0 + Z * , Z 0 − Z * , −2Z 0 ) while the spatial propagation direction k remains arbitrary. From eqs. (2.10) and (2.31) we see that Z 0 and Z * must satisfy −5/4 ≤τ (Z 0 ± Z * ) andτ Z 0 ≤ 5/8 for the pressure to be non-negative . The notation Z 0 ± Z * simultaneously encodes two quantities Z 0 + Z * and Z 0 − Z * and the inequality must be fulfilled for both. We adopt the more stringent bound −5/4 ≤τ (Z 0 ± Z * ) ≤ 5/8. Defining Q = q 10 /q 20 the set of vector equations now reads The expression between rounded brackets in the last line of eq. (2.38) still couples the different components of the magnetic field but it is easier to deal with. We analyze this issue and the decoupling of the different components of the vector quantities in the following subsection.

Fixing the propagation and background orientations
The system of equations (2.36), (2.37), (2.39), (2.40) and (2.38) which describes the linear dynamics of the vector sector is written using a coordinate system in which Z 0ij = diag (Z 0 + Z * , Z 0 − Z * , −2Z 0 ) and, without loss of generality Z 0 , Z * > 0 and Z * < 3Z 0 , and the propagation vector k remains arbitrary. Since eq. (2.38) mixes and couples the different spatial components of the magnetic field, here we show how to decouple them. We begin by setting a rotation transformation O to rotate and fix the propagation direction such that k → k = O k with k = (0, 0, k). Let θ and ϕ be the usual polar and azimuthal angles which describe the arbitrary propagation directionk = (cos ϕ sin θ, sin ϕ sin θ, cos θ). In fact we have to transform all the other vectors and tensors, for instance the anisotropic pressure tensor transforms as Hereafter a 'prime' denotes a transformed quantity. We are only interested in the components of the dynamical variables perpendicular to the propagation direction. Since for the new rotated quantitesk =ẑ , we consider just the x − y directions. It turns out that the tensor which mixes the directions in eq. (2.38) now reads Z 0ij −k ik l Z 0lj . This is not a symmetric tensor but it is straightforward to see that the x − y block does satisfy the symmetry condition, so we can diagonalize it to obtain the eigenvalues with the corresponding normalized eigenvectors For the complete tensor Z 0ij −k ik l Z 0lj there is a third eigenvalue λ 0 = 0 with its corresponding left eigenvectork . We can build a new rotated orthonormal basis made of {ε + ,ε − ,k }. Therefore considering only the components of the vector dynamical variables which belong to the x − y plane, the mixing term in eq. (2.38) can be rewritten in the new basis as a diagonal tensor, namely (2.43) Figure 1. Graphical representation of the propagation directionk and the transverse components ε + andε − of the magnetic field in an extreme case where we consider a particle beam propagating in the directionÎ perpendicular tok =ẑ. To accomplish this condition we fix the orientation of Z 0ij (and thus that ofÎ ) by fixing the angles defined in Section 2.3 such that ϕ = π/2 and 0 ≤ θ ≤ π.
where λ (i) denotes the eigenvalues λ ± corresponding to the directionsε ± respectively. The brackets around the subindex (i) indicate that there is no sum over i. The eq. (2.38) for thē ε ± directions then becomes Finally, the system of equations for the vector dynamical variables is composed by (2.36), (2.37), (2.39), (2.40) and (2.44). These relations must be understood as equations for the rotated (primed) quantities in the directionsε + andε − of the new basis.
In figure 1 we show a graphical representation of the vectorsk ,ε + andε + regarding a concrete example as we explain in the next paragraphs.
Since we define Z 0ij = diag(Z 0 + Z * , Z 0 − Z * , −2Z 0 ) with Z 0 , Z * > 0 and Z * < 3Z 0 , the expressions for the distribution function (2.2) and (2.3) determine that most of the particles propagate through the directionx. Considering an extreme case we could have a propagating beam in the directionÎ withÎ x. We express all our results for rotated (primed) quantities, in this way the particle beam propagates through the directionÎ = OÎ = Ox, where O is the matrix rotation defined at the beginning of Section 2.3, and the magnetic field propagates in k = Ok =ẑ.
To be concrete let us consider that the propagation direction of the magnetic fieldk is perpendicular toÎ , for instance we take ϕ = π/2 and 0 ≤ θ ≤ π. In this casek =ẑ and I = −ŷ for the entire range of θ. Moreover the largest unstable component of the magnetic fieldε − (see Section 2.4 and Figure 5) corresponds toε − =k ×Î =x, while the component ε + =Î = −ŷ as we observe in Figure 1. Furthermore in this scheme, the largest unstable growing component of the magnetic fieldε − coincides with the growing component of the magnetic field related to transverse plasma waves reported in [77].

Dispersion Relation
In this section we compute the dispersion relation of the system made up of eqs. (2.36), (2.37), (2.39), (2.40) and (2.44). We consider the dynamical equations for the vector components related to the directions of the eigenvectorsε ± . All the vector and tensor quantities correspond to the rotated (primed) ones.
We define s = στ , κ = k τ , = πe 2τ 2 q 20 Υ 2 0 and then the sought relation reads 0 = s 5 ± + 2s 4 ± + 1 + where subindex +(−) corresponds to the dynamics of the vector components along the direction of the eigenvectorε + (ε − ). Observe that we have two sets of solutions, each corresponding to either + or −. Instabilities emerge as long as there exists a positive real solution to the equation (2.45) and/or a complex solution with positive real part. In Section 5 we discuss the explicit positive roots of (2.45) with fixed parameters corresponding to the cosmological scenario to be described in Section 4. So far we set the magneto-hydrodynamic framework to study the primordial Weibel instability. We must now obtain the values of the different parameters of the fluid from a concrete particle physics model.

Building the Background
We want to study the evolution of vector magneto-hydrodynamic perturbations of a scalar plasma during Reheating, i.e. the epoch between the end of Inflation and the establishment of the bulk radiation field.
Let us consider a domain of size L D where the state of the fluid may be regarded as mostly homogeneous, though anisotropic. The values of the hydrodynamic variables will be different for each domain and therefore we can only regard them as stochastic. Our line of work is to find the first and second moments of these variables and then assume that we are looking at a typical domain. From Section 2 we see that the fluid variables are the thermal potential α, the temperature T , the four-velocity u µ and the dissipative currents Z µ and Z µν . For the fluid under consideration, it is natural to assume that Z µ = Z µν = 0. For the other variables we have α = α + α (s) , T = T + T (s) and u µ = u µ + u µ (s) , with ũ µ = ũ 0 U µ , where U µ = 1, 0 . Using the constraintsû 2 = −1 andû µ Z µ =û µ Z µν = Z µ µ = 0 we obtain The dimensionless conserved conformal currents are eqs. (2.8) and (2.10).
It is clear that the quantum state of the field is rotationally invariant and so will be the expectation values derived from it. Moreover, the initial state is charge-neutral, and the gravitational particle creation process does not pick a preferred charge sign. Consequently the expectation values must be even under charge reversal. Let us callĴ µ (Qλ) andT µν (Qλ) the Heisenberg current and the EMT averaged over a scale λ. We consider the simplest matching procedure, namely the Landau prescription: subindex S means stochastic, subindex Q stands for quantum mechanical and subindex L D indicates spatial average over a region with typical size L D .
We work under the hypothesis that the fluid is conformal and so the trace of the energy momentum tensor may be neglected. From the point of view of the formalism it is most convenient however to subtract it explicitly before computing the quantum expectation values. Neglecting the trace of the energy momentum tensor means we cannot include the effects of bulk viscosity in our analysis, we intend to address this issue in future works. To reduce the complexity of the matching we assume that whenever a stochastic variable has a non-zero expectation value its fluctuations are negligible. This means that we shall treat the dimensionless temperature Υ andũ 0 as true (but L D -dependent) constants. In this way, from eqs. (2.10) and (2.29) we obtain our first nontrivial equation and the symmetry requirement By the same logic, we shall not consider these quantities when looking at second moments. The matching conditions we need then are For the last correlation, by demanding symmetry we can write Observe that due toũ 0 (s) = 0, the constraints on the hydrodynamic variables are v (s)i Z i S = v (s)i Z ij S = 0, while Z i i = 0 identically, not just in expectation value. We shall not consider stochastic components in Z 0 nor in Z 0i . The previous matching relations give (3.14) from whose trace we determine The current-EMT correlation vanishes because q 1 S = q 1 v i (s) S = q 1 Z i S = 0. From the momentum density correlation we get which determines the rms value of the charge imbalance.

For the EMT-EMT correlation, we begin by observing that under the assumption of Gaussian fluctuations we have
Therefore the trace-free velocity correlation is We therefore make the assumptions v i (s) Z jk S = 0 (3.20) and get To end this section we must mention that it is straightforward to evaluate the r.h.s. of correlation (3.9) using the corresponding components of expr. (4.3) together with the mode functions (A.15). It is obtained that effectively

(3.23)
Therefore to study the instability we can consider only with N ij ij given in expr.

Anisotropic Background Pressure
In a certain coordinate system the anisotropic background pressure tensor reads We must find the statistics for Z 0 and Z * , i.e., we must calculate Z 2 0 S , Z 2 * S and Z 0 Z * S . We begin by asking the statistics to be invariant under the change of sign of Z * , which is just a relabelling of the coordinates, and then it must be that Z 0 Z * S = 0. We also demand that all representations have the same statistics. If we define then Z ij remains unchanged (modulo a rotation of coordinates). Demanding that the statistics is the same means Z 2 0 S = Z 2 0 S and Z 2 * S = Z 2 * S . With the above correspondence between representations we get that a necessary condition to have the same statistics is

Charged Scalar Field in a Flat FRW Universe
As stated above, we assume the presence of a spectator charged scalar field coupled to gravity, which is in its vacuum state during Inflation. At the onset of the subsequent epoch to Inflation the vacuum becomes a particle state and in consequence a scalar plasma is formed. As the quantum state is rotationally invariant and the gravitational particle creation process is insensitive to the charge signs, the mean anisotropic pressure of the plasma and the charge excess will both be zero. However due to quantum fluctuations the rms deviations from those mean values are non-vanishing in any finite domain. As stated above the details of the calculation to find the different expressions presented in this section can be found in Appendix A.
A charged scalar field is mathematically described by the complex fields [86] Φ(x µ ) and Φ † (x µ ), which can be written in terms of real fields as Φ( where Φ 1 and Φ 2 commute and satisfy the Klein-Gordon equation for a massive real spin-0 field arbitrarily coupled to gravity. The Energy-Momentum tensor (EMT) in curved spacetime of each real component of the charged scalar field of dimensionless mass µ = m/H and coupled to gravity through a coupling constant ξ is [87,88] The symbol R µν is the Ricci tensor and R the scalar curvature. The values ξ = 0 and ξ = 1/6 correspond to the minimal and conformal coupling respectively. Rescaling the scalar field as where η µν is the Minkowski metric tensor. The electric four current is [86] After replacing (4.2) for each field we have that J µ = a −2J µ with (4.5)

Computing the Energy Density and Noise Kernels
To compute the energy density and the noise kernels of the quantum scalar gas needed to build the fluid description, we consider that during Inflation the scalar field is in its vacuum state and after the transition to Reheating only adiabatically regularized modes will be present [88]. This means that observers in the local vacuum state will detect a bath of particles built up from super-horizon modes.
In Section 3 we construct the correspondence between a quantum gas and its fluid description assuming the simplest possible matching procedure. We only need three functions in order to set the background of the hydrodynamic theory, namely the VEV T 00 Q and the traces {T ij (r 1 ),T ij (r 2 )} Q and {Ĵ 0 (r 1 ),Ĵ 0 (r 2 )} Q related to the anisotropic pressure and the electric charge excess respectively. In the next subsection we compute these functions.

Computing T 00
In this section we evaluate the VEV ofT 00 at η = 0 from expr. (4.3). Recall that the factor 2 is due to the fact that a complex scalar field consists of two real fields, and expr. (4.3) corresponds to only one of those fields. Replacing the expansion (A.11) and evaluating the expectation value in an adiabatic vacuum state we obtain We numerically integrate this expression and in Figure 2 we show the value of T 00 as a function of ξ for several values of the field mass. In the case of small masses µ 1 this VEV is insensitive to µ, but it shows a strong dependence on the coupling to gravity ξ.

Anisotropic pressure noise kernel N ij
ij (x, y) In this section we compute the traceless part of the pure spatial sector of the energymomentum tensor which carries on the information of the anisotropic pressure. The corresponding noise kernel is the vacuum expectation value of the contracted symmetric two-point function Replacing the Fourier mode decomposition of the scalar field we obtain To evaluate the momentum integrals we can retain only the leading term in the expansion of the scalar field mode functions (see Appendix A) and write (4.9) We must now apply a spatial average over a scale of interest L D to filter the kernel N ijij (r). We consider a top-hat window function and write with r c ∼ 1 the coordinate at which both parts of N ijij (r) cross. Replacing exprs. (4.10) and (4.11) into (4.12) we have that the filtered kernel reads (4.13) Taking r c ∼ 1, we finally have for scales larger than the horizon at the end of Inflation (4.14)

Charge excess noise kernel J 00 (x,ȳ)
To compute the rms value of the electric charge we define a kernel for the electric current in a similar way as we did in the previous section, namely the VEV of the symmetric two-point function J 00 (x, y) = 0| J 0 (x),J 0 (y) |0 (4.15) Replacing expr. (A.11) and (4.5) into (4.15) we obtain that the only terms which contribute to the kernel are To estimate the momentum integrals we again consider two regimes r 1 and r 1. For r 1 we obtain while for r 1 we get (4.18) We now proceed with the average over a scale λ using again a top-hat window function. Since r c is the coordinate where both expressions match we have Finally assuming r c ∼ 1, for large L D we obtain (4.20)

Primordial Weibel Instability
In this section we study the dispersion relation from eq. (2.45). We begin by building the functions that describe the background fluid, Z 0zz , q 10 and q 20 , and computing the parameterŝ τ and Υ 0 . From eqs. (3.32) and (4.14) we read and from the relation (3.30) we have which is consistent with the constraint Z * < 3Z 0 we assume before.
To estimate q 2 10 we use eq. The number density q 20 can be obtained from eq. (3.25) and reads The numerical value of the numerator in the r.h.s. of this equality must be read from Figure  2 according to the chosen value of ξ. The temperature in eqs. (5.3) and (5.4) can be fixed from the constraint (2.15) which becomes a polynomial equation for Υ 0 , namely The dimensionless temperature Υ 0 turns out to be a function of ξ and L D . This equation has only one real positive root shown in Figure 3 as a function of the averaging domain size L D for different values of ξ.
A simple comparison between the values of the temperatures from Figure 3 and the values of the masses in Figure 2 shows that Υ 0 µ is satisfied for a wide range of small masses that include Axion-Like-Particles (ALPs) [89][90][91][92]. In addition, since the scalar curvature of a Universe that expands as radiation-dominated one is vanishing, the fluid of scalar particles can indeed be considered as conformal [14]. The remaining parameter, the relaxation timeτ , is determined by the scattering cross section of the self-interaction between the fluid particles. Assuming that the short range electromagnetic interactions provide the main scattering mechanism, we get with α f s = e 2 /4π 1/137. This relaxation time sets the lifespan of the anisotropic pressure background, i.e. the time during which the instability can evolve. In Figure 4 we show the dimensionless relaxation timeτ as a function of the averaging domain size L D . Our one-fluid model of the instability is valid during the intermediate stage between the end of Inflation and the establishment of the bulk radiation content of the Universe that sets the end of the so-called Reheating era and the beginning of the radiation-dominated era [7,93]. Afterwards the charges will strongly couple to the bulk radiation and the relaxation time (5.6) will no longer be the main mechanism of relaxation or thermalization due to the inclusion of the interactions with radiation. Then the optimal scenario for the development of the instability is the one for which the relaxation time is of the order of or larger than the duration of Reheating.

Study of the Unstable Modes
Now that we have all the parameters of the model, we study the unstable solutions of the dispersion relation (2.45). In principle we have two sets of solutions, one corresponding to the s + modes and the other to the s − modes. We numerically solve the eq. (2.45) using the set of parameters defined at the beginning of this section.
We start analyzing the strength of the instability regarding both the magnetic field component (related to the ± modes) and the orientation of the propagation direction with respect to the anisotropic pressure background (related to the angles θ and φ). We find that s − (k) > s + (k) for any choice of ξ and L D regardless of the scale k and its orientation (θ, φ). In Figure 5 we show s + and s − as a function of k/k D , with k = κ/τ , k D = 2π/L D and L D 20, for three different values of gravity coupling as representative examples.
In view of these results we focus on the s − solution. The dependence of the instability exponent s − with respect to the orientation of the propagation direction is shown in Figure 6. In particular we show the values of s − as a function of θ with fixed ϕ = 0 for three different gravity coupling values and, in this example we also fix the scale k = k D and L D 60. Clearly the alignment (anti-alignment) of the propagation direction and the orientation of the anisotropic pressure tensor, i.e. θ = 0 (θ = π), produces the strongest instabilities. Hereafter we focus on the s − modes with θ = ϕ = 0.  . Normalized instability exponent s − as a function of θ regarding the orientation between the propagation direction and the anistropic pressure tensor, for different gravity coupling values. As a representative example we take L D = 60, k = 2π/L D and ϕ = 0. We observe that the largest values of the instability exponent s − arises for θ = 0 or θ = π, i.e. when the propagation direction is aligned to the orientation of the anisotropic pressure tensor. We also see a range of angles for which the perturbations are stable.

Conformal and unstable evolution of the magnetic field
To compute the evolution of the amplitude of the magnetic field we begin by rewriting the scale factor of the universe a(η) in terms of e-folds according to the conventions and definitions in [93]. We assume that Reheating begins at η = 0 and the scale factor, normalized to a(η = 0) = 1, evolves as in a radiation-dominated universe since the beginning of Reheating [7,93]. In consequence the scale factor reads a (η) = (1 + η) = e Nη , (5.7) with N η the number of e-folds at the time η. The evolution of the magnetic field for η <τ will be given by where B I is the magnetic fields at the end of Inflation. We consider Reheating scenarios whose duration η f satisfies η f ≤τ . In these cases the instability is operative during all the Reheating stage. From Figure 4 we may estimate that η f ≤ 5 × 10 4 which implies that the number of e-folds at the end of Reheating, i.e. N η f = N , must fulfill N 10.
Our aim is to compute the evolution of the magnetic field until the end of Reheating from the expression (5.8). In order to determine all the hydrodynamics parameters we first choose the averaging domain size L D such that the average is performed over all the scales that are inside the horizon at η f . Indeed for a given N we have L D = e N which coincides with the comoving horizon at η f implying that In Figure 7 we compare the instability exponent spectrum (σ − (k)η f ) and the conformal dilution exponent (2N ) in (5.8) as functions of the wavenumber scale k/k f with k f < k < 1 for different values of ξ and several Reheating scenarios determined by N . We observe that for N < 4, the magnetic field instability overcomes the conformal dilution for different ranges of scales k depending on the value of ξ. For shorter Reheating periods with larger values of ξ the instabilities can surpass the dilution, while for longer Reheating stages smaller values of ξ are needed to achieve a net amplification. Nonetheless there is limit N ∼ 4 for which the conformal dilution prevails for any value of ξ and k.

Bonus Track: Magnetic Helicity
The features of the evolution of magnetic fields are not only determined by the Maxwell and fluid equations, but also by the topological properties of the field. They are encoded in the magnetic helicity (MH), which is defined as the volume integral of the scalar product between the vector potential and the magnetic field [94], [95], [96]. It is interpreted as the signed number of twists due to individual field lines as well as due to the interlinkages of different lines.
Magnetic helicity is produced by the magnetic field creation mechanism itself, and it is not simple to envisage a primordial magnetogenesis process that creates fields with non-trivial topology. A possibility is the electroweak epoch, when violation of parity and charge+parity occurred [13,[97][98][99][100], and also the inflationary epoch [101,102] due to couplings of the electromagnetic field with other fields as the axion. In addition stochastic magnetogenesis mechanisms can also induce fields with non-vanishing rms values of MH [103,104]. The average of fluctuations, that determines the magneto-hydrodynamic parameters involved in σ − , is performed over a region of the size of the horizon at the end of Reheating, i.e. k D = 2π/L D = k f . We observe that for N 4 the unstable exponential growth compensates the conformal dilution or even produces a net amplification depending on the gravity coupling ξ. As N gets larger the instability can only overcomes the conformal dilution for lower values of ξ. For N 4 however the conformal dilution dominates for all possible ξ.
In this Sub-section we build the expressions to compute the amplification that an rms magnetic helicity would undergo due to the instability studied in this Section. Here we drop the explicit prime symbol in order to keep a minimal notation. We saw that we can build a basis with the three mutually orthonormal vectors {ε + (k),ε − (k),k}, where the first two are the eigenvectors associated to the eigenvalues λ + and λ − respectively and the third is the propagation direction. We assume that {ε + (k),ε − (k),k} is right-handed, which means that the vectors satisfyε + (k) ×ε − (k) =k. The eigenvectors also satisfyε ± (−k) =ε ∓ (k) under parity transformation of the propagation direction. To preserve the correspondence between eigenvectors and eigenvalues, the latter must transform as λ ± (−k) = λ ∓ (k) under parity.
Some time after the beginning of the evolution, when the decaying solutions have substantially faded away, each field of the vector sector can be written as a linear superposition of the unstable solutions according to η) is real. Since we analyze the dynamics only after some time of the evolution, we cannot consider an initial value problem. Instead we assume that where the angular brackets denote ensemble average, and that for not-so-long times, the solutions are still spherically symmetric. Therefore (5.12) Using the properties of the vectorsε ± and its relation with the transverse projector we read from (5.12) that

Magnetic Helicity
As briefly explained at the beginning of this subsection, the magnetic helicity is defined in a volume V as [94] H M (η) = V dV a (r, η) · b (r, η) (5. 16) where ∇ × a = b. As we are dealing with pure vector fields we can work in the Coulomb gauge and write the vector potential a (r, η) as In this way we obtain where we define the window function with V the volume over which the magnetic helicity is computed. The statistically spherical symmetry we assumed at the beginning of the evolution, or equivalently the lack of 'handedness' or preferred orientation, implies that the mean value of the magnetic helicity vanishes as we can check from (5.18). In consequence we directly compute its rms deviation at equal times. We assume that the variables are Gaussian. In view of the correlations given above, the terms that will survive are the ones with the same numbers of 'plus' and 'minus' signs. In addition the averages corresponding to the mean values vanish. Using the correlations (5.14)-(5.15), the parity properties of the eigenvectors and the eigenvalues, we write Observe that both s + and s − determine the evolution of the rms value of the magnetic helicity. To estimate the evolution of the rms value of H M we consider a large integration volume in (5.16) turning the window function (5.19) into a Dirac delta. By virtue of the orthogonality properties and the closure of the basis vectors, which shows that the Weibel instability also amplifies an initial rms value of the magnetic helicity. It is interesting to observe that while the magnetic field grows mainly with the largest eigenvalue σ − (k), the growth exponent of the magnetic helicity is determined by the sum of both roots σ − (k) + σ + (k). Indeed the rms value of the magnetic field is given by where the last expression stems from the fact that σ + (k) < σ − (k).
We might define a coherence scale for the helical fields as the ratio of the rms values of the magnetic helicity to B 2 , which according to eqs. (5.22) and (5.24) will evolve as As σ + (k) < σ − (k) for all modes (see eq. (2.41) and Fig. 5) it is always satisfied that σ + (k) + σ − (k) < 2σ − (k) so the denominator will grow faster than the numerator. This means that the proposed coherence length for the helical field will shrink as the instability evolves.

Conclusions
We study the feasibility of the development of a relativistic Weibel instability during the Reheating stage in the very early universe prior to the establishment of the main bulk radiation. We particularly use the second order magneto-hydrodynamic description developed in Ref. [74]. We assume that during Inflation there exists an spectator charged scalar field arbitrarily coupled to gravity in its de Sitter vacuum state. Gravitational particle creation occurs at the transition from Inflation to the following stage and the field state becomes an excited many particles one which admits a hydrodynamic description. Quantum fluctuations of the scalar field produce rms deviations of T µν and J µ . We average these quantities over super-horizon regions at the end of Inflation whose comoving size is L D . The quantum rms values are translated into non-equilibrium fluid parameters through the simplest matching procedure between quantum and statistical fluctuations. In particular the energy density, related to the VEV of T 00 , the anisotropic pressures, related to the rms VEV of the energy-momentum tensor, and the charge excess, related to the rms VEV of J 0 , are the relevant ones for the unfolding of the vector instabilities.
We find that for Reheating scenarios with N 4 e-folds, which correspond to lapses of time shorter than the relaxation time τ of the fluid, the unstable exponential growth of the magnetic field compensates the conformal dilution or even produces a net amplification for sub-horizon scales at the end of Reheating depending on the values of the gravity coupling (see Figure 7). Moreover, the instability also amplifies the magnetic helicity, while concentrating it at small scales, thus improving the chances for an helical field to survive.
Although the instability produced by the inflationary initial quantum fluctuations is efficient only at small scales compared with the scales of astrophysical interest, it might be relevant to produce strong helical magnetic fields that work as seeds for posterior amplification processes during Reheating itself or the radiation-dominated stage.
Finally we remark the relevance of the Second Order magneto-hydrodynamics to describe the dynamics of non-ideal processes during the Very Early Universe, in this case by capturing the anisotropic pressure tensor produced by the inflationary quantum fluctuations in a consistent magneto-hydrodynamic framework. support through project UESC/PROPP N • 00220.1300.1876.

A Charged Scalar Field
In this section we give a nut-shell review of charged scalar fields in curved space-time and compute the main functions that to characterize the quantum gas as a real relativistic plasma: the VEV of T 00 Q , the Hadamard two-point function of the traceless part of T ij Q and of the charge density J 0 Q . A charged scalar field is mathematically described by the complex fields Φ(x µ ) and Φ † (x µ ), which can be written in terms of real fields as where Φ 1 and Φ 2 commute and satisfy the Klein-Gordon equation for a massive field arbitrarily coupled to gravity Φ − µ 2 + ξR Φ = 0 .

(A.3)
We assume that all fields and parameters are already dimensionless. The Energy-Momentum tensor (EMT) for each of the real scalar fields is [87,88] (A.4) The symbol represents the relativistic D'Alambertian, R µν the Ricci tensor and R the scalar curvature. The dimensionless mass of the field is µ = m/H and ξ is the coupling to the curvature. The values ξ = 0 and ξ = 1/6 correspond to the minimal and conformal coupling respectively. The operator applied to Φ 2 can be written using the Klein-Gordon equation as Φ 2 = 2∂ α Φ∂ α Φ + 2µ 2 Φ 2 + 2ξRΦ 2 . (A.5) The EMT for each real field is then given by eq. (4.1), with the trace Rescaling the scalar field as Φ = Φ a , (A. 8) expr. (4.1) becomes T µν = a −2T µν withT µν given by eq. (4. 3) The electric four current of a charged scalar field is defined as [86] which in terms of the real fields is After replacing (A.8) for each field we have that J µ = a −2J µ withJ µ given in eq. (4.5)