The charged inflaton and its gauge fields: preheating and initial conditions for reheating

We calculate particle production during inflation and in the early stages of reheating after inflation in models with a charged scalar field coupled to Abelian and non-Abelian gauge fields. A detailed analysis of the power spectra of primordial electric fields, magnetic fields and charge fluctuations at the end of inflation and preheating is provided. We carefully account for the Gauss constraints during inflation and preheating, and clarify the role of the longitudinal components of the electric field. We calculate the timescale for the back-reaction of the produced gauge fields on the inflaton condensate, marking the onset of non-linear evolution of the fields. We provide a prescription for initial conditions for lattice simulations necessary to capture the subsequent nonlinear dynamics. On the observational side, we find that the primordial magnetic fields generated are too small to explain the origin of magnetic fields on galactic scales and the charge fluctuations are well within observational bounds for the models considered in this paper.

In this paper we re-visit particle production in locally gauge invariant models with Abelian and non-Abelian fields coupled to charged scalar fields. In these models we assume that a component of the charged scalar field plays the role of the inflaton condensate. Care has to be taken with gauge fields because they have to satisfy certain constraint equations (along with the usual evolution equations). The natural gauge redundancy can lead to complications in quantization or to spurious gauge modes during numerical evolution. The non-zero vacuum expectation value (vev) of the inflaton condesate during inflation and reheating changes some of the common results for gauge fields coupled to scalar fields with zero vevs. Moreover, certain common gauge choices become ill-defined when the inflaton condensate starts oscillating at the end of inflation. Finally, if significant particle production occurs, back-reaction becomes important and classical lattice simulations are needed to fully explore the nonlinear dynamics of the fields. Initial conditions for such lattice simulations can be nontrivial because of the constraints on the different variables that must be satisfied. We pay special attention to all of these issues in this work. While we restrict ourselves to "minimal" models consistent with local gauge invariance, our techniques and results should carry over to more complicated scenarios such as [85,86].
We calculate particle production during inflation using gauge invariant variables and with proper accounting for the constraints (Section 3). The use of gauge invariant variables naturally avoids issues with spurious gauge degrees of freedom and makes the quantization and subsequent evolution of perturbations particularly transparent. We provide power spectra for the electric and magnetic fields at the end of inflation, along with simple analytic estimate for their shape. In Appendix A we explain their shape via approximate analytic calculations.
While useful during inflation, gauge invariant variables become ill-defined when the inflaton starts oscillating. We argue for the use of well defined Coulumb gauge variables for analysing nonperturbative particle production during preheating at the end of inflation. In Section 4, we carry out a Floquet like analysis for the resonant production of the gauge fields. Here, we point out some minor discrepancies in the literature regarding the productions of the gauge invariant longitudinal component of the electric field. We then estimate the end of preheating by calculating the time when back-reaction of the resonantly produced gauge fields becomes important (Section 4.2). In Appendix B we provide technical details for quantizing and calculating the back-reaction in the Coulomb gauge.
Once nonlinear effects become important, simulations become essential. Nonlinear simulations with gauge fields (especially non-Abelian fields) can be challenging because of the large number of components and the necessity of satisfying constraint equations. In Section 5 of this work we provide a simple prescription for setting up initial conditions for such lattice simulations which can be applied in any gauge. In our prescription, the lattice initial conditions naturally satisfy the necessary gauge constraints. We note that our initial conditions accurately account for metric perturbations, field interactions and gauge constraints to linear order. We provide an example of lattice initial conditions in temporal gauge which is a common choice for simulations. We arrive at these initial conditions via gauge invariant variables; this serves as a model to set up initial conditions in any gauge. We will carry out actual lattice simulations in future work.
We have tried to make the paper self-contained, so that it can be used for future reference easily. To this end, we provide the necessary equations for the perturbations of metric and matter fields (scalar and gauge fields) in position and Fourier space for gauge invariant variables. While we work with gauge invariant variables as far as possible, we also provide a dictionary to translate our results to other popular gauges.
For pedagogical purposes, we carry out the analysis for Abelian fields first (Section 2). We show in the later half of the paper (Section 6) how the non-Abelian analysis can be reduced to an analysis of multiple copies of the Abelian case in the linear regime. Hence, the analysis for the Abelian case, including the setting up of the lattice initial conditions, can be easily carried over to the non-Abelian case. We show how to apply the developed techniques to a SU (2) model and its extension: a SU (2) × U (1) model.
We discuss the observational consequences of a charged inflaton and its gauge fields in Section 7. We summarize our results in the Conclusions section and discuss future directions.

The Abelian model
We consider an action with matter minimally coupled to gravity where R is the Ricci scalar, g is the determinant of the metric and m Pl is the reduced Planck mass. The matter action contains a complex scalar field ϕ and the gauge field A µ : where the field tensor F µν for the gauge fields and the gauge-covariant derivative D µ ϕ are given by In the above equations, ∇ µ is the usual Levi-Civita connection. The action S m is invariant under local U (1) gauge transformations where α(x ν ) is an arbitrary real function of space and time. The total action is also invariant under space-time differomorphisms. The gauge symmetry implies that not all of the components of the 4-vector A µ and the real and imaginary parts of the scalar field are physical degrees of freedom (dof). We remedy this redundancy by working in the appropriate set of gauge invariant variables or by fixing the gauge. The redundancy due to space-time diffeomorphisms is handled in a similar fashion.
We will present our answers as power spectra of the electric and magnetic fields, which are defined in the usual way [95]: 1 (2.5)

U (1) gauge invariants
When the field ϕ = 0, it can be written in polar co-ordinates as Under the local U (1) gauge transformation (see eq. (2.4)) ρ → ρ and Ω → Ω − α. It is convenient to work in local U (1) gauge invariant variables given by the following five fields: In these variables the matter action becomes where V (ρ) = V(|ϕ|). It is worth noting that the variable Ω appearing in eq. (2.6) does not make an appearance in the above action. There is no need to worry about the U (1) gauge redundancy when working with gauge invariant variables. Note that the E and B fields are invariant with respect to U (1) transformations and their expressions in terms of G µ are identical to those in terms of A µ :

Diffeomorphism invariants
We will work in a perturbed Friedmann-Roberston-Walker space-time with the metric:

10)
1 ijk is the antisymmetric Levi-Civita tensor. It does not raise or lower (spatial) indices.
where φ(x σ ), B(x σ ), ψ(x σ ), E(x σ ) are scalar perturbations, S i (x σ ), K i (x σ ) are divergence-free 3-vector perturbations, andh ij (x σ ) is a traceless transverse 3-tensor perturbation. In this perturbed space-time, we define the perturbations of the following U (1) invariant variables: where G 0 (x σ ) and G (x σ ) are scalars and G ⊥ i (x σ ) is a divergence-free 3-vector. Note that the gauge fields vanish at the background level: the spatial components are zero from isotropy of the Friedmann-Robertson-Walker (FRW) background and the equations of motion will set the background temporal component to zero. Hence we will work at the linear level in G µ .
For our scenario, there are two physical scalar metric perturbations, with scalar field perturbation and scalar parts of the gauge field adding three more. We choose to work with the following five diffeomorphism invariant combinations: where H = ∂ τ ln a. The first two are the standard Bardeen variables. G 0 and G are diffeomorphism invariant since the gauge field vanishes at the background level. Similarly, for the vectors we chose to work with the following diffeomorphism invariant com- whereṼ and G ⊥ are divergence free. The traceless, transverse 3-tensor perturbationh ij is already diffeomorphism invariant. Similarly, the electric and magnetic fields defined in eq. (2.9) are already diffeomorphism invariant. It is convenient to split the electric and magnetic fields into divergence-free and curl-free parts (2.14) Note that B = 0 from the definition of B i in eq. (2.9).

Equations of motion
The general equations of motion for the matter and metric fields in curved space-time take the following form: where G µν is the Einstein tensor and the energy momentum tensor is given by In terms of the U (1) gauge invariant variables defined in Section 2.1, the above equations become Equation (2.18) implies the following definition of the conserved 4-current: The equivalent of Maxwell's equations for the electric and magnetic fields are Next, we write down the equations of motion for the background (space independent) fields and linearized perturbations around these background fields in terms of gauge invariant variables.

Background
Assuming the scalar field plays the role of the inflaton and treating the gauge fields as perturbations, the evolution ofρ(τ ) can be determined from eq. (2.17): where H is given by the 00 background Einstein equation The electric and magnetic fields vanish at the background level.

Linearized perturbations in position space
From eq. (2.17) and eq. (2.18) we get the equations of motion for diffeomorphism and U (1) gauge invariant scalar perturbations: where eq. (2.26) is the linearized version of the Gauss constraint. Note that the scalar components of the gauge fields G 0 and G do not depend on the metric perturbations. The evolution of the U (1) gauge and diffeomorphism invariant vector perturbations involving matter fields can be obtained from The perturbations G ⊥ i do not couple to metric perturbations. As mentioned earlier, the tensor perturbations are also decoupled from the matter fields at the linear level. We shall not consider tensor perturbations any further in this paper.
We now turn to the Einstein equations. The energy-momentum tensor in eq. (2.19) is quadratic in the G µ (with G µ = 0 at the background level). Hence, to linear order in the perturbations, the energy-momentum tensor depends only on the perturbations in ρ and the metric. Also T i j = 0 for i = j; there is no anisotropic stress. The linearised Einstein equations (for scalar perturbations) yield Φ = Ψ , (2.28) Note that the gauge fields are completely decoupled from the scalar metric perturbations.
The picture for vector perturbations is even simpler -the linearized Einstein equations for the vector perturbations involve only metric perturbations: which do not affect the matter vector perturbations. At the linearized level, the equations of motion for the electric and magnetic fields defined in eq. (2.9) are simple 2 On the right-hand sides of the last two equations, we have defined the scalar perturbations and divergence-free vector perturbations in the 3-current density respectively.

