Thermalized Axion Inflation

# Thermalized Axion Inflation

Ricardo Z. Ferreira 111rferreira@icc.ub.edu,   Alessio Notari 222notari@ub.edu
###### Abstract

We analyze the dynamics of inflationary models with a coupling of the inflaton to gauge fields of the form , as in the case of axions. It is known that this leads to an instability, with exponential amplification of gauge fields, controlled by the parameter , which can strongly affect the generation of cosmological perturbations and even the background. We show that scattering rates involving gauge fields can become larger than the expansion rate , due to the very large occupation numbers, and create a thermal bath of particles of temperature during inflation. In the thermal regime, energy is transferred to smaller scales, radically modifying the predictions of this scenario. We thus argue that previous constraints on are alleviated. If the gauge fields have Standard Model interactions, which naturally provides reheating, they thermalize already at , before perturbativity constraints and also before backreaction takes place. In absence of SM interactions (i.e. for a dark photon), we find that gauge fields and inflaton perturbations thermalize if ; however, observations require , which is above the perturbativity and backreaction bounds and so a dedicated study is required. After thermalization, though, the system should evolve non-trivially due to the competition between the instability and the gauge field thermal mass. If the thermal mass and the instabilities equilibrate, we expect an equilibrium temperature of where is the effective gauge coupling. Finally, we estimate the spectrum of perturbations if is thermal and find that the tensor to scalar ratio is suppressed by , if tensors do not thermalize.

Departament de Física Quàntica i Astrofísica i Institut de Ciències del Cosmos (ICCUB)

Universitat de Barcelona, Martí i Franquès, 1, E-08028, Barcelona, Spain

## 1 Introduction

Inflation is a successful paradigm for the Early Universe, which provides consistent initial conditions for the radiation era and a spectrum of cosmological perturbations in agreement with observations. Its most common realization is thought as a cold state, with a very long stage of quasi exponential expansion due to a scalar field slowly rolling down a flat potential. Such a dynamics exponentially dilutes any remnant but should still be able to reheat the universe at the end of inflation and provide the hot radiation era through some couplings of the scalar field to the Standard Model fields.

A plethora of slow-roll models has been studied within this paradigm, with more or less agreement with observational data. Here we wish to explicitly show the viability of a physically more rich possibility, namely that of a hot plasma already present during inflation, therefore unifying inflation and reheating, and opening thus a new class of inflationary models with its own peculiar predictions. One main purpose of this paper is to construct a working model in which the field can indeed dynamically create such a plasma.

Such a possibility has also been considered by [1, 2] under the name of warm inflation, by invoking a dissipation term due to a coupling to a thermal bath of particles333See also [3, 4] for related earlier work.. An obvious difficulty is the need of an exponential production of radiation in order to overcome the exponential dilution without spoiling the slow-roll stage. To achieve this goal we propose a thermalized axion inflation model in which we simply couple gauge-fields to an axion-like field (which we will think of as the inflaton), through an axial coupling. The phenomenology of such a coupling during inflation has been frequently studied in the literature (see [5, 6, 7, 8, 9, 10, 11] for an incomplete list of references).

The interest about this coupling lies in the instability it triggers on the gauge fields equation of motion in presence of a constant field velocity , leading to strong particle production, that starts at wavelengths of , slightly smaller the horizon size. This instability is present already at linear order in ; the deep reason behind this fact is that the Lagrangian couples to a CP-odd (and thus T-odd) term. Another interesting feature of such a coupling is that the gauge field production can become so large that it backreacts on the background and dynamically generates slow-roll even in absence of a flat potential [5, 11]. This happens when the parameter that controls the instability, , is large enough. It is unclear how to reliably compute the behavior of perturbations in such a backreacting regime, since the gauge field can also backreact on perturbations. In absence of backreaction there are bounds on : it has been shown that at large the gauge fields can leave a large non-Gaussian effect in the curvature perturbation of cosmic microwave background [6, 12, 8]. Another important constraint comes from requiring perturbativity in the loop expansion [9]. Other constraints were derived from the overproduction of black holes by the same mechanism, assuming a given evolution for as a function of time and assuming to know the behavior of perturbations in the backreacting regime [13, 7, 14]. In principle one could also think of generating extra tensor modes [15] through this mechanism, but this becomes difficult due to both non-Gaussianity constraints and the requirement of perturbativity [9, 10].

The main idea of the present paper is that since the instability is able to produce an enormous amount of gauge fields during inflation, the cross sections for gauge field scatterings are largely enhanced by the occupation numbers and their rates are able to overcome the exponential dilution, thus, naturally leading to thermalization and formation of a hot plasma. The fact that all the dynamics is generated by axion-like particles and gauge fields, both having strong protecting symmetries, makes the whole setup well protected. In particular, the axial coupling respects the shift symmetry of and so all induced quantum and thermal corrections should involve derivatives of and so cannot affect the axion potential even when the field thermalizes444In the U(1) case this is exact, while in the non-abelian case non-perturbative effects can break this symmetry and only leave a discrete shift symmetry.. This fact is of great importance as it means that this setup does not introduces new -problems555Note that this the scenario is still sensitive to possible corrections coming from irrelevant operators, as much as, for example, the standard natural inflation setup [16]. Such operators, if present, could become less problematic in a strongly backreacting case since inflation could happen on a steeper potential, due to a friction effect [17, 5]. However, the study of such a regime goes beyond the scope of this paper..