Linearized perturbations in Fourier space
For calculational purposes, we move to Fourier space. Fourier space is also particularly convenient for solving the various constraint equations, which essentially become algebraic in the relevant variables. Our Fourier convention is f (x , τ ) = f k (τ )e ik ·x d 3 k . We begin with the equations of motion for the scalar perturbations. From the (Fourier space version of the) constraints, eqns. (2.28), we can substitute the gravitational potential Ψ k and its derivative into the evolution equation for δρ k (cf. eq. (2.24)) to obtain an equation of motion which only involves δρ k : (2.31) The remaining scalar perturbations, G 0 and G , are governed again by an evolution equation, eq. (2.25), and a constraint, eq. (2.26). Before moving to the Fourier transformed versions of these equations we define the longitudinal (i.e. curl-free) component of the space-like part of G µ : The Fourier transform of G L can be expressed in terms of a longitudinal polarisation vector, L k , as follows: where we shall call the scalar G L k , the longitudinal mode. The polarisation vector has the following properties: The Fourier transformed equations, eqns. (2.25) and (2.26), then take the form There is a similarity between the pairs G L k , G 0k and (δρ k , Ψ k ). G L k and δρ k are both dynamical fields, evolved according to second order in time equations of motion, eq. (2.36) and eq. (2.24), respectively. Each of these perturbations has its own auxiliary field, G 0k and Ψ k respectively, determined by a constraint equation. Substituting the auxiliary field (i.e. G 0k ) into the equation of motion we obtain the following expression in terms of G L k only: We now turn to vector perturbations. The matter vector perturbations are decoupled from the metric vector perturbations. Similarly to the longitudinal case, we introduce a pair of transverse polarisation vectors T ± k which satisfy The divergence-free perturbations in terms of these polarization vectors are with the equations of motion The fact that Hubble friction does not appear in the evolution equations for the transverse modes is because of conformal invariance of massless gauge fields. However, the longitudinal components (which exist when the gauge field is effectively massive) do feel Hubble friction. One can also rewrite the electric and magnetic fields, and the 3-current current, in terms of longitudinal and transverse modes, e.g.
which in terms of the polarization vectors in Fourier space becomes Similar expansions hold for j i and B i with the exception B L k = 0. Below we give the expressions used in the subsequent sections to calculate the primordial power spectra of the longitudinal and transverse modes of E and B (2.43) The charge and current densities can also be expressed in terms of the G field's longitudinal and transverse modes (2.44) The first identity is simply −kE L k = a 2 j 0k , i.e. the Gauss constraint. In a consistent quantum analysis of the perturbations during and after inflation, one cannot set j 0k and j L k to zero by hand (this differs from the treatment in [62], [96]). This is because the vacuum fluctuations in G L k can be enhanced due to horizon crossing during inflation or non-adiabatic particle production during preheating. We shall see this aspect in detail in the subsequent sections.

Gauge transformations
In the upcoming sections we will use either the gauge invariant variables discussed above or work in some particular gauge depending on which approach is best for the problem at hand. In this short section we provide the relationships between variables in some of the popular gauges used in the literature and the gauge invariant ones. These relationships can also be used to move from one gauge to another. When the variables are well defined, we can use the equations of motion and the solutions in the gauge invariant case to recover the corresponding expressions in our gauge of choice. Occasionally, some of the variables become ill-defined or the relationships between variables require a patching up of co-ordinate maps (for example during reheating). Such cases can be dealt with on an individual basis, or one simply derives the equations of motion and solutions directly from the equations of motion themselves.
In different gauges (but not the gauge invariant case) the complex scalar field is represented as ϕ = (ϕ 0 + iϕ 1 )/ √ 2, and the gauge fields by A µ . For perturbations, we will always work around an FRW background. Using the global U (1) invariance, we set the homogeneous, imaginary part of ϕ,φ 1 = 0. That is ϕ = (φ 0 + δϕ 0 + iδϕ 1 )/ √ 2. For the homogeneous partφ 0 = ±ρ with the ± sign accounting for the change in sign during an oscillation through zero. In Fourier space, the transverse modes and the scalar perturbations along the direction of motion of the homogeneous field (in all the gauges discussed below) are related to the gauge invariant variables as follows: with the equations of motion for δϕ 0 k and A T ± k being identical to the gauge invariant case with ρ →φ 0 .
Coulomb gauge: In this gauge ∂ i A i = 0. In Fourier space we get Unitary gauge: In this gauge ϕ 1 = 0, which yields The equations in the Unitary gauge are identical to those in the gauge invariant one.
Temporal gauge: In this gauge A 0 = 0. However, the theory is still invariant under the time-independent transformation Hence, we completely fix the gauge by choosing an α such that at some moment of time, τ = τ in , δϕ 1 k (τ in ) = 0. With this condition, we have Lorenz gauge: In this gauge ∇ µ A µ = 0. In this case where G(τ, η) is the Greens function of the linear operator L τ ≡ ∂ 2 τ +2H(τ )∂ τ +k 2 . Arriving at the above form of the relationship between variables requires a bit of explanation. In Fourier space the Lorenz gauge condition translates to ∂ τ A 0k + 2HA 0k − kA L k = 0, and the equation governing . This equation yields the particular solution above only if we can set the complementary solution to zero. This can always be done since there is a residual degree of freedom χ such that under the transformations A µ → A µ + ∇ µ χ and δϕ 1 → δϕ 1 −φ 0 g A χ/2, the theory remains invariant; provided χ obeys ∇ µ ∇ µ χ = 0. In Fourier space we get L τ χ k = 0. Since the operator evolving χ k is identical to the one evolving A L k , we can always choose χ k such that the complementary part of A L k vanishes, thus arriving at the particular solution provided above.

Inflationary dynamics
The background dynamics are relatively straightforward during inflation. At the phenomenological level, with an appropriate choice of the potential V and initial conditions we can easily arrange for −∂ t H/H 2 = −a∂ τ (H/a)/H 2 1 for sufficient number of e-folds. For simple models, this corresponds to H = (H/a) ≈ const andρ ≈ const during inflation. Inflation ends when ∂ τ H = 0, when accelerated expansion stops and the field starts rolling quickly. Assuming such a background solution has been found, we focus on the quantum fluctuations around this classical background. Quantization of constrained systems, like the problem at hand, can be tricky. We find that by working in Fourier space with gauge invariant variables and substituting the constraints before quantizing, the process becomes straightforward. Once the appropriate quantized solutions for the scalar and gauge fields are available, we construct the power spectra of the electric and magnetic fields at the end of inflation.

Quantized scalar and vector perturbations
For the purposes of quantization, it is convenient to write down the action for the Fourier components of the dynamical perturbation variables left after imposing the constraints. The equations of motion for these variables δρ k , G L k and G T ± k were provided in the previous section (see eqns. (2.31), (2.37) and (2.40)). The total quadratic action of these variables naturally splits into four parts: where S I are the quadratic actions for the perturbations in the I-th variable f I with The explicit forms of b I (k, τ ) and ω I (k, τ ) will be provided below for the different components. We first outline the general quantization procedure common to all of the components. The conjugate momentum density of the I-th variable where we have made use of f I * k = f I −k . The field operatorsf I k along with their conjugate momenta operators π I k must obey the standard equal time commutators: Note that our fields are not canonically normalized. The non-vanishing commutator and eq. (3. 3) The field operatorsf I k can be expanded in terms of operatorsâ I k and mode functions u I k (τ ) aŝ where each of the mode functions u I k (τ ) satisfies the corresponding field equations of motion obtained by varying the action S I (same equations as those satisfied by f I k ) The final ingredient needed for evolving the mode function u I k (τ ) (and hence the field operators), are the initial conditions for the mode functions. Their normalization will in turn also determine the commutation relation for the operatorsâ I k . We will determine the initial conditions by constructing WKB solutions for the mode functions satisfying eq. (3.7) at very early times.
To proceed to the WKB solutions for the initial conditions we need the explicit forms of b I (k, τ ) and ω I (k, τ ), which are provided below. (3.9) One can check that extremising the action in eq. (3.1) with respect to each of the field perturbations gives the corresponding equations of motion from the previous section, namely eqns. (2.31), (2.37) and (2.40). At early enough times during inflation as a → 0, a given k mode of interest will be deep inside the horizon k H and will dominate all other physical scales (for example, k (ρg A a)) as a → 0. At such early enough times during inflation ∂ 2 τ ln b I , (∂ τ ln b I ) 2 ω 2 I . This hierarchy can be verified by noting that for each component ω 2 With this information at hand, the WKB solution of eq. (3.7) at early enough times during inflation is 3 where for the WKB solution to be valid The time independent normalization of the mode functions ( (2π) 3/2 √ 2 −1 ) was chosen so that the usual commutation relation for the creation and annihilation operators is satisfied. That is, using the above mode functions in eqns. (3.5) and (3.6) implies that at early times Since these operators are time-independent, these relationships remain true at all times. Thus by starting with the initial conditions for the mode functions u I k in (3.10), we can evolve them using eq. (3.7) to any later time and obtain the necessary power spectra. The forward time evolution will take the initially sub-horizon and/or ultra-relativistic (k H and/or k ρg A a) solutions to superhorizon and non-relativistic ones.
Before ending this section we make some general comments about the early time solutions for the mode functions u I k during inflation written explicitly below: (3.13) The early time solution for the mode function u ρ k of δρ k reflects the fact that at early enough times, we are simply dealing with an effectively massless (k 2 a 2 ∂ 2 ρ V ) field on subhorizon k H scales where metric perturbations are negligible. The early time solution is identical to the mode functions for the Minkowski vacuum apart from the trivial a(τ ) scaling. The mode function u L k of the longitudinal component of the gauge fieldsĜ L k is related to u ρ k in a simple way: The rescaled field (ρg A /2k)Ĝ L k behaves as a massless scalar field, which is a manifestation of the Goldstone Boson Equivalence Theorem. The scalefactor does not appear in the early time mode functions u T ± k of the transverse gauge field modesĜ T ± k . This reflects the fact that massless transverse modes of gauge fields are conformally invariant.
Finally, we note that the above early solutions were constructed assuming slow roll inflation. However, for large enough k these solutions remain valid even at the end of inflation as is to be expected, albeit with slightly different conditions on k. The conditions for these solutions to be valid are k 2 H 2 , a 2 ∂ 2 ρ V for u ρ k ; k ∂ τ (aρ)/(aρ), aρg A /2 and k 2 |∂ 2 τ (aρ)/(aρ)| for u L k ; and k aρg A /2 for u T ± k . These are useful for checking the ultra-relativistic solutions at the end of inflation.

Inflationary power spectra
Let us now use the developed formalism to compute the power spectra of the matter fields at the end of inflation -the first moment when ∂ τ H = 0. 4 We will calculate the power spectra of the gauge invariant electric and magnetic fields.
The correlation functions of the electric and magnetic fields can be written in terms of the longitudinal and transverse mode functions as follows: where from (2.43) The power spectra of the electric and magnetic fields are defined as 5 For concreteness, we compute these power spectra at the end of inflation for the chaotic inflation scenario with V (ρ) = (1/2)m 2 ρ 2 . However, our results hold more generally. For an analytic understanding of the features in the power spectra we solve the equations of motion in de Sitter space to zeroth order in slow-roll (i.e. we ignore time derivatives of the inflaton) in the Appendix A. We assume that inflation lasts 60 e-folds and take the inflaton mass to be m = 10 −6 m Pl . The longitudinal and transverse fields,Ĝ L k andĜ T ± k are evolved numerically from very early times when they are deep inside the horizon with k H, k aρg A /2 and k 2 ∂ 2 τ (aρ)/aρ) , and have corresponding electric and magnetic WKB power spectra (cf. eqns. (3.13), (3.15), (3.16)) (3.17) 4 Which is equivalent to the standard expressionä(t) = 0, where t is cosmic time. 5 The two transverse modes would have different power spectra if there was axion coupling which we do not consider here. The way we have split the transverse modes into states ± with well-defined helisities makes our analysis easily extendible to the case of charged Higgs with an additional axion-like interaction.
The power spectra of the electric and magnetic fields at the end of inflation are shown in Fig. 1.
To understand these plots, a 'Compton wavenumber' k C corresponding to the effective mass is particularly important: where the time dependent quantities on the right-hand side are all evaluated at the time of interest (usually at the end of inflation). The spectra behave differently based on the relative size of k C and the Hubble scale H. When the Compton wavenumber ofĜ L k andĜ T ± k is subhorizon, i.e. k C H, we have As is evident from the above scalings, the magnetic field and both the transverse and longitudinal electric field power spectra have double power-law forms, with k C setting the break in all three of them. This is indeed expected for ∆ 2 cares only about k C (cf. eq. (2.40)). However,Ĝ L k is affected by the expansion rate as well. From its equation of motion we would expect to see something in ∆ 2 E L near k C and the Hubble scale, H. The reason why there are no features in ∆ 2 E L near H in the strong coupling regime is given in Appendix A. On the other hand when k C H, we find that (3.20) In this regime, the transverse modes are still unaffected upon crossing the Hubble horizon. ∆ 2 is again a double power-law. However, this time the break is at k = k C (k C /H). The additional suppression of k C /H is unexpected, at least if one looks at the equation of motion ofĜ T ± k , eq. (2.40). For the interested reader, we explain this in Appendix A. Also note that in this regime ∆ 2 B T ± is a single power-law, given by deep subhorizon Minkowskian power spectrum. For 0.1 < k C /H < 1 there is a crossover behaviour from a double power-law of the k C H regime to the single power-law of k C H. The power spectrum for the longitudinal component of the electric field ∆ 2 E L exhibits a power excess on super-Hubble scales, with the correction factor found numerically to be 1 < T k < 2 at the end of inflation. On super-Compton scales, k < k C , T k ≈ 2. We also do not see any features in the range k C < k < H in ∆ 2 E L because the k-dependent prefactor multiplying ∂ τĜ L k in the definition ofÊ L k , eq. (2.43), cancels the additional k-dependence.   Figure 1. The power spectra of the transverse and longitudinal electric and magnetic fields at the end of inflation for m 2 |ϕ| 2 inflation (solid lines) for different strengths of the couplings g A between gauge fields and the inflaton condensate. Power spectra at the end of inflation ignoring the effects of expansion (dashed lines) are shown for comparison. We show three different spectra: the longitudinal electric field, ∆ 2 E L (green), the transverse electric field, ∆ 2 E T ± (black), and the magnetic field, ∆ 2 B T ± (orange). The mode functions of the fields responsible for these spectra can have two distinct forms depending on the magnitude of the scale k C ≡ aρg A /2 whereρ is the magnitude of the inflaton field at the end of inflation. For the upper row with k C H, ∆ 2 E L (green) and ∆ 2 B T ± (orange) are nearly single power-laws, whilst ∆ 2 E T ± (black) is a double power-law with a break k ∼ k C (k C /H). For the lower row with k C H, all three power spectra have double power-law forms, with a break at k ∼ k C . Note that one can also use these plots to read-off the power spectrum of the charge density, ∆ 2 j0 = (k/a) 2 ∆ 2 E L . More detailed expressions for the spectra are provided in eqns. (3.19) and (3.20). An analysis of these spectra in de Sitter space (which explains the different power-laws) is provided in Appendix A. Note that the Hubble parameter H (and its conformal counterpart H) at the end of inflation was used to make the relevant quantities dimensionless on the horizontal and vertical axes.
Note that T k is time dependent and T k → 1 in de Sitter space-time (see Appendix A for the evaluation of the power spectra in de Sitter space-time). The deviation of T k from 1 is observed towards the end of inflation. It is thus important to solve for the mode functions in a background that deviates from de Sitter in order to capture this effect and obtain the correct power spectra at the end of inflation.
We also give the power spectra calculated under the assumption that H = 0 andρ = const (hence k C = const), cf. dashed lines in Fig. 1. These are 'quasi-Minkowski' conditions in the sense that the effects of expansion are ignored. We can then calculate the evolution of the mode functions and obtain the following expressions for the power for electric and magnetic fields: When k C H expansion is important for superhorizon modes k < H. However, if k C ≥ H the effects due to expansion are negligible and eq. (3.21) is a good approximation to the power spectrum in the electric and magnetic fields.

Preheating dynamics
At the end of inflation the inflaton field begins to oscillate around the minimum of its potential. The linearised equations describing the evolution of the fluctuations in the inflaton and gauge fields will then contain oscillating terms. Such oscillating terms can lead to an exponential growth of matter fluctuations. This preheating period begins soon after the end of inflation and its end is marked by the back-reaction of the fluctuations on the inflaton condensate and/or the moment when the linearised analysis of the fluctuations stops being a good approximation [9].
The power spectra of electric and magnetic field at the end of the preheating era are shown in Fig. 2. These spectra are calculated numerically (for quadratic inflation) and properly account for the quantum nature of the fields, effects of expansion and metric perturbations. A detailed understanding of these spectra is the main goal of this section. As we will see, calculating these spectra using gauge invariant variables (and Unitary gauge variables) becomes somewhat unwieldy (in particular for the longitudinal components of the gauge fields). The Coulomb gauge turns out to be well suited for this calculation.
In Section 4.1, we will use Floquet theory to explore the instabilities in matter perturbations using gauge invariant variables and gauge dependent variables. In Section 4.2, a Hartree approximation is used to understand the effects of back-reaction on the homogeneous condensate and quantify a time when preheating ends.