When thermalization is reached, energy moves from the horizon to smaller scales thus completely changing the predictions of this scenario. In fact, we will show that at large the inflaton perturbations become thermal inside the horizon, therefore changing the standard vacuum prediction at horizon crossing and, as a consequence, the predictions for cosmological observables. However, as we will see, this happens when backreaction and higher loop corrections become important and so a dedicated study is needed.

Nevertheless, gauge field thermalization can also happen in a perturbative and non-backreacting regime. That is the case if one considers a, probably more realistic, scenario where the gauge fields belong to the Standard Model (SM). In this case the system can easily reach a partially thermalized state, at lower , in which gauge fields are thermal, while the inflaton perturbations become thermal only at higher . In this case reheating of the universe is completely unified with inflation and no extra ingredient is needed: radiation era simply starts when the potential driving inflation becomes subdominant compared to radiation energy density [18, 19, 11]. Thus, an interesting question to ask is whether the phenomenological constraints on due to non-Gaussianity can be alleviated. In the thermal regime, because of this transfer of energy to smaller scales, we have good reasons to think that a phenomenologically viable window can exist, although in the present paper we will only give qualitative arguments for this to be the case and postpone a full analysis to future work.

But this is not the end of the story. Even though thermalization can be obtained from the initially large occupation numbers the subsequent evolution can be very non-trivial as a result of the competition between the instability, which also becomes less efficient after thermalization, and the generation of thermal masses for the gauge field. Although this subsequent dynamics requires a dedicated study including all these important effects, we argue why the instability and the presence of thermal masses should balance each other and either lead to some periodic behavior or to some stationary stage at an equilibrium temperature . Interestingly, in the latter case we derive very interesting predictions for the spectrum of scalar and tensor perturbations.

The structure of the paper is as follows: in section 2 we summarize the features of the model; in section 3 we study the onset of thermalization, by writing Boltzmann-like equations for scatterings and decays, both with the axial coupling and the SM couplings; in section 4 we solve numerically the set of Boltzmann equations to verify our expectations; in section 5 we study the phenomenological constraints and predictions of the thermalized system; in section 6 we discuss how the presence of thermal masses affects the subsequent evolution of the system; finally in sec. 7 we draw our conclusions. In appendix A we present some additional material related to the numerical solutions as well as some additional derivations.

## 2 The axial coupling

Let us consider an axion-like scalar field with a potential , coupled to a gauge field , with field strength , via an axial coupling, as described by the action

 S=∫d4x√−g[M2p2R+12∂μϕ∂μϕ−V(ϕ)−14FμνFμν−ϕ4fFμν~Fμν]. (2.1)

Here is the dual of . Note that the above coupling induces a periodic potential, of period , only in the non-abelian case, while in the abelian case can be completely unrelated to . Actually, in the context of this paper the only relevant ingredient is the axial coupling to gauge fields, and so our considerations apply also to the case of a non-periodic potential.

We split the scalar field into a spatially homogeneous background value, which drives inflation, and a perturbation, and we approximate the metric with a de-Sitter form in conformal coordinates , where and goes from (past) to (future). Here is the Hubble constant during inflation and a prime denotes the derivative taken with respect to conformal time . In such a background, the gauge field of comoving momentum satisfies the following equation of motion in Coulomb gauge (see [20] and references therein):

 A′′±+(k2±2kξτ)A±=0,ξ≡˙ϕ2fH, (2.2)

where denotes the positive and negative helicity. It is immediate to see that the axial coupling triggers an instability at low in one of the gauge field polarizations, , controlled by the dimensionless parameter . Instead, the other polarization, , gets a different dispersion relation although still positive.

We consider the simple case , which is a good approximation if is slowly-rolling down its potential. It is also a good approximation if we are in the regime of strong backreaction of the gauge fields on  [5, 11]. In both cases , where is the first slow-roll parameter. The equation of motion then has analytical solutions that can be written for example by a Whittaker function:

 A+(k,τ)=1√2keπξ/2W−iξ,1/2(2ikτ), (2.3)

which, in rough terms, connects the flat space oscillatory regime deep inside the horizon666Note however that even at the mode functions always have a logarithmic time-dependent phase. with an exponential growth for , until the solution approaches a constant well outside the horizon, when .

## 3 Thermalized Axion Inflation

The instability created by the axial coupling leads to a strong production of gauge fields , that starts already at momenta larger than . The goal of this section is to point out that the scatterings and decays involving the gauge field are strongly enhanced by the large particle number and can, in some cases, change dramatically the predictions studied so far in the literature. We will analyze two cases: (1) the gauge field is only coupled to , i.e. the gauge field is a “dark photon”, not part of the SM; (2) the gauge field belongs to the SM (either abelian or non-abelian) with known interactions to charged particles. The latter case is probably more realistic, because to have successful reheating we would anyway need to couple the system to the SM.

Crucially the scattering probabilities are Bose-enhanced if the number of is large, which will happen when . As a consequence, if the particle number overcomes some threshold, which depends on and , the scatterings will lead to an equilibrium distribution and as a result the modes will be redistributed from to a new scale given by the temperature . Roughly speaking, when the rates associated to scatterings are larger than we expect the fields to be in a thermal distribution. Such a thermalization involves both polarizations of the gauge field, and , as well as or other charged particles (depending on the size of the cross sections) and would lead, a priori, to a Bose-Einstein (BE) distribution defined as

 NBE(k)=(e−ω(k)aT−1)−1, (3.1)

where is the comoving energy of the particle. Note, however, that , due to the imaginary dispersion relation modes and the associated sourcing of low momenta modes , would be infinitely populated if one could wait an arbitrarily long time. In our inflationary background each mode spends a finite amount of time in the instability region and so it never reaches an infinite occupation number but, as a remnant of such dynamics, we expect relevant deviations from a BE distribution, such as the presence of a peak at low momentum. On the other hand, should still be described by a BE distribution with a modified dispersion relation with , where while the fluctuations, if thermalized, are instead described by a massless BE distribution.