Floquet analysis
We ignore metric perturbations for the Floquet analysis of the instabilities in the matter perturbations. This is a plausible approximation since the vector and tensor metric perturbations are decoupled from matter anyway, while the scalar metric perturbations are suppressed on subhorizon scales. We also ignore expansion at the background level since the universe does not expand much during the short period of preheating. For notational simplicity we drop quantum operators from our expressions while carrying out Floquet analysis. 6 The linearized equations of  Figure 2. The power spectra of the transverse and longitudinal electric and magnetic fields after a few inflaton oscillations at the end of inflation: the longitudinal electric field, ∆ 2 E L (green), the transverse electric field, ∆ 2 E T ± (black), and the magnetic field, ∆ 2 B T ± (orange). Each plot was evaluated for a different coupling g A , and the result is presented at a fixed time t br ≈ 10 2 m −1 after inflation. t br is the time of back-reaction for the specific coupling g A m Pl /m = 10 3 (lower right corner). The effects of non-adiabatic particle production during preheating is clearly visible in the lower row (large coupling) for k/H m/H k C /H ∼ k C /H. The boundary between small and large coupling g A m Pl /m ≈ 1, can be understood via a Floquet analysis of the instabilities. A requirement of broad resonance [9], leads to the upper bound on k. The features in the top row (small couplings) can be also understood in terms of the inflaton oscillations. We have rescaled ∆ 2 E T ± and ∆ 2 B T ± by a 4 , and rescaled ∆ 2 E L by a 5 (with a = 1 at the end of inflation) to roughly separate the effects of resonant particle production from the expected red-shifting of the fields. We also plot the power spectra before resonant particle production begins, i.e. right at the end of inflation (dashed lines) for comparison. For ease of comparison with Fig. 1, H (and its conformal counterpart H) at the end of inflation was used for constructing the relevant dimensionless ratios on the horizontal and vertical axes. Note that H br ≈ 0.25H, that is k/H br ≈ 4k/H, so the main features in the lower row are all significantly subhorizon at the time of evaluation. motion for the Fourier modes of the matter fields have the general form (cf. Section 3.1, but now ignoring expansion, quantum operators and metric perturbations) where ω I (k, τ ) and b I (k, τ ) are periodic because of the oscillating inflaton. Floquet theory tells us that eq. (4.1) has solutions of the form where P I k ± (τ ) are periodic functions having the same period T as the inflaton, and µ I k is the Floquet exponent. The Floquet exponents and the periodic functions can be calculated by solving the equation of motion for f I k twice in the interval τ 0 ≤ τ ≤ τ 0 + T with f , ∂ τ f I(2) k = (0, 1) as initial conditions. The real part of the Floquet exponents is given by (cf. e.g. [12] for details) where f and ∂ τ f are evaluated at τ = τ 0 + T . The effects of the more general initial conditions (for example, those from the end of inflation) can then be incorporated by appropriately scaling the periodic solutions.

Gauge invariant analysis
For the gauge invariant variables, the specific forms of the coefficients b I and ω I in eq. (4.2) are quite simple in flat spacetime. For the inflaton perturbations f ρ k = δρ k , (cf. eq. (3.8), with a = 1, Similarly, for the transverse components of the gauge fields (4.5) In both of these cases, one can calculate the respective Floquet exponents µ ρ k and µ T ± k onceρ(τ ) is obtained from the background dynamics using the usual techniques. 7 For a chaotic inflation potential V (ρ) = (1/2)m 2 ρ 2 , [µ ρ k ] = 0 whereas [µ T k ] is shown in Fig. 3.
Longitudinal mode: The analysis of the longitudinal mode of the gauge field f L k = G L k is a bit more subtle. For this case (4.6) 7 From the background equation of motion forρ near the minimum of V (ρ), eq. (2.22), one might infer thatρ is oscillatory with positive and negative values. This seems apparently at odds with the positive definiteness of ρ as implied by its definition as a "radial" variable. This is of course a minor inconvenience due to our choice of "radial" gauge invariant variables and can be understood in terms of a discontinuous jump of the angular variable at the origin which we have ignored. There is the other issue of the definition Gµ whenρ = 0, cf. eq. (2.7), though it remains well defined away from this point.
To deal with the singularity,ρ = 0, in the coefficient ∂ τ ln b L of ∂ τ G L k (the "damping" term) appearing in eq. (4.1), we make the following change of variables (4.7) Thereby, we lose the singular damping term in eq. (4.1) for f L k =Ǧ L k , at the expense of the singularity inȟ(ρ(τ )). The explict form of the coefficients after the transformation arě . (4.8) We note that the second term inω 2 L , proportional to ∂ 2 τρ (τ )/ρ(τ ), is never singular provided ρ −1 ∂ρV (ρ) is well-defined at the originρ = 0 8,9 , which is usually the case. Hence the solutionš are well behaved with no singularities present in the Floquet exponents or in the periodic functions. Reverting back to the original variables, cf. eq. (4.7), we have (4.10) The Floquet exponents, µ L k , remain unaffected and singularity-free and are shown in Fig. 3. Only the periodic functions, P L k ± (τ ), are singular (or have 'spikes') at the zeroes ofρ, because of theȟ(ρ(τ )) factor. Importantly, the physical observables: the longitudinal electric field E L k and the charge and current densities, j 0k and j L k , are immune to the singularities because of thē ρ-dependent pre-factors in their definitions cancel appropriately (cf. eqns. (2.43) and (2.44)).
We note that an identical analysis holds for the Unitary gauge ( [ϕ] = 0) as well. In that gauge, the equation of motion for the A L k has b L and ω L given in eq. (4.6) withρ → [φ] (cf. Section 2.4). Hence A L k grows exponentially as well with periodic spikes occurring whenever [φ] = 0. This still leads to spike-free exponential growth of the longitudinal electric field and charge and current densities. We also confirm this result below by carrying out the calculation in the Coulumb gauge where no singularities are present in the equations of motion. Our conclusion is different from that of [87], where the authors argued that the coefficient, We also note that the authors in [97] argued that there should be strong resonance in the longitudinal modes (stronger than the transverse modes), though they did not pursue this issue in detail. We find that the resonance in both longitudinal and transverse modes is comparable. 8 Recall that ∂ 2 τρ /ρ = −ρ −1 ∂ρV (ρ). 9 If we are in a broken-symmetry, static state with ∂ρV (ρ)|ρ =0 = 0 and ∂τρ = 0, the effective equations of motion G L k (seen more easily withǦ L k ) and G T ± k are equal. Equal effective masses and similar behaviors for G L k and G T ± k were assumed for example in [85,86]. However, since during preheating ∂ρV (ρ)|ρ =0 = 0 and ∂τρ = 0 and hence G L k and G T ± k have different effective masses (although the masses still remain comparable). Note that there is no qualitative difference between the longitudinal and transverse modes. As the universe expands, a rough estimate of the amount of particle production in a given Fourier mode can be made by flowing across this plot towards the origin along the white lines (momentum k on the horizontal axis redshifts as a −1 and the inflaton field amplitude decays on the vertical axis as ∼ a −3/2 ). Note that the "local" Floquet exponent µ k is large compared to the instantaneous expansion rate H: [µ k ]/H ∝ m Pl /m in the light colored bands, with m Pl /m ∼ 10 6 from CMB observations. However, whether a given Fourier mode passes through these bands depends on g A . To see this note that at the beginning of preheating |φ| ≈ 0.1m Pl , hence g A |φ|/m 10 5 g A . Hence, for g A 10 −5 , most of the low momentum modes start "below" the lowest resonance bands and never experience significant growth. For g A 10 −5 the unstable momentum modes pass through many Floquet bands as the universe expands and this can lead to significant particle production.

Coulomb gauge analysis
One might be troubled by the presence of singularities in the intermediate steps of the gauge invariant analysis, especially if one wishes to carry out numerical calculations. In the Coulumb gauge the calculation of the physical observables E L k , j 0 k and j L k , is 'clean' throughout, i.e. no singularities arise whatsoever in the intermediate steps.
For more details on the well defined field content and the equations of motion in the Coulomb gauge see Section 2.4 and Appendix B.
The equations of motion for perturbations δϕ 0 k and A T ± k in the Coulumb gauge are identical to the ones for δρ k and G T k respectively in the gauge invariant scenario (withρ →φ 0 ). For a chaotic inflation potential Fig. 3 and are identical to the related Floquet exponents in the gauge invariant case.
In the Coulumb gauge we need to analyse δϕ 1 k instead of the longitudinal mode A L k . For f 1 k = δϕ 1 k in eq. (4.2), the coefficients are Both coefficients are non-singular (assuming φ 0 −1 ∂φ0V (φ 0 ) is well-defined at the originφ 0 = 0 as is usually the case). We can calculate the Floquet exponents based on this equation, however, we shall change variables and remove the damping term for reasons that will become clear below.
We make the non-singular transformation which yieldsb 1 =b L andω 1 =ω L withρ →φ 0 (cf. eq. (4.8)). This identification immediately yields δφ 1 k =Ǧ L k whereǦ L k was provided explicitly in eq. (4.9). Finally undoing our transformation of variables in eq. (4.12), we have the solution with P 1 k ± (τ ) =ľ(φ 0 (τ ))P L k ± (τ ). The Floquet exponent µ L k is identical to the one in eq. (4.9). A numerical computation of this exponent assuming a chaotic inflationary potential is shown in Fig. 3. Also note that now the periodic functions, P 1 k ± (τ ), have no troublesome spikes sinceľ is non-singular everywhere. This is an improvement over the gauge invariant analysis. In that case, G L k =ȟǦ L k , had singular spikes becauseȟ was singular. Hence, the Coulomb gauge allows us to calculate δϕ 1 k (τ ), in a safe, singularity-free way. It is well suited for numerical computations. There is no problem with calculating physical variables such as E L k , B T ± k and E T ± k . For example the longitudinal electric field (cf. Appendix B) 14) has no singularities in it. Before moving on to back-reaction, let us revisit the features seen in Fig. 2 in light of what we now understand from our Floquet analysis. The resonance structure of the transverse mode, eq. (4.5), is identical to that of χ in the popularg 2 ϕ 2 χ 2 toy model [81]. Making use of well-known results about this model (see for example [98]), the parametric excitation of transverse modes is expected to be most efficient in the broad resonance regime g A |φ|/m ∼ g A m Pl /m 1 for modes in the range k m g A |φ|/m ∼ m g A m Pl /m. This should be true for the longitudinal mode as well based on the Floquet plots shown in Fig. 3. The lower three plots in Fig. 2 show how for g A m Pl /m 1, the transverse and longitudinal modes in the range k/H m/H k C /H ∼ k C /H are significantly enhanced, in agreement with the analysis of this section (once expansion is appropriately included).

Back-reaction and end of preheating
So far we have extensively relied on a linear analysis of perturbations. As we have seen, during preheating, the covariant derivative coupling between the gauge fields and the inflaton causes resonant Fourier modes of electric and magnetic fields to grow exponentially fast. This growth cannot proceed forever. The resonance is eventually shut-off by the gauge field's back-reaction on the oscillating inflaton condensate, which ultimately leads to the condensate's fragmentation and our linear analysis stops holding any more.
To estimate the time of back-reaction, we shall investigate the effective equation of motion of the inflaton condensate in the Hartree approximation. We shall work in the non-singular Coulomb gauge, and include the background expansion as well (but ignore metric perturbations). Assuming that the slow-roll inflation happens along the realφ 0 axis, the effective condensate equation in the Hartree approximation becomes During preheating the occupation numbers become much greater than one, and the fields are expected to approach the classical limit. In this limit, the commutators of non-commuting variables vanish and the ambiguity regarding operator ordering is not significant. Moreover, note that the expectation values of non-commuting operators above can be complex. However, taking the real par is a reasonable approximation since in the classical limit the imaginary part becomes negligible compared to the real one (we have verified this numerically as well). Our next step is to derive expressions for the expectation values in eq. (4.15) in terms of mode functions for the different fields. In the Coulomb gauge, and in Fourier space, the temporal component of the gauge field can be written in terms of δϕ 1 k (see Appendix B): and its first derivative, after taking into account the equations of motion, conveniently reduces to These then translate to the mode functions, u 1 k (τ ) and u A 0 k (τ ), of the quantised fluctuations

18)
This immediately yields the following expectation values: (4.20) The only term which involves transverse modes is defined in the Coulomb gauge aŝ (4.22) With these expressions, we can now calculate back-reaction effects iteratively. We first evolve the mode functions from vacuum initial conditions in eq. (3.13), using the background solution ϕ 0 of eq. (4.15) with the correction terms (quadratic in the perturbations) ignored. These then allow us to get the necessary expectation values explicitly. Next, we re-evaluate the background solution, now including the correction terms (quadratic in the perturbations) in eq. (4.15). The moment when the new background solution deviates from the uncorrected one (see the left panel in Fig. 4) back-reaction has become important. For large couplings, g A > 10 −4 in the Chaotic inflation model (V = m 2 ρ 2 /2, m = 10 −6 m Pl ), we get the following expression for the number of e-folds of expansion after the end of inlation when gauge fields back-react on the condensate ∆N br ≈ 1.72g −0.1 A . We caution that this scaling is to be taken as a very rough guide in the range of 3 × 10 −4 < g A < 2.5 × 10 −3 (see the right panel in Fig. 4), the true dependence is likely non monotonic. Below the lowest value of g A we have considered, the Hubble expansion quickly redshifts all modes below the lowest resonance band and we see negligible gauge field production.

Initial conditions for lattice simulations
In this section we give a concise description of initial conditions for numerical lattice studies of reheating. Our linearlized calculations including expansion, metric perturbations and the quantum nature of the fields were sufficient during inflation and preheating. However, after preheating the dynamics of matter fields on subhorizon scales can become highly non-linear. Since (i) the matter field occupation numbers grow rapidly during preheating, (ii) predominantly on subhorizon scales, while (iii) the metric perturbations remain small, one can evolve the classical  There are two distinct parts to setting up initial conditions on the lattice. First, for making the transition from quantized fluctuations to classical ones in the continuum limit, a prescription is needed. The classical field fluctuations must also satisfy the necessary constraints. The second is discretizing these initialized continuous fields on the lattice. We first focus on getting from quantized fluctuations to classical ones, with an accounting for the constraints. We will move from gauge invariant variables to variables in the temporal gauge (a popular choice of lattice simulations), but the prescription is applicable in any gauge.
Working in the local U (1) gauge invariant variables, the matter fields with the quantized fluctuations at the end of preheating, written in terms of the mode functions in Fourier space (in the continuum limit) are as follows. In the expressions below, the time τ is the time of specification of the initial conditions on the lattice (this could be the end of inflation or some time before back-reaction).
We also need time derivatives of the above expressions. In that case, the mode functions u I k get differentiated; the creation and annihilation operators (â I k ,â I † k ) are time independent. By switching from the quantum two point function in real space 0|f I (x)f J (x + r)|0 to the classical 2-point correlation functions in real space we switch from a quantum to a classical description. 10 We now treat the creation and annihilation operators as complex numberŝ with phases distributed uniformly on [0, 2π), and with a Rayleigh distributed amplitudes set in the following manner [99]: Here X I k is a uniform deviate on (0, 1) and Y I k is a uniform deviate on [0, 1). The amplitude is governed by the requirement of consistency with the quantum 2-point function. The index I denotes the field under consideration, i.e. I = {ρ, L, T + , T − }. Note that the complex numbers a L k in eq. (5.1) and eq. (5.4) are the same. So far the analysis have been carried out in a gauge invariant framework. It is not difficult to find the initial conditions for lattice simulations in a particular gauge. For example, in the A 0 = 0 (temporal) gauge -the common gauge choice for numerical simulations, the conversion between gauge invariant variables and gauge dependent fields is given in eqns. (2.45) and (2.48). For convenience one can choose τ in = τ in eq. (2.48), so that δϕ 1 k (τ ) = 0 and A L k (τ ) = G L k , but We note that despite the random nature of the complex numbers in eq. (5.6), the Gaussian constraint on the lattice is automatically satisfied to linear order in the perturbations, since the expression for G 0 (x , τ ) in eq. (5.4), in terms of the complex a L k , takes care of the first order terms in eq. (2.18).
If one wishes the Gauss constraint to be met to machine precision on the lattice, a simple correction has to be made. After initialising δρ(x , τ ), Note the δρ appearing on the right-hand side of the last equation. We then use ∂ τ G L,corr j (x , τ ) in place of ∂ τ G L j (x , τ ) to generate the necessary gauge dependent variable in eq. (5.7). We remind the reader that all expressions so far in this section are for fields on a continuous space-time manifold. The reason is transparency of the transition from quantum to classical fields accounting for the constraints. For completeness we outline the discrete lattice counterparts to the continuous variables. Let ∆ be the separation between neighbouring lattice points and L the length of the lattice. Then the following substitutions should be used to define the fields on the lattice x → x n = n∆ , where n = [n x , n y , n z ] with −L/(2∆) ≤ n x,y,z ≤ L/(2∆). For example There is a minor subtlety in the definition of the polarisation vectors. Consider a finite difference lattice code in which the discrete divergence of a vector is given by 11 Then the continuous directional definitions of the longitudinal and transverse polarisation vectors from eq. (2.34) and eq. (2.38) should be modified to These modifications hold for the spatial discretization stencil in eq. (5.12). For other stencils one can derive similar expressions for the orientation of the polarisation vectors. We stress that if the spatial discretization is not accounted for by the polarization vectors, there will be a violation of the Gauss constraint.
11 E.g. in the standard Wilsonian approach to lattice gauge theories, gauge fields live on space-time links, Gµ(xµ +μ ∆ 2 ) and scalar fields live on the nodes of the space-time lattice, ρ(x µ ). Hereμ is the space-time unit vector and ∆ is the separation between two neighbouring points on the space-time lattice.

The non-Abelian models
We now show that the developed techniques in the previous sections can be used for more complicated non-Abelian models. The main purpose here is to reduce calculations of these complicated models (as far as possible) to the ones we have carried out in the Abelian case. We shall first consider the case of SU (2) non-Abelian gauge fields and then extend to SU (2) × U (1), describing the Electroweak sector of the Standard Model.