In absence of collisions the energy density in the gauge fields is given by [21], due to the continuous excitation of modes, as described by eq. (2.2). A simple estimate of the expected temperature at thermalization is roughly given by . Note that in our system three scales of interest are present: the horizon scale , the instability scale (each mode starts getting excited when its momentum satisfies ) and the temperature at thermalization . If there is a hierarchy of scales: . We will consider in the range 110: in fact, if it gets very large backreaction will start and one should consider the dynamics of the background together with the mode evolution [5, 11], which we postpone to future work. In any case, even in presence of backreaction should be related logarithmically to the parameters of the potential and typically be at most  [5].

After thermalization is reached, however, the system can evolve in a non-trivial way. In fact, having an interacting plasma typically implies that gauge fields have thermal masses, which should screen the instability and, as a result, quench the particle production, so that the temperature is expected to decrease. Moreover, the energy extraction from the scalar field is proportional to [5, 11]; such a quantity should be suppressed if and are in equilibrium, since they tend to compensate each other when averaged over a thermal distribution. One could imagine that the system could reach a stationary configuration at a temperature smaller than , or perhaps an oscillatory behavior. We postpone the study of the full evolution after thermalization to future work, although we will comment on some expected features in section 6.

### 3.1 Boltzmann equation and particle numbers

In order to study the dynamics that we have just described we define an effective -dependent particle number for each field and derive the associated Boltzmann-like equations. Written in this form we can then insert the standard flat space scattering terms, multiplied by the appropriate scale factors777This should be a good approximation as long as the scattering rates are larger than the expansion rate , which is precisely the regime we are interested in..

We will only consider subhorizon modes because we expect superhorizon modes to become frozen and not to participate in the collisions. This simplified treatment should give an accurate order of magnitude estimate for the parameters and such that thermalization happens, while we postpone a more rigorous treatment for future work.

For a given field , with mode functions , we define a comoving particle number as the ratio between energy density and energy per particle :

 1/2+NX(k)≡k2|Xk|2+|X′k|22ω(k). (3.2)

Deep inside the horizon, in the Bunch-Davies vacuum, the fields have a clear particle interpretation, since they behave as at , with vacuum particle number . In the scalar field case the previous definition will be valid for the canonically normalized field , which satisfies

 u′′k+(k2−2τ)uk=0, (3.3)

where we neglected, for simplicity, slow-roll corrections. In appendix A.4 we also show that the definition of matches the more standard one given in terms of the Bogolyubov coefficient .

As we mentioned in the beginning of this section, the frequency of a mode inside the instability band is non-trivial and so we need to make an extra assumption in the above definition eq. (3.2), when specifying the form of . In fact, the frequency that appears in the equations of motion is in principle time-dependent and can be imaginary: for the gauge fields and for the scalar. While in the scalar case inside the horizon and so the notion of particle is well defined, that is not the case for . To overcome this problem, while keeping the treatment simple, we assume in the above definition that . For the polarization this is a good approximation, especially deep inside the horizon, since as an order of magnitude . However, as we discussed before, for that might not be a good approximation when , since . Although each comoving mode is redshifted and so only stays a short amount of time in such regions, to circumvent this problem one can imagine that at each point in time, when the possible thermalization of the system is to be analyzed, is instantaneously driven to zero, either by slowing down or by increasing . In that case the particle number becomes well defined for all modes . In fact in our numerical simulation we obtained very similar results by using this approach: inserting a large particle number for the gauge fields, by using as an initial condition the solution of the equation of motion with a source but without collisions, and then evolving the system for short time scales in absence of the source term.

Using the above definition of effective particle number we are able to rewrite the equations of motion, eq. (3.3) and eq. (2.2), as an equation for the number of right-handed gauge fields, , and particles, :

 N′γ+(k) = −4kξτRe[gA(k,τ)]|gA(k,τ)|2+k2(Nγ+(k)+1/2), (3.4) N′u(k) = 4τ2Re[gu(k,τ)]|gu(k,τ)|2+k2(Nu(k)+1/2), (3.5) N′γ−(k) = 0, (3.6)

where and . We have also included the left-handed gauge fields, , which are not sourced and so conserved in the absence of collisions. This set of first order Boltzmann-like differential equations is exact and it has a suitable form to include collision terms. Note, however, that in the presence of collisions we would also need to know what happens to the evolution of and . We assume to be the ones given by the free solution in the absence of scatterings888Note that in the case of free solutions and so . Therefore, what makes the source to be non-zero is precisely the modified dispersion relation. In the case of this effect is negligible inside the horizon and so we could even have neglected the source for the processes we are interested in.. This approximation is well justified if we want to capture the onset of thermalization. For we use the well-known exact positive frequency solution of eq. (3.3)

 u=eikτ√2k(1+ikτ),gu=i(k2τ2+ikτ−1)τ(kτ+i), (3.7)

while for we can use the exact solution eq. (2.3).

### 3.2 Collision terms