SU (2) gauge fields
We consider the following action for a charged scalar doublet The action is invariant under the local SU (2) transformation where the non-Abelian gauge fields are A µ = A a µ σ a /2, and β(x ν ) = β a (x ν )σ a /2, with σ 1 , σ 2 , σ 3 being the three Pauli matrices. The repeated indices are summed over (both for field and spacetime indices).
It is possible to write the above action in terms of local SU (2) invariant fields (analogous to the Abelian case). We will show below that these gauge invariant fields decouple from each other at the linear level and the problem reduces to three identical copies of the Abelian model. The quantisation scheme, preheating analysis as well as setting up of lattice initial conditions discussed in the previous sections then carries over without any difficulty.
We begin by writing the scalar doublet as [98] where ρ will be our real scalar field which forms a condensate during inflation (e.g. the inflaton), and M is a unitary matrix of unit determinant This {ρ, Ω a } decomposition is similar to the polar one give in eq. (2.6) for the Abelian model. In a now familiar way, cf. eq. (2.7), we proceed to define SU (2) invariant non-Abelian fields More explicitly, G µ = G a µ σ a /2, which in component form is The action in eq. (6.1) then simplifies to where the interactions between the gauge bosons are hidden in the definition of the field tensor We will again work at the linear level in G a µ , since the arguments for the vanishing backgrounds of the gauge fields used in the Abelian model apply to this case as well. The inflaton fluctuations (δρ) are decoupled from the gauge fields, G a µ . Furthermore δρ is the only field coupled to the metric perturbations. At linear order we can ignore the interactions between the gauge bosons, i.e. the last term in (6.9). Thereby, the three gauge fields, G a µ , are decoupled from each other and each of them is treated as the gauge field from the Abelian model. The quadratic action for the matter perturbations splits into (cf. eq. (3.1)) Each S I is of the general form given in eq. (3.2). The coefficients b ρ (k, τ ), ω ρ (k, τ ) are those in eq. (3.8), and b aL (k, τ ), ω aL (k, τ ) and b aT ± (k, τ ), ω aT ± (k, τ ) are equal to the longitudinal and transverse ones in eq. (3.9), respectively. Therefore, the inflationary power spectra will be those of three identical copies of longitudinal modes and six identical copies of transverse modes (i.e. one longitudinal and two transverse modes per gauge field). The quantization procedure is identical to the Abelian case with no additional subtleties. The inflaton oscillations during preheating again present a problem for the gauge invariant fields. The gauge invariant fields, G a µ , are ill-defined when Ω a are ill-defined which happens every time ϕ = 0. In analogy with the Abelian case, we will work in the Coulomb gauge. In terms of its "cartesian" components, where we have used the global SU (2) invariance of the action to rotate the internal ϕ axes to align with the direction of motion of the homogeneous field. We have taken this direction to be alonḡ ϕ 0 , with all the other homogeneous components set to zero. 12 In the Coulomb gauge for each 12 Note that Ω a 1, and to linear order in perturbations, the scalar doublet in eq. (6.4) can be written as Let us consider a more realistic scenario in which the scalar has SU (2) × U (1) charges where and F µν (B) and F µν (A) are the field tensors defined in eqns. (2.3) and (6.2), respectively. The coefficients in the covariant derivative are defined in this particular way, so that one can refer the scalar doublet to the Standard Model Higgs field, whose hypercharge is −1/2. The action is invariant under the local U (1) and SU (2) transformations 16) where U was defined in eq. (6.3).
We can write the SU (2) sector (but not the U (1) sector) in terms of gauge invariant variables, ρ and G a µ , by repeating the initial steps (eqns. (6.4)-(6.7)) from Section 6.1. We could have defined SU (2) and U (1) invariant fields. However, following the standard treatment of the Electroweak sector of the Standard Model of Particle Physics, we will fix the gauge for the U (1) sector. With an eye towards the preheating analysis, we will choose Coulomb gauge for the U (1) fields, where ∂ i B i = 0 (in Fourier space, B L k = 0). In addition, instead of proceeding to quantization and then to a preheating analysis, we will work with certain linear combinations of G and B fields to make contact with the Standard Model. To this end, we define the charged W bosons as (cf. eq. (6.7)) The Z boson and the massless photon, A, emerge after rotating in the

17)
This explains the choice of signs and labeling of components in ϕ. The terms on the right-hand side in eq. (6.6) reduce to Recalling eq. (6.7), one arrives at the familiar expression for the gauge invariant fields at linear order in perturbations (cf. eq. (2.7)): G a µ = A a µ + ∇µΩ a . This relation, along with δϕ a = g Aρ Ω a /2 andρ =φ 0 , then allows us to easily derive the equations of motion for A a µ and δϕ a from the equations for G a µ , since G a µ = A a µ + ∇µ(2δϕ a /g Aρ ).
through the Weinberg angle (6.18) The action from eq. (6.14) becomes 19) where the interactions between the gauge bosons are included into the definitions of the field tensors Note that ρ, A µ , Z µ and W ± µ are invariant under SU (2) transformations, eq. (6.3). However, A µ and W ± µ change under a U (1) transformation, eq. (6.16), as follows: That is why, to remove this redundancy, we choose B L k = 0, which implies A L k = Z L k tan θ W , cf. eq. (6.17). This is also a good place to point out that W ± µ are a conjugate pair of complex fields. The '±' should not be confused with the two transverse polarisation states. We work at the linear level in δρ, W ± µ , Z µ , A µ (the arguments for the vanishing backgrounds of the gauge fields used in the previous cases apply here as well). Once again the inflaton fluctuations δρ are decoupled from the rest of the matter fields, and the only ones coupled to the metric perturbations. Neglecting the interactions between the gauge bosons, i.e. the second order terms in eq. (6.20) and expressing W ± µ in terms of G 1 µ , G 2 µ , the second order in perturbations matter action splits into, cf. eqns. (3.1) and (6.10), 13 The S I have the form shown in eq. (3.2). For I = ρ, b ρ (k, τ ) and ω ρ (k, τ ) are given in eq. (3.8). For I = aL, aT ± (with a = {1, 2}), the coefficient pairs {b aL (k, τ ), ω aL (k, τ )} and {b aT ± (k, τ ), 13 While defined W ± µ because of our desire to connect to the Standard model, we write the action for the perturbations in terms of G 1 µ and G 2 µ , because it brings part of the action in a form which looks like multiple copies of the Abelian case. ω aT ± (k, τ )} are the longitudinal and transverse coefficients in eq. (3.9), respectively. The other coefficients are

(6.24)
Notice that only the transverse modes of the photon, A, contribute to the action. The longitudinal modes of A do not contribute since the photon is massless -a manifestation of the Ward identity.
We can now quantise the transverse modes of A and the longitudinal and transverse modes of Z, G 1 , G 2 as done in Section 3.1.
We remind the reader that the procedure outlined above is well-defined during inflation, when ρ 2 /2 = ϕ † ϕ = 0. Similarly to the Abelian and the SU (2) cases, G µ , are ill-defined when ρ = 0 and this can generally happen during preheating. Once again the right prescription is to work in well-defined variables in a non-singular gauge. In the Coulomb gauge, B L = 0 and A aL = 0, the longitudinal modes G 1L , G 2L and Z L can be expressed in terms of the welldefined perturbations of ϕ: δϕ 1 , δϕ 2 and δϕ 3 respectively. In this gauge, the transverse modes are even easier to handle: G 1T ± and G 2T ± reduce to A 1T ± and A 2T ± , respectively; similarly Hence we can keep the expressions for the transverse modes since they have no singularities in them with ρ →φ 0 . As usual, using the SU (2) global invariance, we have rotated our internal axes along the direction of motion of the inflaton.
For purposes of comparison with the literature, we provide the Floquet charts for the longitudinal and transverse modes of the W and Z fields which are excited by the oscillations of the Higgs condensate for the case when V (|ϕ|) = λ ϕ † ϕ 2 = λρ 4 /4, see Fig. 5. The longitudinal modes are excited parametrically during reheating. The inflaton decay to transverse modes of the gauge fields mimics the well studied ϕ 2 χ 2 scalar field model [100]. Without the longitudinal modes, the back-reaction of the transverse gauge fields would be identical to the one in ϕ 2 χ 2 . However, recent numerical experiments [77] indicate that there are small differences between the behaviour of the inflaton decaying to scalars and to Abelian gauge fields. This might be due to the longitudinal modes being excited and playing a role in the back-reaction process. In general, calculations during (p)reheating and estimation of the effects of non-Abelian interaction terms are convenient in the Coulomb gauge (or some other welldefined gauge during preheating). The traditional unitary gauge is not a good choice because the inflaton oscillates through the origin. Although as we saw in the Abelian case, the electric and magnetic fields can be well defined even when using gauge invariant variables, the intermediate steps leading up their calculation often involve singular behavior which at the very least leads to numerical unpleasantness. The same comment holds for the Unitary gauge since the equations of motion are identical to those in gauge invariant variables.  Figure 1 in [87]. Our analysis shows a similar qualitative behavior for the longitudinal modes as well (right). For this plot the scalar field potential is taken to be V (ρ) = λ(ϕ † ϕ) 2 = λρ 4 /4. We use the same notation as in [87]: κ = k/ 2λ|φ| 2 = k/ λρ 2 , q W = g 2 A /4λ and q Z = (g 2 A + g 2 B )/4λ. The amplitude of the oscillating scalar field and the momentum redshift as a −1 , rendering κ redshift-independent. That is why the flow lines from Fig. 3 do not appear here. As the universe expands, Fourier modes do not cross through multiple bands. This is a manifestation of the conformal nature of the quartic potential after the end of inflation.

Observational consequences
Particle production during and after inflation in models with a charged inflaton coupled to gauge fields may leave potentially observable signatures such as primordial magnetic fields, charge fluctuations as well as scalar and tensor metric perturbations. We schematically discuss each one in turn.