Thermalization is triggered by the instability in the gauge fields, which populates the phase space and enhances the collision rates, represented by scatterings and decays. Instead, the instability in the scalars, eq. (3.3), does not play a big role in the thermalization but it is crucial to freeze the perturbations around horizon crossing. We insert, thus, collision terms to the right hand side of eqs. (3.4), (3.6) and (3.5) including scatterings and decays , as:

 N′γ+(k) = −4kξτRe[gA(k,τ)]|gA(k,τ)|2+k2(Nγ+(k)+1/2)+S+++S+ϕ+D+ϕ+S+−, (3.8) N′u(k) = 4τ2Re[gu(k,τ)]|gu(k,τ)|2+k2(Nu(k)+1/2)−S+ϕ−D+ϕ, (3.9) N′γ−(k) = −S+−, (3.10)

where the superscript in the collision terms denotes which particles are involved in the process. For completeness we should also include processes involving and . However, since are not produced by the source, such processes should be not relevant to reach thermalization and so we neglect them.

We summarize here all the assumptions used in deriving the above Boltzmann equations:

• the function in eq. (3.8) is given by the solution in absence of scatterings, which is accurate at least before the onset of thermalization;

• the particle number entering in the collision terms is given by eq. (3.2) with , as discussed above;

• backreaction effects of gauge fields on are neglected for simplicity.

#### Scatterings

In the scatterings we use the flat space cross sections for the canonically normalized fields and use the massless free propagator for both the gauge fields and for . Note that, when using the canonical field the axial vertex gets rescaled by . We do not treat exactly such factors in the diagrams but we simply multiply each scattering operator by an overall and assume no other change in the computation. If thermalization happens relatively fast compared to the expansion rate this approximation should be accurate. We, anyway, never run our numerical codes for more than e-fold of expansion. Since the typical energy is , we always assume , so that we are below the cutoff of the theory. For the same reason one could also require because after thermalization modes of energy are populated; however if this requirement is violated it simply means that one is also exciting modes of other fields, belonging to a UV completion of our effective model (involving for example heavy fermions).

We consider the scatterings shown in fig. 1. Each scattering term has the form

 Sk=1ω(k)∫4∏i=2(d3→ki(2π)3(2ωi))|Mi|2(2π)4δ(4)(kμ+kμ2−kμ3−kμ4)B(k,k2,k3,k4), (3.11)

where are the momenta of the external legs, , , is the matrix element of the process (given in appendix A.3) and are the phase space factors

 B(k1,k2,k3,k4) = N1(k1)N2(k2)[1+N3(k3)][1+N4(k4)]−(k1↔k3,k2↔k4), (3.12)

where will depend on the particle in the process. We also assumed CP-invariance (which is true, since we work at tree level).

From the above phase space factors it is clear why collision rates are enhanced by the particle numbers in the initial and final states. For this reason the dominant process should be the one that involves only external .

#### Decays

In flat space the decay rate of a massive particle of mass with an axial coupling is simply given by , so in the massless limit there should be no decay. Here however things are more complicated: we expect a nonzero decay (and inverse decay) rate due to the tachyonic instability of the or, in other words, due to the energy transfer from the background to the fields.

We estimate such a decay rate by looking at the 1-loop correction to the 2-point function , with gauge fields running in the loop. Using the fact that the two point function is related to the mode functions via , we can relate the loop correction to an increase in the particle number and thus to an inverse decay rate999Note that using the Bogolyubov coefficients to define the particle number, as explained in sec. A.4, one also finds that the two point function is related to the particle number as where is the other Bogolyubov coefficient.

To compute the loop correction we use the in-in formalism by considering the interaction Hamiltonian . After Fourier transforming and integrating over the delta functions according to appendix A.1 we get

 ⟨ukuk⟩loop(τ) = a2⟨δϕkδϕk⟩loop(τ) = a2∫τ−∞dτ′(2π)12(2f)2∫τ′−∞dτ′′∫d3q∣∣→e+(→q)⋅→e−(→k−→q)∣∣2× ×⟨[A′q(τ′)A|→k−→q|(τ′)|→k−→q|δϕk(τ′),[A′|→k−→q|(τ′′)Aq(τ′′)qδϕk(τ′′),δϕk(τ)δϕk(τ)]]⟩+ +perm.,

where we omitted a factor of , are the polarization vectors defined in A.4 and we used . Here […] stands for a commutator and is an expectation value. Then, by taking the time derivative of this expression and by making several approximations, described in appendix A.2, we arrive at

 dNu(k)dτ≃−8ξb(2π)3(2f)2a2k× (3.14) ×∫dqq3|→k−→q|min(q,|→k−→q|)Nu(k)(1+Nγ+(q))(1+Nγ+(|→k−→q|))−Nγ+(q)Nγ+(|→k−→q|)(1+Nu(k)),

where is a number of order related to the angular integrals. Written in this form we can now identify the right hand side with the decay term of eq. (3.9).

### 3.3 Analytical estimates

Before proceeding to the numerical evaluation of the Boltzmann-like system of equations it is helpful to get some analytical estimates for the results. When solving numerically the system further assumptions will need to be made and so it is important to have some estimations to compare with. Of course we should keep in mind that the analytical estimations should be seen as order of magnitude estimates for , and so could have corrections. For the analytical approximation we will use the following assumptions:

• only subhorizon modes participate in the collisions, which is possible if ;

• the collision integrals peak where the is maximal, which means at energies ;

• thermalization happens when the collisions are faster than the Hubble expansion;

For the gauge fields the dominant scatterings are because they are enhanced by the most powers of . After integrating the delta functions over the angles the scattering term has the form (we omit here factors of ):

 S++=1ω(k1)∫dk2dk3cSf4[Nγ+,3Nγ+,4(1+Nγ+,1)(1+Nγ+,2)−Nγ+,1Nγ+,2(1+Nγ+,3)(1+Nγ+,4)], (3.15)

where and the coefficient is a dimension 4 combination of the energies involved in the process. The scattering of is only non-zero in the s-channel, whose matrix element is easy to compute and is given in eq. (A.11). In order to estimate when scatterings become relevant note that the Bose-Einstein factors in eq. (3.15) can be also rewritten as:

 Nγ+,1Nγ+,2+Nγ+,1Nγ+,2Nγ+,3+Nγ+,1Nγ+,4Nγ+,2− −Nγ+,3,Nγ+,4Nγ+,2−Nγ+,1Nγ+,3Nγ+,4−Nγ+,3Nγ+,4. (3.16)

We assume the integrals to peak at energies , somewhere between and , where the modes are more densely populated, with . In this limit the above expression goes as , hence, we estimate

 S++≈ω5βSf4N3γ+. (3.17)

Here is a numerical factor that comes from the collision integrals and the cross sections. The scatterings become relevant when this term is comparable with the left-hand side of eq. (3.8) which is of order (or than the source term, which is of the same order). Therefore, the condition in the particle number to have thermalization is

 Nγ+≫√βSHf4ω5. (3.18)

In the next section we compare this estimate with the numerical results. This can also be translated into a bound on , for a given , by using the fact that, in absence of thermalization, the particle number depends exponentially on . We can estimate the total particle number in the band by numerically integrating , using the exact mode functions and fitting with an exponential. This yields . Using this expression in the above condition for thermalization, and setting , then gives

 ξ≳0.45ln(fH)+2.7, (3.19)

which we will compare with numerical results. From this estimate we see that, even if the scattering efficiency is suppressed when , thermalization can still happen for sufficiently large .

Now we apply the same type of analysis for the decays, which have the form

 D+ϕ = 1ω1∫dk2cDf2[Nu3(1+Nγ+,1)(1+Nγ+,2)−Nγ+,1Nγ+,2(1+Nu3)]), (3.20)

where is a dimension 3 coefficient. Similarly to the scatterings, the decay term can be approximated by

 D≈ω3βDf2N2γ+, (3.21)

where we assumed . Here is the inverse of the pre-factor of eq. 3.14. Comparing with gives

 Nγ+(ω)≫βDHf2ω3⟹ξ≳0.45ln(fH)+4.3, (3.22)

which shows that they are subdominant with respect to the scatterings, eq. (3.19). So there may be a regime in which scatterings are in equilibrium but decays are not, which will lead to the presence of a chemical potential, as we will discuss below.

Finally we can estimate the threshold of thermalization for and . In the case of the collision terms can be estimated as , with , which should be compared to the left-hand side of the Boltzmann equation, of order . This leads to the condition

 Nγ+≫√β−Hf4ω5⟹ξ≳0.45ln(fH)+2.2, (3.23)

which is actually realized even more easily than eq. (3.18), since . The case of is also analogous and leads to

 Nγ+≫√βϕHf4ω5⟹ξ≳0.45ln(fH)+2.4, (3.24)

where and so it is also realized almost at the same time as the above eq. (3.23).

Among the above reactions all the scatterings conserve particle number and so they can lead to kinetic equilibrium, but possibly with a nonzero chemical potential ; only the decays can change the particle number and lead to chemical equilibrium, driving to 0 (we are not considering here other processes such as scatterings, which could also do the same). Therefore, we expect the following behavior if, for a fixed value of , we decrease : (1) first eqs.(3.23) and (3.24) would be satisfied, which means that and should start tracking ; (2) then should reach kinetic equilibrium when condition (3.18) is fulfilled, so that all species should reach a Bose-Einstein distribution, possibly with a nonzero ; (3) finally also the decays go in equilibrium driving to zero, reaching blackbody distributions.

Such a picture would be modified in presence of other interactions (as in the case of SM interactions): in this case gauge fields can reach both kinetic and chemical equilibrium more easily, because of fast processes, as well as number changing processes. However, this might not happen to , if its only interaction is the axial coupling, and we may have the situation in which gauge fields are fully in equilibrium, but not .

### 3.4 Standard Model couplings

So far we have only considered processes involving the axial coupling between and the gauge fields, which are the only relevant ones if the gauge field is a “dark photon”, not belonging to the SM. In this section we consider the case in which the gauge field in the axial coupling belongs instead to the Standard Model (SM), either the U(1) hypercharge or a non-abelian SU(2) or SU(3) gauge boson 101010Depending on the couplings of the Higgs to gravity and in a very high energy regime the Higgs vev could be zero and so all gauge bosons could be massless during inflation, otherwise one can adapt these estimates to other cases, where for instance only photons and gluons are massless.. In fact this is the most interesting situation for two reasons: first it gives a natural way to reheat the universe into SM particles and second such a coupling of the inflaton to the SM may make the model directly testable.

The purpose of this section is to study whether the SM self-interactions and scatterings involving gauge bosons and SM charged fields could help in thermalizing the system. The interesting feature of such interactions is that they are not suppressed by powers of like the previous diagrams. Instead, they are proportional to gauge coupling constants and so, if , they are the relevant interactions for thermalization. Moreover, this makes the scenario more predictive since there is only one parameter () that sets the thermalization condition for the gauge bosons. However, as we mentioned before, in this case the perturbations might not be thermalized, depending on the value of .

We will not include the SM interactions in the numerical evaluation of the Boltzmann-like equations, in the next section, since the purpose of this work is just a proof of principle for thermalization. We leave a deeper and more precise study with all such interactions for future work. Nevertheless, we can still estimate the onset of thermalization due to scatterings with SM particles and compare it with the interactions generated from the axial coupling.

#### Photon-Fermion scatterings