Magnetic fields
At the end of inflation the magnetic field power spectrum on superhorizon scales is blue and given by a power-law, see eqns. (3.19) and (3.20). If the gauge fields are 'lighter' than the comoving Hubble scale during inflation, i.e. for k C H we have ∆ B T ± ∝ k 2 . If the gauge fields are "massive" enough, k C H, the power spectrum is much steeper: ∆ B T ± ∝ k 5/2 . After inflation, the oscillating inflaton can resonantly amplify the gauge fields. For the m 2 |ϕ| 2 model considered in Section 3.2, the resonance is broad and effective for k C /a ∼ g A |φ| m, whereφ is the vev of the inflaton during the early stages of preheating at the end of inflation. 14 The above considerations show that the optimal range for production of magnetic fields on superhorizon scales is k C H during inflation and k C /a m during preheating. However, for m 2 |ϕ| 2 inflation, k C /(am) 1 during preheating implies k C /H 1 during inflation, thus excluding the possibility of broad resonance and a relatively shallow power law for the magnetic 14 In the more popular notation, the parameter determining the strength of the resonance, q, is given by q = k C /(am). The resonance is broad when q 1 and narrow or absent when q 1.
field spectrum. We reached the same conclusion via a more detailed calculation whose results were summarized in Fig. 2. We have to settle for two separate regimes: ∆ B T ± ∝ k 2 which receives no resonant amplification during preheating, or ∆ B T ± ∝ k 5/2 which can be resonantly amplified until it backreacts. In the contemporary universe, observations indicate a magnetic field strength of ∼ 10 −7 G on 1 Mpc scales (for example, see [67]). If a seed magnetic field B seed 10 −25 G is present on comoving scales corresponding to 1 Mpc at the time of matter-radiation equality, the galactic dynamo mechanism can potentially amplify B seed to the observed values today [67]. As we show below, the seed field generated in the models under consideration are far too small compared to what is required observationally.
For the case when k C H during inflation, the magnetic field has a shallow spectrum ∝ k 2 on superhorizon scales but is not resonantly amplified during preheating. At the time of decoupling we obtain B seed ∼ ∆ dec B T ± ∼ (2π) −1 (k/a dec ) 2 ∼ 10 −50 G (for example, see [62]) where we used a dec ∼ 10 −3 , k = few Mpc −1 ∼ 10 −38 GeV and 1 GeV 2 ∼ 10 20 G . When k C H, the magnetic field can be resonantly amplified, with an amplification factor that can be as large as ∼ 10 4 − 10 5 . The maximum value is set by back-reaction considerations. However, there is a suppression factor k/k C (see eq. (3.19)) which turns out to be more important. The net effect is that B seed is suppressed compared to the no-resonance scenario. For the case when k C = 10 3 H at the end of inflation (shown in Fig. 2), B seed ∼ 10 −58 G. For this estimate we assumed that the back-reaction takes place 2.5 e-folds after inflation, and subsequently the magnetic field redshifts in the usual way: a 4 ∆ 2 B T ± = const for another 57.5 e-folds. We note that while broad resonance does not happen for k C H for the m 2 |ϕ| 2 models, it can occur for k C < H for steeper potentials (e.g. λ|ϕ| 4 ). It might be possible to boost the amplitude of the seed field up to 10 −45 G in such cases. However, this is still too small to be observationally relevant. Successful magnetogenesis from reheating has been recently discussed in models different from ours (see for example, [89,90]).

Charge fluctuations
Just like the magnetic field, the charge fluctuations also have a blue power spectrum at the end of inflation. Recalling Gauss law, −kE L k = a 2 j 0k , we arrive at the useful expression ∆ j 0 = (k/a)∆ E L . Indeed eqns. (3.19) and (3.20) imply that on superhorizon scales ∆ j 0 ∝ k 2 and k 5/2 for k C H and k C H respectively, reminiscent of the magnetic field spectrum. After inflation, parametric resonance can amplify the charge fluctuations in the m 2 |ϕ| 2 model, provided we are in the broad resonance regime g A |φ| m (i.e. k C H with a spectrum ∝ k 5/2 on superhorizon scales). Note that unlike the magnetic field case, we see some production of charge fluctuations on superhorizon scales at the end of inflation even when m g A |φ|. This fluctuation production is related to the oscillation of the inflaton, but cannot be analysed with simple Floquet analysis since it occurs on superhorizon scales where expansion plays a significant role.
If we imagine that the gauge field in question is the electromagnetic field and put g A = 2e one can compare the charge fluctuations to the observational bounds from vorticity and cosmological magnetic fields on the electrical charge asymmetry of the universe [101], namely that ∆ j 0 /(en B ) 10 −26 on co-moving scales corresponding to 10 2 kpc today, where n B ∼ 1.5×10 −10 T 3 is the number of baryons and T is the photon temperature. For k C < H, we can assume the universe reheats immediately after the end of inflation and the charge fluctuation is transferred to a charge asymmetry in the ordinary matter. 15 The charge asymmetry will then redshift as a 3 ∆ j 0 = const, akin to number density. Then the fractional charge asymmetry ∆ j 0 /en B ∼ 10 9 × (k/a rh T rh ) 2 (g A |ϕ inf |/T rh ) ∼ 10 −36 , where we used |ϕ inf | ∼ m Pl , k = 10 Mpc −1 , T end ∼ 10 16 GeV and a end T end ∼ T 0 = 2.5 × 10 −4 eV. 16 For the broad resonance regime, k C H, the amplitude of the charge fluctuations can be parametrically amplified up to 10 3 − 10 4 . However, there is also a suppression factor which dominates on the scales of interest. Explicitly, for k C = 10 3 H discussed in Fig. 2, if we assume that the universe reheats immediately after back-reaction has taken place (2.5 e-folds after the end of inflation) and then the universe expands for another 57.5 e-folds, the suppression factor is roughly ∼ exp(−2.5)k/k C , implying ∆ j 0 /(en B ) ≈ 10 −45 .
Again, we recall that in the m 2 |ϕ| 2 model, broad resonance does not happen for k C H. For steeper power-law potentials, parametric resonance is efficient for light gauge fields and the excess of charged particles per baryon can go up to 10 −32 . We once again caution the reader that the calculations here are meant to be schematic.

Metric perturbations
The production of gauge fields may affect the primordial metric perturbations, giving rise to (for example) non-gaussianities and gravitational waves. Since the gravitational production of gauge fields during inflation is not significant, they are not expected to give rise to any interesting signals. However, if the post-inflationary dynamics include parametric amplification of fields then interesting signatures are possible. For example, light scalar fields might develop a non-zero vev across the observable universe today along with perturbations around this vev within each Hubble patch at the time of preheating. The details of the following non-linear stage might depend on the vev in each separate universe. For instance, in [102][103][104] it was shown that as the inflaton decays to a light χ field in φ 2 χ 2 models, the back-reaction depends on the initial vevχ i within the Hubble patch. This in turn affects the expansion history within the Hubble patch and therefore may lead to non-gaussianities.χ i also affects the gravitational waves produced during the non-linear evolution within the separate patches. This may lead to low-multiple corrections to the stochastic gravitational wave background [105,106].
For k C H one can show that the gauge field sector has a light scalar with a scale invariant power spectrum -it is δϕ 1 in the Coulomb gauge, see eq. (2.46) and Appendix A. This implies that δϕ 1 develops a non-zero vev across the sky with deviations from it within each Hubble patch at the time of preheating. Its back-reaction should lead to similar effects to the ϕ 2 χ 2 models, since the resonance structure is similar. Note that for the m 2 |ϕ| 2 model we have been considering, backreaction cannot happen because resonance is inefficient if k C H (which was necessary for developing a vev for δϕ 1 during inflation). The absence of resonance for k C H is a feature of m 2 |ϕ| 2 model. For steeper potentials (e.g λ|ϕ| 4 ), broad resonance can occur for k C H. 15 Note that at the end of inflation this charge would correspond to longitudinal electric fields of enormous strength. This may lead to Schwinger pair production, which can annihilate the charge on this scales. For k C = H, e∆ E L ∼ 10 −2 GeV m electron . However, this estimate is not reliable, since the vev of the Standard Model Higgs is largely uncertain at the end of inflation. It is not clear what the masses of the fermions are. 16 We caution the reader that this is a very rough estimate.