Let us first consider gauge boson-fermion scatterings, through pair production. We consider the case of electron-positron production by photon scattering where the differential cross section is given, in the center-of-mass frame and in the high energy limit, by

 dσdΩγγ→e−e+=2πα2(4ω)2(1+cos2θsin2θ), (3.25)

where , is the electric charge and is the center-of-mass energy. The total cross section has a divergence for small angles which can be regulated by imposing an IR cutoff, in our case given by the Hubble scale , or, in the case of thermalization, by a thermal mass of order . Now we can compare this expression with photon-photon scattering mediated by , which in the center-of-mass frame is

 σϕγγ→γγ=|M|264π2(2ω)2≃14(ωf)4164π2(2ω)2. (3.26)

Comparing the two cross sections we have

 σγγ→e−e+σϕγγ→γγ≃2πα2(4ω)214(ωf)4164π2(2ω)2=128π3α2f4ω4, (3.27)

where we assumed the logarithmic term to give a coefficient of order one and the typical energy of the process to be . It is clear that, if , the cross section for pair production is much larger than the scatterings mediated by . In particular if we extrapolate the Standard Model to high energies the relevant coupling is the U(1) hypercharge, whose changes from at GeV to at GeV. Therefore, the ratio of cross-sections is always above one for any , which is anyway required for the validity of effective field theory.

On the other hand, pair production is less Bose-enhanced than photon-photon scattering by one power of . To estimate the threshold for thermalization we compare with where is roughly the particle number density. Therefore, the thermalization condition is

 Nγ+≫8πα2, (3.28)

where we again assumed . Note that this is quite rough, since we have neglected powers of in and . Moreover, we have just considered one diagram, while in reality we have to also sum over all particle-antiparticle pairs of charged fermions in the standard model. Having this is mind and taking into account fractional charges we get an extra factor of approximately111111Here we only consider pair production terms, while adding also Compton scatterings probably increases the result by roughly another factor of 2, at large occupation numbers. 85. Compared to eq. (3.18), this means that thermalization will happen faster through SM interactions than through the axial coupling for any value of . In particular, using the value of derived in the previous subsection, , and the value of at GeV the system would thermalize when .

We can also compare these scattering rates with the rate of particle production by the Schwinger effect given by  [22, 23, 24], where is the coherent electric field and is the electron (or any charged fermion) mass. If we estimate the electric field as and consider the regime , we find a similar threshold for thermalization in terms of .

#### Gauge-field self interactions

Another case in which thermalization of gauge fields can happen quite easily is if the gauge bosons belong to or . In this case, there will be cubic and/or quartic self-interactions already at tree-level, which are boosted in the Boltzmann equation by the number density of positive helicity gauge bosons. Consider for example a quartic tree level interaction, with cross section , with the gauge group coupling constant. The cross-section of gluon-gluon scattering is121212We neglect factors of order 1 in the total cross-section.

 σgg→gg≃9παs2(4ω)2, (3.29)

where . Thus, the condition for thermalization is

 Nγ+≫1αs. (3.30)

For , whose coupling constant runs from at GeV to at GeV scale 131313Note that this analysis is only valid for gauge-fields in the perturbative regime where the propagator and the equation of motion remains unmodified by self-interactions., this means that gluon scatterings are the most efficient process for thermalization and the system would thermalize also around .

## 4 Numerics

This section is devoted to the numerical evaluation of the Boltzmann-like system of equations described in the previous section. Due to the several approximations we have used, the numerical results do not aim at being precise but instead at giving a rough description of the phenomena. Therefore, the results should be seen as an order of magnitude estimate. We start by listing the main approximations and simplifications used:

• We work with a finite set of momenta. For this reason we need to discretize the -integrals appearing in both the scatterings and decays terms in eqs. (3.8), (3.10) and (3.9). Therefore, we substitute

 ∫dk→Δk∑k,

where is the mode spacing. We typically use 10 discrete modes, multiples of a given .

• We only consider subhorizon modes, where the flat space result and the particle interpretation are a good approximation. We start our simulation when the largest mode is subhorizon, using , and stop it when such mode goes superhorizon ().

• We give as initial condition for the solution of the equation of motion eq. (2.2). For and one can use the vacuum solution and eq (3.7), respectively, although the final results will basically be insensitive to such choice.

• To estimate the temperature of the system at thermalization, , one would need to ensure that the chosen window of modes covers the bulk of the energy, which consists of modes of momenta . However, grows exponentially with and the instability band consists of modes of between and . Thus, for large we would need a very large window of momenta to extract correctly the temperature. For numerical reasons we do not consider such a large box but only a band of modes between and approximately. For this reason we cannot evaluate correctly when is large, although the system still approaches a Bose-Einstein distribution in the chosen window.

• We neglect any backreaction on the scalar field equation of motion. This would be relevant when becomes of order of the derivative of the potential and would require to follow the treatment of [5, 11] and to include in the system the background equation for . We postpone this richer situation to future work.

• As already stressed in section 3.1, we assume the functions , defined in eqs. (3.4) and (3.5), to be always given by the solution of the equations of motion in absence of collisions. This assumption might fail after thermalization, but should be valid until there.

• The expressions for the scatterings are based on flat space cross sections, which is justified if thermalization is faster than the Hubble expansion. The decays are instead based on translating results from the two-point function of in the in-in formalism. A more rigorous treatment that combines non-equilibrium field theory with Boltzmann equations is of course desirable, but we think that our treatment should give a correct picture and a reliable order of magnitude estimate.

Note that we use a massless dispersion relation for all the species we consider. This implies that we will get massless BE distributions, while in reality the gauge fields have nontrivial dispersion relations due to the axial coupling. In fact, as we already stressed in sec. 3, we expect deviations from this limit for modes such that , but we postpone the analysis of such deviations to a forthcoming publication.

• Our criterion for thermalization is to verify if the average difference to a BE distribution is smaller than . The average difference is defined as

 ΔNN≡1Ntot∑kNnorm(k)−Neq(k,T)Neq(k,T), (4.1)

where is the total number of modes, is the Bose-Einstein distribution computed at a temperature extracted from the energy density of the system and is the particle number normalized by the ratio for one given value of . The reason to use is that, for large , the mode is outside the box which means that the estimated temperature is not accurate and so there would be an offset between the equilibrium distribution and the numerical result. Note that if the system thermalizes but decays are not efficient the distribution would have a chemical potential (), since particle number is conserved at the time of thermalization. We find that the chemical potential is small compared to the temperature and it is only visible in the smallest mode for such cases (see appendix A.5). For this reason we discard the first mode in the average difference.

• Some of the plots include regions with , which are close to the UV cutoff. However, here we only want to show that our estimates agree well with the numerical results in a region of parameter space where the numerics works well, so that we can then extrapolate for larger values of and , where the numerical treatment becomes more difficult to perform.

On generic grounds, the numerical results confirm our expectations, i.e., that thermalization occurs when the particle number, controlled by , is larger than a given threshold. If decreases, the threshold becomes lower and the system thermalizes more easily. The thresholds for thermalization are in agreement with eq. (3.19). We separate the results in 2 different cases of interest. First we consider a system with modes only. This is the case where thermalization is the most efficient because the gauge field number is larger. Then, we include scatterings with and and numerically evaluate the threshold for thermalization. In appendix A.5 we also add some plots which show the complete distributions as a function of momenta, at a given time.

Let us start by considering only gauge fields and their scatterings, with matrix elements given in eq. (A.11). In fig. 2 we show the total particle number and the average difference to a thermal distribution , defined in eq. (4.1), for fixed and varying and vice-versa, evaluated at the end of the simulation (when the longest mode is horizon size). We assume the system has thermalized when . When , decreases sharply to roughly for , while for that happens for . This is in agreement with the plots on the right which show that the analytical estimate for the threshold of thermalization, using eq. (3.18), intersects the total particle number at roughly the same values of and . In the upper right plot we also see that the particle number depends exponentially on and is well approximated by , as expected from the analytical estimates in sec. 3.1.

In the second case we consider all scatterings and decays involving and . The presence of decays and of more degrees of freedom leads to a decrease in the particle number which implies less efficient thermalization, i.e, it happens for smaller (larger) value of (). In fig. 3 we show the same plots as in fig. 2. By looking at the average differences to a Bose-Einstein distribution (left plots) one can see that for , lower plots, the two gauge field polarizations thermalize at while for all species thermalize when . Comparing the average differences, on the left, with the particle numbers and estimates, on the right, we see that thermalization requires a value of about 0.6 larger and a value of about two times lower than the estimate.

Such results seems to have parallel in the results obtained in [25] where the authors also found that the power spectrum of both polarizations of the gauge field and of the inflaton perturbations follow each other while, at the same time, there is a transfer of energy to the UV (fig. 7 of [25]).

Finally, we also evaluate numerically the threshold for thermalization by collecting the values of and such that . The threshold is well fitted by

 ξ=0.44log(fH)+3.4, (4.2)

which is close to the estimated threshold in eq. (3.19), up to an offset of 0.7 in . We will use this threshold equation in the next section when discussing the phenomenology.

## 5 Phenomenology

After thermalization is reached, however, the evolution is expected to be non-trivial. Indeed the temperature should decrease and the system could also depart from a thermal spectrum at low momenta. In fact, having a temperature normally generates thermal masses for the gauge fields, which tend to screen the instability and, as a result, the source will be less efficient and the temperature is expected to decrease. Moreover, the energy density extracted from the scalar field is given by , averaged over a thermal distribution; such a quantity should be suppressed if and are in equilibrium, since the scatterings tend to restore parity symmetry. We do not try to estimate such quantity exactly here because, while is easy to compute in equilibrium, the same is not true for , due to its complex frequency . We postpone to future work a refined numerical treatment, including such effects, but we anticipate in this section some of the expected features, under the assumption that thermalization is successful. Generically we can anticipate how the scalar and tensor curvature perturbation, and respectively, should be affected by thermalization:

• Superhorizon conservation:

On superhorizon scales the gauge field mode function goes to a constant, so its energy density decreases as . Because the sourcing of adiabatic curvature perturbation is proportional to this quantity, we expect a negligible isocurvature sourcing of the scalar curvature perturbation in the uniform density gauge, . In appendix A.7 we elaborate on this point. Therefore, it should be a good estimate to evaluate correlators of at horizon crossing. We also disregard here the possibility that a backreaction regime, which could change the predictions on perturbations, can be present at any time during inflation.

• Loop corrections:

In the absence of thermalization strong constraints on were derived from non-Gaussianities [6, 12, 8] and the requirement of perturbativity [9, 10] due to loop corrections involving the gauge fields. We expect these corrections to become smaller in the thermal regime by the fact that energy moves from the horizon size to smaller scales, possibly , and so at horizon crossing the effect on correlators should be smaller.

• Parity symmetry:

On top of the above suppression, should also be suppressed in a thermal environment due the tendency of the scatterings to restore parity symmetry. Because each vertex introduces one power of , loop corrections from the axial coupling to odd correlators are proportional to and so they will also be suppressed.

### 5.1 Power spectrum and tensor to scalar ratio

In this subsection we estimate the power spectrum of curvature perturbation assuming that a thermal regime is reached and under the assumptions stated above. We will have different cases depending on whether perturbations are thermalized or not.

Let us briefly review, first, the standard vacuum case. In this case the mode functions associated with the canonically normalized field are given, in the de-Sitter limit, by

 |uk|2=12k∣∣∣1−ikτ∣∣∣2≃−kτ→012k3τ2, (5.1)

which means that the power spectrum of the scalar curvature perturbation in the uniform density gauge, , in the superhorizon limit, , is [26]

 Pvacζ≡|ζk|2k32π2=∣∣∣Hk32π2a˙ϕuk∣∣∣2=H44π2˙ϕ2. (5.2)

Alternatively, one could have evaluated the subhorizon expression for at horizon crossing, , and plugged it into the definition of , which is constant on superhorizon:

 |uk|2≃−kτ→112k⇒Pvacζ=H2k32π2˙ϕ212ka2∣∣ ∣∣−kτ→1=H4∗4π2˙ϕ2∗, (5.3)

where denotes quantities evaluated at horizon crossing. This last result is actually more general because and are evaluated at horizon crossing and so it also holds in the quasi de Sitter case.

If the inflaton perturbations are not thermalized they will still be given at zero order by the above equations, but they will be sourced at one-loop by thermal gauge field fluctuations. This case is, however, more complicated and we also postpone it to a future analysis.

If the inflaton perturbations are instead thermal, the expected outcome is much simpler by following the same strategy of evaluating at horizon crossing. In this case is related to the particle number through eq. (3.2), which should be reliable slightly inside the horizon. Then, using , we have

 ∣∣uthermk∣∣2=1k(12+Nk)∣∣∣−kτ→1≃1kT∗H∗, (5.4)

which then implies that the curvature perturbation is simply obtained by replacing the vacuum particle number (1/2) with the thermal value :

 Pthermζ=T∗H∗H4∗2π2˙ϕ2∗. (5.5)

Note that this prediction naturally has similarities with those of warm inflation [2, 27] from the fact that both use the thermal particle number. However we do not consider any friction term, induced by the thermal bath, in the equation of motion for . Regarding the tilt of the spectrum, although both cases give a nearly scale invariant power spectrum, in the thermal case the departure from scale invariance has a different functional dependence on the slow-roll parameters. In the vacuum case the spectral index is due to the time variation of and at horizon crossing, while here the ratio at horizon crossing also matters.

Indeed the spectral index is given by:

 ns−1≡dlnPthermζdlnk = −6ϵH+2η+dln(T∗/H∗)dlnk, (5.6)

where we have used the following definitions for slow-roll parameters:

 ϵH≡−˙HH2,η≡ϵH−¨ϕH˙ϕ. (5.7)

The reason why we introduced is to keep the analysis fully general: in the backreaction dominated case 141414Note, however, that in the backreaction case one might also worry that could have a non-negligible contribution from gauge field fluctuations at horizon crossing, which is not taken into account by the above eq. (5.5)., in fact, this parameter differs from the previous definition, .

We can also compute the tensor-to-scalar ratio. Tensor perturbations () have couplings to gauge fields smaller than the scalars and so are more difficult to thermalize. Therefore, while we assume here to be thermal, we analyze first the most likely case of a standard vacuum spectrum 151515We neglect here possible 1-loop contributions to . In principle, in fact, there could also be a region where the tensors are non-thermal but have relevant 1-loop corrections to their 2-point function., for tensors. In this case the tensor to scalar ratio is suppressed, because scalar perturbations are enhanced, giving:

 r≡PhPthermζ=16ϵH∗2T∗. (5.8)

This case is phenomenologically very interesting because it would help polynomial large field models (, with ), but also others models of inflation where and are comparable, to agree with observations. In fact the observational bound usually puts a quite stringent bound on , while in our case if is sufficiently large the bound is relaxed and those models are not anymore in tension with observations.

We also analyze for completeness the case of thermalized tensor modes. Since they behave as massless scalar fields in de-Sitter, the result would be similar to eq. (5.5), i.e., enhanced by over the vacuum case

 Pthermalh=H∗T∗π2M2p. (5.9)

In this case the tensor-to-scalar ratio would remain unchanged, , while the tensor tilt would be non-trivial:

 nT≡dlnPthermhdlnk = −2ϵH+dln(T∗/H∗)dlnk. (5.10)

This regime would actually be very interesting. In fact, normally the tensor tilt is necessarily red (), reflecting the fact that the Hubble constant decreases during inflation (and so ); here, instead, we could find a blue spectrum depending on the behavior of . It is however difficult to have thermalization of tensors, due to their weaker coupling, as we will see in section 6.

### 5.2 Constraints on f and ξ

We start by briefly reviewing the main observational constraints on derived in the literature and then discuss how are they affected by thermalization.

If the inflaton perturbations are not thermal, and in absence of backreaction, several constraints on have been derived. Non-Gaussianity on CMB scales constrains [6, 12, 8]. Moreover, perturbativity on the loop expansion for the cosmological correlators requires [9]

 H2f2e2πξ16π2l<1, (5.11)

where is some loop factor. For and imposing , the previous bound constrains although specific 1-loop computations suggest a weaker bound [9, 10].

Another important threshold corresponds to the backreaction of gauge fields, in which the standard slow-roll regime driven by gravitational friction does not apply anymore. This happens if

 V′(ϕ)≃3H˙ϕ≪⟨F~F⟩f.