Conclusions
In this paper, we carried out a self-contained analysis of particle production during and after inflation in models with charged scalars coupled to Abelian and non-Abelian fields. We calculated the power spectra of the produced gauge fields in the regime where the equations of motion could be linearized. We also provided a prescription for setting up initial conditions for lattice simulations to further evolve the fields nonlinearly.
To make our treatment self-contained, we provided the necessary equations for linear perturbations (including metric fluctuations) in terms of gauge invariant variables, and provided a dictionary to relate these variables to those in some common gauges. This should allow results to be translated between gauges quite easily. For pedagogical purposes, we carried out the initial analysis for the Abelian model. In the later half of the paper we provided the explicit generalization to the non-Abelian cases.
We carried out the quantization and evolution of the perturbations during inflation in terms of gauge invariant variables. After substituting the constraint equations into the original action, the quantization and subsequent evolution was carried out in a straightforward manner. Gauge invariant variables have the advantage of automatically dealing with the correct number of degrees of freedom; we never having to worry about spurious gauge modes. We numerically calculated the electric and magnetic field power spectra at the end of inflation for a simple m 2 |ϕ| 2 inflation (see Fig. 1). We provided an understanding of the shape of the spectra (blue tilt and broken power law behavior -see eqns. (3.19) and (3.20)) in terms of two scales: the co-moving Compton wavenumber of the gauge fields k C ≡ a|φ|g A /2 (φ is the vev of the inflaton condensate during inflation, g A is the coupling strength and a is the scalefactor) and the conformal Hubble scale H. We found that the transverse modes are always dominant over the longitudinal ones on sub-Compton scales, k k C . For gauge fields with k C H the longitudinal modes are as energetic as the transverse ones on super-Compton scales, k k C , whereas when k C H the longitudinal modes dominate on super-Compton scales. The longitudinal mode of the electric field is directly proportional to charge fluctuations, and hence charge fluctuations are also generated during inflation. While our numerical calculations assumed a particular model, these results are expected to hold more generally for potentials where the inflaton is rolling slowly. As a further check, we carried out a semi-analytic analysis in de Sitter universe and confirmed the shapes of the electric and magnetic field spectra.
While gauge invariant variables are suited for calculations during inflation, they become illdefined when the inflaton starts oscillating after inflation. After inflation, the Coulomb gauge turns out to be a particularly well-suited for preheating studies. We carried out an explicit numerical calculation of resonant gauge field production during preheating and provided power spectra of the corresponding electric and magnetic fields at the end of preheating (see Fig. 2). For an analytic understanding of the main features of the spectra, we carried out a Floquet analysis of the instabilities with the gauge constraints included (see Fig. 3). In both calculations, we found that the longitudinal modes are as prone to resonant excitation during preheating as the transverse ones. The characteristic wavenumber range of the instability is k m(g A m Pl /m) 1/2 , whereas the growth rate of fluctuations is determined by g A m Pl /m. We also estimated that the gauge fields back-react on the scalar condensate within about 10 oscillations in the case of a m 2 |ϕ| 2 potential and for couplings, g A ∼ > 10 −3 (see Fig. 4). For g A 10 −4 gauge fields are not produced significantly enough to back-react.
If back-reaction becomes important, the occupation numbers in the gauge fields are usually very high. Further evolution, which is usually highly nonlinear, can be investigated by numerical lattice simulations within the classical approximation. We provided a scheme for initializing such simulations with the gauge constraints being violated only at second order in the perturbations. We also provide a work-around enabling the gauge constraints to be satisfied more precisely at the expense of insignificant (second order in perturbations) modification to the spectrum of the gauge fields. Our prescription for initial conditions is not restricted to a particular gauge; we provide an explicit example of setting up such conditions in the temporal gauge which is commonly used for simulations (see Section 5). The derived power spectra as well as the initial conditions prescription might be useful for lattice simulations of Abelian and non-Abelian fields at the end of inflation (for example, [50, 70-72, 74, 76-78, 91, 94]).
We considered two non-Abelian models. In the SU (2) model we showed that to linear order in perturbations, the problem conveniently reduces to three decoupled identical replicas of the Abelian model. We then extended the local group to SU (2)×U (1), with the Electroweak sector of the Standard Model of Particle Physics in mind. Again, to linear order in the perturbations, the problem splits into three copies of the Abelian model. There are two identical copies, describing the evolution of the massive W bosons, and a similar one for the massive Z boson, the only difference being the coupling constants. The massless photon A remains decoupled from the other sectors of the theory. The framework describing the longitudinal and transverse modes in the Abelian model can be straightforwardly applied to the W and Z sectors, including the scheme for initialising lattice simulations. We provided Floquet charts to capture the instability of these fields during preheating, assuming the scalar condensate to be characterised by λ(ϕ † ϕ) 2 self-interaction (see Fig. 5). In this case, we again confirm that the longitudinal and transverse mode can be excited at a comparable level.
Our analysis allowed us to estimate the magnetic field and charge fluctuations from inflation and preheating in a self-consistent manner. We found that the charged inflaton produces magnetic field that cannot exceed ∼ 10 −50 G on 1 Mpc scales today (consistent with [81]) for m 2 |ϕ| 2 potential. Parametric resonance in models with steeper potentials, e.g. λ|ϕ| 4 may boost the magnetic field by a factor of 10 5 , which is still not enough to seed an efficient galactic dynamo mechanism. The charge fluctuations are shown to be also below the observational bounds, if the inflaton is assumed to be electrically charged. The excess charged particles per baryon are found to be at most 10 −36 on 0.1 Mpc scales for m 2 |ϕ| 2 potential. Preheating in models with steeper potentials, e.g. λ|ϕ| 4 may amplify this by a factor of 10 4 . This is well within the observational bound of 10 −26 [101]. The possibility of an electrically charged inflaton is not ruled out.
Our analysis also reveals that the asymuthal degree of freedom of the charged inflaton develops a scale-invariant power spectrum if the gauge fields are lighter than the Hubble scale. This gives us the possibility of generating non-Gaussianities and gravitational waves [50,[102][103][104][105][106] from the non-linear stage of reheating which we leave for future work.
in Fig. 1 (cf. eq. (3.19)). Secondly, the difference in the k C dependence implies that at late times (k k C ) there will be more power stored in the form of transverse electric field and less in the form of magnetic field. The product of the transverse electric and magnetic power spectra, however, is unchanged.
We repeat the same procedure for the magnetic and transverse electric fields in the weak coupling regime, k C /H < 1 2 , i.e. 1 4 > z 2 > 0. Initially, when k k C , u T ± k ∼ e −ikτ / √ k and ∂ τ u T ± k ∼ √ ke −ikτ . Later on, k k C , u T ± k ∼ 1/ √ k. That is why in this coupling regime the magnetic power spectrum is described by a single power-law. The late time behaviour of the ∂ τ u T ± k is slightly more intriguing. For k 2 C /H k k C , ∂ τ u T ± k ∼ √ k. However, when k k 2 C /H, ∂ τ u T ± k ∼ k C / √ k, implying a double power-law for the transverse electric field, at a scale defined by the square of k C . The actual power spectra again agree with what we found in Section 3.2 for the m 2 ρ 2 /2 inflation (cf. eq. (3.20)) The longitudinal mode is harder to approach analytically. It can be solved analytically only for k k C and k k C . In the former (subhorison) limit the equation of motion reduces to The exact solution is given by where we have used the WKB initial conditions, cf. eq. (3.13). Note u L k has two sorts of harmonic terms: ones of constant amplitude and ones that are linear in −kτ (recall k C = −g Aρ /2Hτ in de Sitter space-time). However, ∂ τ u L k has only terms of the latter kind 18 -this will be important when discussing the weak coupling power spectrum of ∆ 2 E L on super-Hubble scales. In the superhorison limit, k k C , the longitudinal mode is governed by This equation has a general solution in terms of hypergeometric functions which is not very illuminative and we shall not give here. It is more straightforward just to solve the full equation of motion eq. (A.1) for u L k and from its evolution to infer ∆ 2 E L . We did the same thing with the transverse modes. 1 We again start with the strong coupling limit, k C /H 1 2 . At early times, k k C , u L k ∼ e −ikτ √ k/k C and ∂ τ u L k ∼ e −ikτ k 3/2 /k C , as shown in Fig. 6. See why the two are similar in eq. (A.6). At late times, k k C , u L k ∼ 1/ k C and ∂ τ u L k ∼ k C . Rewriting the longitudinal electric field power spectrum as ∆ 2 E L ∼ k 3 |∂ τ u L k | 2 /[1+(k/k C ) 2 ] 2 , we recover the familiar k-scalings (cf. eq. (3.19)) The final case we consider is the weak coupling regime of the longitudinal modes, k C /H 1 2 . At the beginning, when k H, u L k ∼ e −ikτ √ k/k C . However, after the k-mode crosses out the Hubble horizon, but is still shorter than the Compton wavelength, i.e. k C k H, u L k ∼ 1/ √ k. This transition upon Hubble horizon exit for sub-Compton modes is not observed for ∂ τ u L k . Instead for all k k C , ∂ τ u L k ∼ e −ikτ k 3/2 /k C , cf. eq. (A.6). Later on when k k C , u L k ∼ 1/ √ k and ∂ τ u L k ∼ k C / √ k. The power spectrum of the longitudinal electric field than becomes Although, the k-scaling of ∆ 2 E L in the weak coupling regime (cf. eq. (3.20)) is accounted for by our calculations in de Sitter space-time, the power excess on super-Hubble scales seen in Fig. 1 is not. This requires a more general consideration in which the time dependences of H and ρ are included.

B The Abelian model in Coulomb gauge
In Coulomb gauge ∂ i A i = 0, the background variables we consider are (φ 0 ,φ 1 ) = (φ 0 (τ ), 0) (using the U (1) gauge symmetry). The linearised perturbations in Fourier space are δφ 0 k , δϕ 1 k , A 0k , A T ± k . Note that we have not chosen a particular space-time slicing, i.e. we are working in diffeomorphism invariant variables. E.g. δφ 0 is given by δρ from eq. (2.12), withρ →φ 0 . The corresponding equations of motion are as follows.
For the background dynamics, the equations of motion are ∂ 2 τφ 0 + 2H∂ τφ 0 + a 2 ∂φ0V (φ 0 ) = 0 , (B.1) and equations of motion for the perturbations in real and imaginary parts of ϕ are given by: