Thermalized Axion Inflation
TThermalized Axion Inflation
Ricardo Z. Ferreira ∗ , Alessio Notari † Departament de F´ısica Qu`antica i Astrof´ısica i Institut de Ci`encies del Cosmos (ICCUB)Universitat de Barcelona, Mart´ı i Franqu`es, 1, E-08028, Barcelona, Spain
Abstract
We analyze the dynamics of inflationary models with a coupling of the inflaton φ to gauge fields ofthe form φF ˜ F /f , as in the case of axions. It is known that this leads to an instability, with exponentialamplification of gauge fields, controlled by the parameter ξ = ˙ φ/ (2 fH ), which can strongly affect thegeneration of cosmological perturbations and even the background. We show that scattering ratesinvolving gauge fields can become larger than the expansion rate H , due to the very large occupationnumbers, and create a thermal bath of particles of temperature T during inflation. In the thermalregime, 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 StandardModel interactions, which naturally provides reheating, they thermalize already at ξ (cid:38) .
9, beforeperturbativity 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 ξ (cid:38) . ξ (cid:38)
6, which is above the perturbativity and backreaction bounds andso a dedicated study is required. After thermalization, though, the system should evolve non-triviallydue to the competition between the instability and the gauge field thermal mass. If the thermal massand the instabilities equilibrate, we expect an equilibrium temperature of T eq (cid:39) ξH/ ¯ g where ¯ g is theeffective gauge coupling. Finally, we estimate the spectrum of perturbations if φ is thermal and findthat the tensor to scalar ratio is suppressed by H/ (2 T ), if tensors do not thermalize. Inflation is a successful paradigm for the Early Universe, which provides consistent initial conditions forthe radiation era and a spectrum of cosmological perturbations in agreement with observations. Its mostcommon realization is thought as a cold state, with a very long stage of quasi exponential expansion dueto a scalar field slowly rolling down a flat potential. Such a dynamics exponentially dilutes any remnantbut should still be able to reheat the universe at the end of inflation and provide the hot radiation erathrough 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 agreementwith 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 purposeof this paper is to construct a working model in which the field can indeed dynamically create such aplasma. ∗ [email protected] † [email protected] a r X i v : . [ a s t r o - ph . C O ] N ov uch a possibility has also been considered by [1, 2] under the name of warm inflation , by invokinga dissipation term due to a coupling to a thermal bath of particles . An obvious difficulty is the needof an exponential production of radiation in order to overcome the exponential dilution without spoilingthe slow-roll stage. To achieve this goal we propose a thermalized axion inflation model in which wesimply couple gauge-fields A µ to an axion-like field φ (which we will think of as the inflaton), throughan axial coupling. The phenomenology of such a coupling during inflation has been frequently studied inthe literature (see [5–11] for an incomplete list of references).The interest about this coupling lies in the instability it triggers on the gauge fields equation of motionin presence of a constant field velocity ˙ φ , leading to strong particle production, that starts at wavelengthsof O (( ξH ) − ), 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 solarge that it backreacts on the background and dynamically generates slow-roll even in absence of aflat potential [5, 11]. This happens when the parameter that controls the instability, ξ ≡ ˙ φ/ (2 f H ), islarge enough. It is unclear how to reliably compute the behavior of perturbations in such a backreactingregime, since the gauge field can also backreact on perturbations. In absence of backreaction there arebounds on ξ : it has been shown that at large ξ the gauge fields can leave a large non-Gaussian effectin the curvature perturbation of cosmic microwave background [6, 8, 12]. Another important constraintcomes from requiring perturbativity in the loop expansion [9]. Other constraints were derived from theoverproduction of black holes by the same mechanism, assuming a given evolution for ξ as a functionof time and assuming to know the behavior of perturbations in the backreacting regime [7, 13, 14]. Inprinciple one could also think of generating extra tensor modes [15] through this mechanism, but thisbecomes 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 amountof gauge fields during inflation, the cross sections for gauge field scatterings are largely enhanced by theoccupation numbers and their rates are able to overcome the exponential dilution, thus, naturally leadingto 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 wellprotected. In particular, the axial coupling respects the shift symmetry of φ and so all induced quantumand thermal corrections should involve derivatives of φ and so cannot affect the axion potential even whenthe field thermalizes . This fact is of great importance as it means that this setup does not introducesnew η -problems .When thermalization is reached, energy moves from the horizon to smaller scales thus completelychanging the predictions of this scenario. In fact, we will show that at large ξ the inflaton perturbationsbecome thermal inside the horizon, therefore changing the standard vacuum prediction at horizon crossingand, as a consequence, the predictions for cosmological observables. However, as we will see, this happenswhen 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-backreactingregime. That is the case if one considers a, probably more realistic, scenario where the gauge fields See also [3, 4] for related earlier work. In the U(1) case this is exact, while in the non-abelian case non-perturbative effects can break this symmetry and onlyleave a discrete shift symmetry. Note that this the scenario is still sensitive to possible corrections coming from irrelevant operators, as much as, forexample, the standard natural inflation setup [16]. Such operators, if present, could become less problematic in a stronglybackreacting case since inflation could happen on a steeper potential, due to a friction effect [5, 17]. However, the study ofsuch a regime goes beyond the scope of this paper. γ , and theinflaton, φ .belong to the Standard Model (SM). In this case the system can easily reach a partially thermalizedstate, at lower ξ , in which gauge fields are thermal, while the inflaton perturbations become thermalonly at higher ξ . In this case reheating of the universe is completely unified with inflation and noextra ingredient is needed: radiation era simply starts when the potential driving inflation becomessubdominant compared to radiation energy density [11, 18, 19]. Thus, an interesting question to askis whether the phenomenological constraints on ξ due to non-Gaussianity can be alleviated. In thethermal regime, because of this transfer of energy to smaller scales, we have good reasons to think that aphenomenologically viable window can exist, although in the present paper we will only give qualitativearguments 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 initiallylarge occupation numbers the subsequent evolution can be very non-trivial as a result of the competitionbetween the instability, which also becomes less efficient after thermalization, and the generation ofthermal masses for the gauge field. Although this subsequent dynamics requires a dedicated studyincluding all these important effects, we argue why the instability and the presence of thermal massesshould balance each other and either lead to some periodic behavior or to some stationary stage at anequilibrium temperature T eq > H . Interestingly, in the latter case we derive very interesting predictionsfor 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; insection 3 we study the onset of thermalization, by writing Boltzmann-like equations for scatterings anddecays, both with the axial coupling and the SM couplings; in section 4 we solve numerically the set ofBoltzmann equations to verify our expectations; in section 5 we study the phenomenological constraintsand predictions of the thermalized system; in section 6 we discuss how the presence of thermal massesaffects the subsequent evolution of the system; finally in sec. 7 we draw our conclusions. In appendix A wepresent some additional material related to the numerical solutions as well as some additional derivations.3 The axial coupling
Let us consider an axion-like scalar field φ with a potential V ( φ ), coupled to a gauge field A µ , with fieldstrength F µν , via an axial coupling, as described by the action S = (cid:90) d x √− g (cid:34) M p R + 12 ∂ µ φ∂ µ φ − V ( φ ) − F µν F µν − φ f F µν ˜ F µν (cid:35) . (2.1)Here ˜ F µν = (cid:15) µναβ F αβ / (2 √− g ) is the dual of F µν . Note that the above coupling induces a periodicpotential, of period f , only in the non-abelian case, while in the abelian case V ( φ ) can be completelyunrelated to f . Actually, in the context of this paper the only relevant ingredient is the axial couplingto 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, anda perturbation, φ = φ ( τ ) + δφ ( (cid:126)x, τ ) and we approximate the metric with a de-Sitter form in conformalcoordinates ds = a ( τ )[ dτ − d(cid:126)x ], where a = ( − Hτ ) − and τ goes from −∞ (past) to 0 − (future).Here H ≡ a (cid:48) /a is the Hubble constant during inflation and a prime denotes the derivative taken withrespect to conformal time τ . In such a background, the gauge field of comoving momentum k satisfiesthe following equation of motion in Coulomb gauge (see [20] and references therein): A (cid:48)(cid:48)± + (cid:18) k ± kξτ (cid:19) A ± = 0 , ξ ≡ ˙ φ f H , (2.2)where ± denotes the positive and negative helicity. It is immediate to see that the axial coupling triggersan instability at low k 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 ξ (cid:39) constant, which is a good approximation if φ is slowly-rolling downits potential. It is also a good approximation if we are in the regime of strong backreaction of the gaugefields on φ [5, 11]. In both cases ξ = (cid:112) (cid:15)/ M p /f , where (cid:15) ≡ ˙ φ / (2 M p H ) is the first slow-roll parameter.The equation of motion then has analytical solutions that can be written for example by a Whittakerfunction: A + ( k, τ ) = 1 √ k e πξ/ W − iξ, / (2 ikτ ) , (2.3)which, in rough terms, connects the flat space oscillatory regime deep inside the horizon with an expo-nential growth for − kτ < ξ , until the solution approaches a constant A k = e πξ / (2 √ πkξ ) well outsidethe horizon, when − kτ (cid:46) (8 ξ ) − . The instability created by the axial coupling leads to a strong production of gauge fields γ + , that startsalready at momenta larger than H . The goal of this section is to point out that the scatterings anddecays involving the gauge field are strongly enhanced by the large particle number and can, in somecases, 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 Note however that even at τ → −∞ the mode functions always have a logarithmic time-dependent phase. γ + is large, which willhappen when ξ (cid:38)
1. As a consequence, if the particle number overcomes some threshold, which dependson ξ and f , the scatterings will lead to an equilibrium distribution and as a result the modes will beredistributed from k (cid:46) ξH to a new scale given by the temperature T . Roughly speaking, when therates Γ s associated to scatterings are larger than H 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 othercharged particles (depending on the size of the cross sections) and would lead, a priori , to a Bose-Einstein(BE) distribution defined as N BE ( k ) = (cid:16) e − ω ( k ) aT − (cid:17) − , (3.1)where ω ( k ) = √ k + a m is the comoving energy of the particle. Note, however, that γ + , due to theimaginary dispersion relation modes and the associated sourcing of low momenta modes k/a < ξH , wouldbe infinitely populated if one could wait an arbitrarily long time. In our inflationary background eachmode spends a finite amount of time in the instability region and so it never reaches an infinite occupationnumber but, as a remnant of such dynamics, we expect relevant deviations from a BE distribution, suchas the presence of a peak at low momentum. On the other hand, γ − should still be described by a BEdistribution with a modified dispersion relation with ω = k + m , where m = − kξ/τ 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 ρ R ≈ − H e πξ /ξ [21],due to the continuous excitation of modes, as described by eq. (2.2). A simple estimate of the expectedtemperature at thermalization is roughly given by ¯ T ≈ ρ / R ≈ . He πξ/ . Note that in our system threescales of interest are present: the horizon scale H , the instability scale 2 ξH (each mode starts gettingexcited when its momentum satisfies k/a < ξH ) and the temperature at thermalization ¯ T . If ξ (cid:38) O (1)there is a hierarchy of scales: H < ξH < ¯ T . We will consider ξ in the range 1 ∼
10: in fact, if it gets verylarge backreaction will start and one should consider the dynamics of the background together with themode 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 O (10) [5].After thermalization is reached, however, the system can evolve in a non-trivial way. In fact, hav-ing an interacting plasma typically implies that gauge fields have thermal masses, which should screenthe instability and, as a result, quench the particle production, so that the temperature is expectedto decrease. Moreover, the energy extraction from the scalar field is proportional to ˙ φ/ (4 f ) (cid:104) F ˜ F (cid:105) = a − (cid:82) kd d/dτ ( | A + | − | A − | ) [5, 11]; such a quantity should be suppressed if γ + and γ − are in equilib-rium, since they tend to compensate each other when averaged over a thermal distribution. One couldimagine that the system could reach a stationary configuration at a temperature smaller than ¯ T , or per-haps an oscillatory behavior. We postpone the study of the full evolution after thermalization to futurework, although we will comment on some expected features in section 6.5 .1 Boltzmann equation and particle numbers In order to study the dynamics that we have just described we define an effective k -dependent particlenumber N X ( k ) for each field X and derive the associated Boltzmann-like equations. Written in this formwe can then insert the standard flat space scattering terms, multiplied by the appropriate scale factors .We will only consider subhorizon modes because we expect superhorizon modes to become frozen andnot to participate in the collisions. This simplified treatment should give an accurate order of magnitudeestimate for the parameters f and ξ such that thermalization happens, while we postpone a more rigoroustreatment for future work.For a given field X , with mode functions X k , we define a comoving particle number as the ratiobetween energy density ρ k and energy per particle ω ( k ):1 / N X ( k ) ≡ k | X k | + | X (cid:48) k | ω ( k ) . (3.2)Deep inside the horizon, in the Bunch-Davies vacuum, the fields have a clear particle interpretation, sincethey behave as X k = e ikτ / √ k at τ → −∞ , with vacuum particle number N X = 0. In the scalar fieldcase the previous definition will be valid for the canonically normalized field u ≡ aδφ , which satisfies u (cid:48)(cid:48) k + (cid:18) k − τ (cid:19) u k = 0 , (3.3)where we neglected, for simplicity, slow-roll corrections. In appendix A.4 we also show that the definitionof N X ( k ) matches the more standard one given in terms of the Bogolyubov coefficient N k = | β k | .As we mentioned in the beginning of this section, the frequency of a mode inside the instabilityband is non-trivial and so we need to make an extra assumption in the above definition eq. (3.2), whenspecifying the form of ω ( k ). In fact, the frequency that appears in the equations of motion is in principletime-dependent and can be imaginary: ω , − = ( k ± kξ/τ ) for the gauge fields and ω u = ( k − /τ )for the scalar. While in the scalar case ω u > ω + . To overcome this problem, while keeping the treatment simple, weassume in the above definition that ω + , − (cid:39) k . For the ( − ) polarization this is a good approximation,especially deep inside the horizon, since as an order of magnitude ω − ( k ) (cid:39) k . However, as we discussedbefore, for γ + that might not be a good approximation when − kτ (cid:39) ξ , since ω + ( k ) →
0. Although eachcomoving mode k is redshifted and so only stays a short amount of time in such regions, to circumventthis problem one can imagine that at each point in time, when the possible thermalization of the systemis to be analyzed, ξ is instantaneously driven to zero, either by slowing down φ or by increasing f . Inthat case the particle number becomes well defined for all modes k . In fact in our numerical simulationwe obtained very similar results by using this approach: inserting a large particle number for the gaugefields, by using as an initial condition the solution of the equation of motion with a source but withoutcollisions, 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, N γ + ( k ), and u particles, This should be a good approximation as long as the scattering rates are larger than the expansion rate H , which isprecisely the regime we are interested in. u ( k ): N (cid:48) γ + ( k ) = − kξτ Re [ g A ( k, τ )] | g A ( k, τ ) | + k (cid:0) N γ + ( k ) + 1 / (cid:1) , (3.4) N (cid:48) u ( k ) = 4 τ Re [ g u ( k, τ )] | g u ( k, τ ) | + k ( N u ( k ) + 1 / , (3.5) N (cid:48) γ − ( k ) = 0 , (3.6)where g u ≡ u (cid:48) /u and g A ≡ A (cid:48) + /A + . We have also included the left-handed gauge fields, N γ − ( k ), whichare not sourced and so conserved in the absence of collisions. This set of first order Boltzmann-likedifferential equations is exact and it has a suitable form to include collision terms. Note, however, thatin the presence of collisions we would also need to know what happens to the evolution of g A and g u . Weassume g u,A to be the ones given by the free solution in the absence of scatterings . This approximationis well justified if we want to capture the onset of thermalization. For g u we use the well-known exactpositive frequency solution of eq. (3.3) u = e ikτ √ k (cid:18) ikτ (cid:19) , g u = i (cid:0) k τ + ikτ − (cid:1) τ ( kτ + i ) , (3.7)while for g A we can use the exact solution eq. (2.3). Thermalization is triggered by the instability in the gauge fields, which populates the phase space andenhances 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 perturbationsaround horizon crossing. We insert, thus, collision terms to the right hand side of eqs. (3.4), (3.6)and (3.5) including scatterings S n and decays D , as: N (cid:48) γ + ( k ) = − kξτ Re [ g A ( k, τ )] | g A ( k, τ ) | + k (cid:0) N γ + ( k ) + 1 / (cid:1) + S ++ + S + φ + D + φ + S + − , (3.8) N (cid:48) u ( k ) = 4 τ Re [ g u ( k, τ )] | g u ( k, τ ) | + k ( N u ( k ) + 1 / − S + φ − D + φ , (3.9) N (cid:48) γ − ( k ) = − S + − , (3.10)where the superscript in the collision terms denotes which particles are involved in the process. Forcompleteness we should also include processes involving γ − and φ . However, since γ − are not producedby 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 g A in eq. (3.8) is given by the solution in absence of scatterings, which is accurate atleast before the onset of thermalization; • the particle number entering in the collision terms is given by eq. (3.2) with ω ( k ) (cid:39) k , as discussedabove; Note that in the case of free solutions g ( k, τ ) = − ik and so Re[ g ] = 0. Therefore, what makes the source to be non-zerois precisely the modified dispersion relation. In the case of u this effect is negligible inside the horizon and so we could evenhave neglected the source for the processes we are interested in. 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 themassless free propagator for both the gauge fields and for u . Note that, when using the canonical field u the axial vertex gets rescaled by 1 /a ( τ ). We do not treat exactly such factors in the diagrams butwe simply multiply each scattering operator by an overall 1 /a ( τ ) and assume no other change in thecomputation. If thermalization happens relatively fast compared to the expansion rate this approximationshould be accurate. We, anyway, never run our numerical codes for more than O (1) e-fold of expansion.Since the typical energy is ω ( k ) (cid:46) ξH , we always assume ξH (cid:28) f , so that we are below the cutoff of thetheory. For the same reason one could also require T (cid:28) f because after thermalization modes of energy ω ( k ) (cid:46) T are populated; however if this requirement is violated it simply means that one is also excitingmodes of other fields, belonging to a UV completion of our effective model (involving for example heavyfermions).We consider the scatterings shown in fig. 1. Each scattering term S k has the form S k = 1 ω ( k ) (cid:90) (cid:89) i =2 (cid:32) d (cid:126)k i (2 π ) (2 ω i ) (cid:33) | M i | (2 π ) δ (4) ( k µ + k µ − k µ − k µ ) B ( k, k , k , k ) , (3.11)where k µi = ( ω i , (cid:126)k i ) are the momenta of the external legs, k i = | (cid:126)k i | , k µ = ( ω, (cid:126)k ), M n is the matrix elementof the process (given in appendix A.3) and B are the phase space factors B ( k , k , k , k ) = N ( k ) N ( k ) [1 + N ( k )] [1 + N ( k )] − ( k ↔ k , k ↔ k ) , (3.12)where N i will depend on the particle in the process. We also assumed CP-invariance (which is true, sincewe work at tree level).From the above phase space factors it is clear why collision rates are enhanced by the particle numbersin the initial and final states. For this reason the dominant process should be the one that involves onlyexternal γ + . Decays
In flat space the decay rate of a massive particle of mass m with an axial coupling is simply given byΓ d ∝ m φ /f , 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 otherwords, 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 (cid:104) u k u k (cid:105) ,with gauge fields γ + running in the loop. Using the fact that the two point function is related to themode functions via (cid:104) u k u k (cid:48) (cid:105) = | u k | (2 π ) δ ( (cid:126)k − (cid:126)k (cid:48) ), we can relate the loop correction to an increase in theparticle number and thus to an inverse decay rate .To compute the loop correction we use the in-in formalism by considering the interaction Hamilto-nian H int = (cid:82) d x √− g δφ F ˜ F / (4 f ). After Fourier transforming and integrating over the delta functions Note that using the Bogolyubov coefficients to define the particle number, as explained in sec. A.4, one also finds thatthe two point function is related to the particle number as (cid:104) u k u k (cid:105) ∝ N k +2Re (cid:2) β k α ∗ k (cid:3) where α k is the other Bogolyubovcoefficient (cid:104) u k u k (cid:105) loop ( τ ) = a (cid:104) δφ k δφ k (cid:105) loop ( τ ) (3.13)= a (cid:90) τ −∞ dτ (cid:48) (2 π ) (2 f ) (cid:90) τ (cid:48) −∞ dτ (cid:48)(cid:48) (cid:90) d q (cid:12)(cid:12)(cid:12) (cid:126)e + ( (cid:126)q ) · (cid:126)e − ( (cid:126)k − (cid:126)q ) (cid:12)(cid:12)(cid:12) ×× (cid:68)(cid:104) A (cid:48) q ( τ (cid:48) ) A | (cid:126)k − (cid:126)q | ( τ (cid:48) ) | (cid:126)k − (cid:126)q | δφ k ( τ (cid:48) ) , (cid:104) A (cid:48)| (cid:126)k − (cid:126)q | ( τ (cid:48)(cid:48) ) A q ( τ (cid:48)(cid:48) ) q δφ k ( τ (cid:48)(cid:48) ) , δφ k ( τ ) δφ k ( τ ) (cid:105)(cid:105)(cid:69) ++ perm. , where we omitted a factor of δ (0), (cid:126)e + , − are the polarization vectors defined in A.4 and we used q ≡ | (cid:126)q | .Here [...] stands for a commutator and (cid:104) ... (cid:105) is an expectation value. Then, by taking the time derivativeof this expression and by making several approximations, described in appendix A.2, we arrive at dN u ( k ) dτ (cid:39) − ξb (2 π ) (2 f ) a k × (3.14) × (cid:90) dq q | (cid:126)k − (cid:126)q | min( q, | (cid:126)k − (cid:126)q | ) N u ( k )(1 + N γ + ( q ))(1 + N γ + ( | (cid:126)k − (cid:126)q | )) − N γ + ( q ) N γ + ( | (cid:126)k − (cid:126)q | )(1 + N u ( k )) , where b is a number of order O (10 − ) related to the angular integrals. Written in this form we can nowidentify the right hand side with the decay term D + φ of eq. (3.9). Before proceeding to the numerical evaluation of the Boltzmann-like system of equations it is helpful toget some analytical estimates for the results. When solving numerically the system further assumptionswill need to be made and so it is important to have some estimations to compare with. Of course weshould keep in mind that the analytical estimations should be seen as order of magnitude estimates for f , and so ξ could have O (1) corrections. For the analytical approximation we will use the followingassumptions: • only subhorizon modes participate in the collisions, which is possible if ξ (cid:38) • the collision integrals peak where the N γ + is maximal, which means at energies H (cid:46) ω (cid:46) ξH ; • thermalization happens when the collisions are faster than the Hubble expansion;For the gauge fields the dominant scatterings are γ + γ + ↔ γ + γ + because they are enhanced by themost powers of N γ + . After integrating the delta functions over the angles the scattering term has theform (we omit here factors of a ): S ++ = 1 ω ( k ) (cid:90) dk dk c S f (cid:2) N γ + , N γ + , (1 + N γ + , )(1 + N γ + , ) − N γ + , N γ + , (1 + N γ + , )(1 + N γ + , ) (cid:3) , (3.15) where N γ + ,i ≡ N γ + ( k i ) and the coefficient c S is a dimension 4 combination of the energies ω ( k i ) involvedin the process. The scattering of γ + is only non-zero in the s-channel, whose matrix element is easy tocompute and is given in eq. (A.11). In order to estimate when scatterings become relevant note that theBose-Einstein factors in eq. (3.15) can be also rewritten as: N γ + , N γ + , + N γ + , N γ + , N γ + , + N γ + , N γ + , N γ + , −− N γ + , , N γ + , N γ + , − N γ + , N γ + , N γ + , − N γ + , N γ + , . (3.16)9e assume the integrals to peak at energies ω i (cid:39) ω , somewhere between H and 2 ξH , where the modesare more densely populated, with N γ + ,i ≡ N γ + (cid:29)
1. In this limit the above expression goes as N γ + ,hence, we estimate S ++ ≈ ω β S f N γ + . (3.17)Here β S = O (10 ) 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) whichis of order N γ + H (or than the source term, which is of the same order). Therefore, the condition in theparticle number to have thermalization is N γ + (cid:29) (cid:114) β S Hf ω . (3.18)In the next section we compare this estimate with the numerical results. This can also be translated intoa bound on ξ , for a given f /H , by using the fact that, in absence of thermalization, the particle numberdepends exponentially on ξ . We can estimate the total particle number in the band 1 < − kτ < ξ bynumerically integrating N γ + ≈ k | A + | , using the exact mode functions and fitting with an exponential.This yields N γ + ≈ − e . ξ . Using this expression in the above condition for thermalization, and setting ω ≈ H , then gives ξ (cid:38) .
45 ln (cid:18) fH (cid:19) + 2 . , (3.19)which we will compare with numerical results. From this estimate we see that, even if the scatteringefficiency is suppressed when f (cid:29) H , thermalization can still happen for sufficiently large ξ .Now we apply the same type of analysis for the decays, which have the form D + φ = 1 ω (cid:90) dk c D f (cid:2) N u (1 + N γ + , )(1 + N γ + , ) − N γ + , N γ + , (1 + N u )] (cid:1) , (3.20)where c D is a dimension 3 coefficient. Similarly to the scatterings, the decay term can be approximatedby D ≈ ω β D f N γ + , (3.21)where we assumed N γ + (cid:29) N φ . Here β D ≈ /ξ is the inverse of the pre-factor of eq. 3.14. Comparing D with N γ + H gives N γ + ( ω ) (cid:29) β D Hf ω = ⇒ ξ (cid:38) .
45 ln (cid:18) fH (cid:19) + 4 . , (3.22)which shows that they are subdominant with respect to the scatterings, eq. (3.19). So there may be aregime in which scatterings are in equilibrium but decays are not, which will lead to the presence of achemical potential, as we will discuss below.Finally we can estimate the threshold of thermalization for φ and γ − . In the case of γ − the collisionterms can be estimated as S + − ≈ ω N γ + N γ − / ( β − f ), with β − = O (10 ), which should be compared to10he left-hand side of the Boltzmann equation, of order N γ − H . This leads to the condition N γ + (cid:29) (cid:114) β − Hf ω = ⇒ ξ (cid:38) .
45 ln (cid:18) fH (cid:19) + 2 . , (3.23)which is actually realized even more easily than eq. (3.18), since β − (cid:28) β S . The case of φ is also analogousand leads to N γ + (cid:29) (cid:114) β φ Hf ω = ⇒ ξ (cid:38) .
45 ln (cid:18) fH (cid:19) + 2 . , (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 kineticequilibrium, but possibly with a nonzero chemical potential µ ; only the decays can change the particlenumber and lead to chemical equilibrium, driving µ to 0 (we are not considering here other processes suchas 2 ↔ ξ , we decrease f /H : (1) first eqs.(3.23) and (3.24) would be satisfied, which means that N γ − and N φ should start tracking N γ + ; (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 fast2 ↔ φ , if itsonly interaction is the axial coupling, and we may have the situation in which gauge fields are fully inequilibrium, but not φ . 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 thissection we consider the case in which the gauge field in the axial coupling belongs instead to the StandardModel (SM), either the U(1) hypercharge or a non-abelian SU(2) or SU(3) gauge boson . In fact this isthe most interesting situation for two reasons: first it gives a natural way to reheat the universe into SMparticles 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 involvinggauge bosons and SM charged fields could help in thermalizing the system. The interesting feature ofsuch interactions is that they are not suppressed by powers of 1 /f like the previous diagrams. Instead,they are proportional to gauge coupling constants and so, if f (cid:29) H , they are the relevant interactionsfor 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 thiscase the φ perturbations might not be thermalized, depending on the value of f /H .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 Depending on the couplings of the Higgs to gravity and in a very high energy regime the Higgs vev could be zero andso all gauge bosons could be massless during inflation, otherwise one can adapt these estimates to other cases, where forinstance only photons and gluons are massless.
Photon-Fermion scatterings
Let us first consider gauge boson-fermion scatterings, through pair production. We consider the caseof electron-positron production by photon scattering where the differential cross section is given, in thecenter-of-mass frame and in the high energy limit, by dσd Ω γγ → e − e + = 2 πα (4 ω ) (cid:18) θ sin θ (cid:19) , (3.25)where α = e / (4 π ), e is the electric charge and 2 ω is the center-of-mass energy. The total cross sectionhas a log divergence for small angles which can be regulated by imposing an IR cutoff, in our case givenby the Hubble scale H , or, in the case of thermalization, by a thermal mass of order eT . Now we cancompare this expression with photon-photon scattering mediated by φ , which in the center-of-mass frameis σ φγγ → γγ = | M | π (2 ω ) (cid:39) (cid:18) ωf (cid:19) π (2 ω ) . (3.26)Comparing the two cross sections we have σ γγ → e − e + σ φγγ → γγ (cid:39) πα (4 ω ) (cid:16) ωf (cid:17) π (2 ω ) = 128 π α f ω , (3.27)where we assumed the logarithmic term to give a coefficient of order one and the typical energy of theprocess to be ω . It is clear that, if f (cid:29) H , the cross section for pair production is much larger thanthe scatterings mediated by φ . In particular if we extrapolate the Standard Model to high energies therelevant coupling is the U(1) hypercharge, whose α changes from 1 /
40 at 10 GeV to 1 /
60 at 10 GeV.Therefore, the ratio of cross-sections is always above one for any f > ω , which is anyway required for thevalidity of effective field theory.On the other hand, pair production is less Bose-enhanced than photon-photon scattering by onepower of N γ + . To estimate the threshold for thermalization we compare H with n γ + σ γγ → e − e + where n γ + ≈ N γ + H is roughly the particle number density. Therefore, the thermalization condition is N γ + (cid:29) πα , (3.28)where we again assumed ω (cid:39) H . Note that this is quite rough, since we have neglected powers of ξ in ω and n γ + . Moreover, we have just considered one diagram, while in reality we have to also sum overall particle-antiparticle pairs of charged fermions in the standard model. Having this is mind and takinginto account fractional charges we get an extra factor of approximately
85. Compared to eq. (3.18), thismeans that thermalization will happen faster through SM interactions than through the axial couplingfor any value of f > H . In particular, using the value of N γ + derived in the previous subsection, Here we only consider pair production terms, while adding also Compton scatterings probably increases the result byroughly another factor of 2, at large occupation numbers. γ + = 10 − e . ξ , and the value of α at 10 GeV the system would thermalize when ξ (cid:38) . S (cid:39) e E π H exp (cid:16) − πm e eE (cid:17) [22–24], where E is the coherent electric field and m e is the electron(or any charged fermion) mass. If we estimate the electric field as E (cid:39) N γ + H and consider the regime E (cid:29) m e , 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 belongto SU (2) or SU (3). 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 ∝ g s , with g s the gauge groupcoupling constant. The cross-section of gluon-gluon scattering is σ gg → gg (cid:39) πα s ω ) , (3.29)where α s = g s / (4 π ). Thus, the condition for thermalization is N γ + (cid:29) α s . (3.30)For SU (3), whose coupling constant α s = g s / (4 π ) runs from 1 /
40 at 10 GeV to (cid:39) ,this means that gluon scatterings are the most efficient process for thermalization and the system wouldthermalize also around ξ (cid:39) . This section is devoted to the numerical evaluation of the Boltzmann-like system of equations describedin the previous section. Due to the several approximations we have used, the numerical results do notaim at being precise but instead at giving a rough description of the phenomena. Therefore, the resultsshould be seen as an order of magnitude estimate. We start by listing the main approximations andsimplifications used: • We work with a finite set of momenta. For this reason we need to discretize the k -integrals appearingin both the scatterings and decays terms in eqs. (3.8), (3.10) and (3.9). Therefore, we substitute (cid:90) dk → ∆ k (cid:88) k , where ∆ k is the mode spacing. We typically use 10 discrete modes, multiples of a given k min . • We only consider subhorizon modes, where the flat space result and the particle interpretation area good approximation. We start our simulation when the largest mode k min is subhorizon, using | k min τ | = 2, and stop it when such mode goes superhorizon ( | k min τ | < We neglect factors of order 1 in the total cross-section. Note that this analysis is only valid for gauge-fields in the perturbative regime where the propagator and the equationof motion remains unmodified by self-interactions. 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 basicallybe insensitive to such choice. • To estimate the temperature of the system at thermalization, ¯ T , one would need to ensure thatthe chosen window of modes covers the bulk of the energy, which consists of modes of momenta k/a (cid:46) ¯ T . However, ¯ T grows exponentially with ξ and the instability band consists of modes of k/a between H and ξH . Thus, for large ξ we would need a very large window of momenta to extractcorrectly the temperature. For numerical reasons we do not consider such a large box but only aband of modes between H and O (10) H approximately. For this reason we cannot evaluate correctly¯ T when ξ is large, although the system still approaches a Bose-Einstein distribution in the chosenwindow. • We neglect any backreaction on the scalar field equation of motion. This would be relevant when (cid:104) F ˜ F (cid:105) /f becomes of order of the derivative of the potential dV /dφ and would require to follow thetreatment of [5, 11] and to include in the system the background equation for φ . We postpone thisricher situation to future work. • As already stressed in section 3.1, we assume the functions g u,A , defined in eqs. (3.4) and (3.5), tobe always given by the solution of the equations of motion in absence of collisions. This assumptionmight fail after thermalization, but should be valid until there. • The expressions for the scatterings are based on flat space cross sections, which is justified ifthermalization is faster than the Hubble expansion. The decays are instead based on translatingresults from the two-point function of δφ in the in-in formalism. A more rigorous treatment thatcombines non-equilibrium field theory with Boltzmann equations is of course desirable, but we thinkthat 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 thatwe will get massless BE distributions, while in reality the gauge fields have nontrivial dispersionrelations due to the axial coupling. In fact, as we already stressed in sec. 3, we expect deviationsfrom this limit for modes such that k/a (cid:46) ξH , but we postpone the analysis of such deviations toa forthcoming publication. • Our criterion for thermalization is to verify if the average difference to a BE distribution is smallerthan 1%. The average difference is defined as∆ NN ≡ N tot (cid:88) k N norm ( k ) − N eq ( k, T ) N eq ( k, T ) , (4.1)where N tot is the total number of modes, N eq ( k, T ) is the Bose-Einstein distribution computedat a temperature T extracted from the energy density of the system and N norm ( k ) is the particlenumber normalized by the ratio N eq ( k ∗ , T ) /N ( k ∗ , T ) for one given value of k ∗ . The reason to use N norm ( k ) is that, for large ξ , the mode k ≈ T is outside the box which means that the estimatedtemperature is not accurate and so there would be an offset between the equilibrium distributionand the numerical result. Note that if the system thermalizes but decays are not efficient thedistribution would have a chemical potential ( µ ), since particle number is conserved at the time ofthermalization. We find that the chemical potential is small compared to the temperature and it14 .0 2.5 3.0 3.50.10.20.51 ξ ( N - N eq ) / N eq Total particle numberThreshold for thermalization3x10 - e ξ ξ P a r t i c l e N u m be r f / H ( N - N eq ) / N eq Total particle numberThreshold for thermalization f / H P a r t i c l e N u m be r Figure 2: Left plots: average normalized difference to a Bose-Einstein distribution (defined in eq. (4.1)).Right plots: Total particle numbers and threshold for thermalization. Upper plots: f = H and varying ξ . Bottom plots: ξ = 3 . f /H .is only visible in the smallest k mode for such cases (see appendix A.5). For this reason we discardthe first mode in the average difference. • Some of the plots include regions with f (cid:39) H , which are close to the UV cutoff. However, here weonly want to show that our estimates agree well with the numerical results in a region of parameterspace where the numerics works well, so that we can then extrapolate for larger values of f and ξ ,where the numerical treatment becomes more difficult to perform.On generic grounds, the numerical results confirm our expectations, i.e. , that thermalization occurswhen the particle number, controlled by ξ , is larger than a given threshold. If f /H decreases, thethreshold becomes lower and the system thermalizes more easily. The thresholds for thermalization arein agreement with eq. (3.19). We separate the results in 2 different cases of interest. First we consider asystem with γ + modes only. This is the case where thermalization is the most efficient because the gaugefield number is larger. Then, we include scatterings with γ − and φ and numerically evaluate the thresholdfor thermalization. In appendix A.5 we also add some plots which show the complete distributions as afunction of momenta, at a given time.Let us start by considering only γ + gauge fields and their scatterings, with matrix elements given ineq. (A.11). In fig. 2 we show the total particle number and the average difference to a thermal distribution∆ N/N , defined in eq. (4.1), for fixed f and varying ξ and vice-versa, evaluated at the end of the simulation15when the longest mode is horizon size). We assume the system has thermalized when ∆ N/N (cid:46) . ξ = 3 .
9, ∆
N/N decreases sharply to roughly 0 . f (cid:46) H , while for f = H that happens for ξ (cid:38) .
7. This is in agreement with the plots on the right which show that the analytical estimate for thethreshold of thermalization, using eq. (3.18), intersects the total particle number N total γ + = (cid:80) k N γ + ( k ) atroughly the same values of f and ξ . In the upper right plot we also see that the particle number dependsexponentially on ξ and is well approximated by N γ + (cid:39) × − e . ξ , as expected from the analyticalestimates in sec. 3.1.In the second case we consider all scatterings and decays involving φ, γ + and γ − . The presence ofdecays and of more degrees of freedom leads to a decrease in the γ + particle number which implies lessefficient thermalization, i.e , it happens for smaller (larger) value of f ( ξ ). In fig. 3 we show the sameplots as in fig. 2. By looking at the average differences to a Bose-Einstein distribution (left plots) one cansee that for ξ = 3 .
9, lower plots, the two gauge field polarizations thermalize at f (cid:46) H while for f = H all species thermalize when ξ (cid:38) .
9. Comparing the average differences, on the left, with the particlenumbers and estimates, on the right, we see that thermalization requires a value of ξ about 0.6 largerand a value of f /H about two times lower than the estimate.Such results seems to have parallel in the results obtained in [25] where the authors also found thatthe power spectrum of both polarizations of the gauge field and of the inflaton perturbations follow eachother 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 f such that ∆ N/N < .
1. The threshold is well fitted by ξ = 0 .
44 log (cid:18) fH (cid:19) + 3 . , (4.2)which is close to the estimated threshold in eq. (3.19), up to an offset of 0.7 in ξ . We will use thisthreshold equation in the next section when discussing the phenomenology. After thermalization is reached, however, the evolution is expected to be non-trivial. Indeed the tempera-ture 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 theinstability 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 dρ φ /dt = ˙ φ/ (4 f ) (cid:104) F ˜ F (cid:105) =˙ φ/f (cid:82) ( d k ) / (2 π ) kd/dτ ( | A + | − | A − | ) /a , averaged over a thermal distribution; such a quantity shouldbe suppressed if γ + and γ − are in equilibrium, since the scatterings tend to restore parity symmetry. Wedo not try to estimate such quantity exactly here because, while N γ − is easy to compute in equilibrium,the same is not true for N γ + , due to its complex frequency ω + . We postpone to future work a refined nu-merical 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 andtensor curvature perturbation, ζ and h respectively, should be affected by thermalization: • Superhorizon conservation:
On superhorizon scales the gauge field mode function goes to a constant, so its energy densitydecreases as a − . Because the sourcing of adiabatic curvature perturbation is proportional to this16 + γ - ϕ ξ ( N - N eq ) / N eq N total N γ + ( β s ( f / H ) ) / ξ P a r t i c l e N u m be r γ + γ - ϕ / H ( N - N eq ) / N eq N γ + ( β s , γ ( f / H ) ) / f / H P a r t i c l enu m be r Figure 3: Left plots: average normalized difference to a Bose-Einstein distribution (defined in eq. (4.1))for each of the degrees of freedom ( φ, γ + , γ − ). Right plots: Total particle numbers and threshold forthermalization. Upper plots: f = H and varying ξ . Bottom plots: ξ = 3 . f .17uantity, we expect a negligible isocurvature sourcing of the scalar curvature perturbation in theuniform density gauge, ζ . In appendix A.7 we elaborate on this point. Therefore, it should be agood estimate to evaluate correlators of ζ at horizon crossing. We also disregard here the possibilitythat a backreaction regime, which could change the predictions on perturbations, can be present atany time during inflation. • Loop corrections:
In the absence of thermalization strong constraints on ξ were derived from non-Gaussianities [6,8,12]and the requirement of perturbativity [9, 10] due to loop corrections involving the gauge fields. Weexpect these corrections to become smaller in the thermal regime by the fact that energy movesfrom the horizon size to smaller scales, possibly k/a (cid:39) T , and so at horizon crossing the effect oncorrelators should be smaller. • Parity symmetry:
On top of the above suppression, (cid:104) F ˜ F (cid:105) should also be suppressed in a thermal environment due thetendency of the scatterings to restore parity symmetry. Because each vertex introduces one powerof F ˜ F , loop corrections from the axial coupling to odd ζ correlators are proportional to (cid:104) F ˜ F (cid:105) andso they will also be suppressed. In this subsection we estimate the power spectrum of curvature perturbation assuming that a thermalregime is reached and under the assumptions stated above. We will have different cases depending onwhether perturbations are thermalized or not.Let us briefly review, first, the standard vacuum case. In this case the mode functions associated withthe canonically normalized field are given, in the de-Sitter limit, by | u k | = 12 k (cid:12)(cid:12)(cid:12)(cid:12) − ikτ (cid:12)(cid:12)(cid:12)(cid:12) (cid:39) − kτ → k τ , (5.1)which means that the power spectrum of the scalar curvature perturbation in the uniform density gauge, ζ , in the superhorizon limit, − kτ →
0, is [26] P vac ζ ≡ | ζ k | k π = (cid:12)(cid:12)(cid:12)(cid:12) Hk π a ˙ φ u k (cid:12)(cid:12)(cid:12)(cid:12) = H π ˙ φ . (5.2)Alternatively, one could have evaluated the subhorizon expression for u at horizon crossing, − kτ (cid:39) ζ , which is constant on superhorizon: | u k | (cid:39) − kτ → k ⇒ P vac ζ = H k π ˙ φ ka (cid:12)(cid:12)(cid:12)(cid:12) − kτ → = H ∗ π ˙ φ ∗ , (5.3)where ∗ denotes quantities evaluated at horizon crossing. This last result is actually more general because H ∗ 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 aboveequations, 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.18f the inflaton perturbations are instead thermal, the expected outcome is much simpler by followingthe same strategy of evaluating u k at horizon crossing. In this case u k is related to the particle number N k through eq. (3.2), which should be reliable slightly inside the horizon. Then, using | u (cid:48) k | (cid:39) k | u k | , wehave (cid:12)(cid:12) u therm k (cid:12)(cid:12) = 1 k (cid:18)
12 + N k (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) − kτ → (cid:39) k T ∗ H ∗ , (5.4)which then implies that the curvature perturbation is simply obtained by replacing the vacuum particlenumber (1/2) with the thermal value T ∗ /H ∗ : P therm ζ = T ∗ H ∗ H ∗ π ˙ φ ∗ . (5.5)Note that this prediction naturally has similarities with those of warm inflation [2, 27] from the fact thatboth use the thermal particle number. However we do not consider any friction term, induced by thethermal bath, in the equation of motion for ζ . Regarding the tilt of the spectrum, although both casesgive a nearly scale invariant power spectrum, in the thermal case the departure from scale invariance hasa different functional dependence on the slow-roll parameters. In the vacuum case the spectral index isdue to the time variation of H and ˙ φ at horizon crossing, while here the ratio T /H at horizon crossingalso matters.Indeed the spectral index is given by: n s − ≡ d ln P therm ζ d ln k = − (cid:15) H + 2 η + d ln( T ∗ /H ∗ ) d ln k , (5.6)where we have used the following definitions for slow-roll parameters: (cid:15) H ≡ − ˙ HH , η ≡ (cid:15) H − ¨ φH ˙ φ . (5.7)The reason why we introduced (cid:15) H is to keep the analysis fully general: in the backreaction dominatedcase , in fact, this parameter differs from the previous definition, (cid:15) = ˙ φ / (2 H M P ).We can also compute the tensor-to-scalar ratio. Tensor perturbations ( h ) have couplings to gaugefields √ (cid:15) smaller than the scalars and so are more difficult to thermalize. Therefore, while we assumehere ζ to be thermal, we analyze first the most likely case of a standard vacuum spectrum , for tensors.In this case the tensor to scalar ratio is suppressed, because scalar perturbations are enhanced, giving: r ≡ P h P therm ζ = 16 (cid:15) H ∗ T ∗ . (5.8)This case is phenomenologically very interesting because it would help polynomial large field models( V ( φ ) ∝ φ n , with n > (cid:15) and η are comparable, to agree withobservations. In fact the observational bound r (cid:46) . (cid:15) , while inour case if T ∗ /H ∗ is sufficiently large the bound is relaxed and those models are not anymore in tension Note, however, that in the backreaction case one might also worry that ζ could have a non-negligible contribution fromgauge field fluctuations at horizon crossing, which is not taken into account by the above eq. (5.5). We neglect here possible 1-loop contributions to P h . In principle, in fact, there could also be a region where the tensorsare non-thermal but have relevant 1-loop corrections to their 2-point function. i.e. , enhanced by 2 T ∗ /H ∗ over thevacuum case P thermal h = H ∗ T ∗ π M p . (5.9)In this case the tensor-to-scalar ratio would remain unchanged, r = 16 (cid:15) , while the tensor tilt would benon-trivial: n T ≡ d ln P therm h d ln k = − (cid:15) H + d ln( T ∗ /H ∗ ) d ln k . (5.10)This regime would actually be very interesting. In fact, normally the tensor tilt is necessarily red ( n T < (cid:15) H > T ∗ /H ∗ . It is however difficult to havethermalization of tensors, due to their weaker coupling, as we will see in section 6. f and ξ We start by briefly reviewing the main observational constraints on ξ derived in the literature and thendiscuss 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 ξ (cid:46) . H f e πξ π l < , (5.11)where l is some loop factor. For l (cid:39) and imposing P vac ζ = 2 . × − , the previous bound constrains ξ (cid:46) . ξ (cid:46) . V (cid:48) ( φ ) (cid:39) H ˙ φ (cid:28) (cid:104) F ˜ F (cid:105) f . (5.12)In the case without collisions we have (cid:104) F ˜ F (cid:105) (cid:39) − H e πξ /ξ [5], which means that backreaction isabsent if f /H (cid:29) × − e πξ /ξ / .However, we stress here that these constraints do not directly apply if the system is thermalizedsimply because the occupation numbers completely change and new constraints should be derived. Infact, as we mentioned before, since thermalization moves particles from the horizon size to smaller scales,we expect the thermal case to be less constrained, as we discuss here.But there are some simple constraints that we can directly apply. Namely, irrespectively of whetherthe spectrum of perturbations is thermal or not, the vacuum spectrum cannot exceed the observationalbound, P vac ζ ≤ P obs ζ ≡ . × − . The reason is that vacuum fluctuations are a lower bound and adiabaticperturbations cannot be erased once they cross the horizon [28]. If we now recall the definition of ξ , this20 �� � � � � � � �� �� � � � � � � � � �� �� � ϕ - � � � � � � �� �� �� �� �� - �������������� � �� - � � � � � � � � �� � �� � ( � � � � � � � � � � � � �� �� �� �� ) � �� � � ��� �� �� ( � � � � � � � � � � � � �� �� �� �� ) ξ f / H Figure 4: Green: gauge fields thermalize through SM interactions. Red: Region where the perturbativeexpansion breaks down in the absence of thermalization (straight line derived from a parametric estima-tion and dashed from an explicit 1-loop calculation [9, 10]). Purple: backreaction region in the absenceof thermalization. Orange: inflaton perturbations thermalize (and also gauge fields thermalize, even inabsence of SM interactions). Blue: Vacuum fluctuations larger (or equal) than the observed value.imposes a lower bound on f ,4 πf ξH = 1 (cid:113) P vac ζ (cid:38) × ⇒ fH (cid:38) × ξ . (5.13)In fig. 4 we plot the different thresholds for thermalization derived in sections 3.3 and 4. Jointlywith the thermalized regions we overlap the perturbative constraint eq. (5.11) (derived in absence ofthermalization) for l = 10 , the backreaction threshold (also in absence of thermalization) eq. (5.12) andthe observational constraint eq. (5.13). In absence of SM interactions the phenomenological allowed regionfor φ to be thermalized is such that one would need to take into account higher order loop correctionsto the propagators as well as backreaction, at least before the onset of thermalization. Therefore, adedicated study is required. However, we stress that in the presence of SM interactions the system canthermalize before backreaction although, in this case, it is also crucial to take into account thermal massesto establish when φ thermalizes and when backreaction is relevant. We address this issue in the nextsection.Finally, we argue why, generically, the constraints on ξ are very much alleviated when the systemthermalizes. The reason is two-folded. First, the fact that energy moves from the horizon scale tothe UV, with a consequent decrease in the particle number, changes the peak of the integrals fromthe horizon scale to the temperature scale. This will affect all correlators. Moreover, odd correlatorsare further suppressed by the following reason. The interaction Hamiltonian associated with the axial21oupling is given by [6, 8] H int = ξ (cid:90) d x ζ F µν ˜ F µν . (5.14)Therefore, all odd correlators will necessarily be proportional to, at least, one power of (cid:68) F ˜ F (cid:69) ∝ d/dτ ( N + − N − ). In a thermal bath N + and N − tend to equilibrate through scatterings and so anyparity asymmetry tends to be suppressed. This asymmetry would be completely erased if the two he-licities had exactly the same dispersion relation. However this is not the case, since the axial couplingchanges the dispersion relations of ω + and ω − in different ways, and so a small asymmetry should remainthere, proportional to ξH , even in thermal equilibrium.Computing loop corrections in the thermal regime requires a dedicated study [29]. Nonetheless, inappendix A.6 we provide a rough derivation of the parametric dependence of the non-Gaussian parameter f NL (cid:39) (cid:10) ζ (cid:11) / (cid:10) ζ (cid:11) with the temperature, considering the case where gauge fields are thermalized but φ is not. We find f NL (cid:39) c ξ P vac ζ O (cid:18) T H (cid:19) , (5.15)where c is a small number containing inverse powers of (2 π ) and the result of the angular integration andso needs to be derived in a more accurate computation. If we consider the case of instantaneous thermal-ization of the plasma, then, the energy density in the gauge fields is simply given by the initial energy ρ γ (cid:39) − H e πξ /ξ [21]. Therefore, the plasma would have a temperature ¯ T (cid:39) ρ / (cid:39) . He πξ/ /ξ / which means that f NL would be proportional to e πξ . Comparing this result with the non-thermalcase where f NL (cid:39) − e πξ P vac ζ /ξ [6] there is a parametrical suppression of e − πξ which translatesinto a weaker constraint for ξ . In particular, assuming 10 − > c > − and using the constraint f NL < O (10) [30] one would find ξ (cid:46) −
7. However such a high temperature ¯ T is actually not whatwe finally expect, as we discuss in the next section. In reality, we expect a much smaller temperatureto be reached, thus relaxing even more the constraints. In appendix A.6 we also argue that, if φ is alsothermalized, f NL can be at most multiplied by an additional power of T /Hf thermal NL (cid:46) d ξ P vac ζ O (cid:18) T H (cid:19) , (5.16)where d is another small coefficient. So far we have studied the conditions for thermalization to occur due to the initially large occupationnumber of gauge bosons. However, it is more subtle to understand what happens after thermalization isestablished.As we discussed in previous sections, due to the presence of an imaginary dispersion relation for γ + ,it is unclear how is the shape of the distribution at low momenta, k/a < ξH . In fact such modes areexponentially produced, although in a finite window of time, and so one expects relevant deviations froma BE distribution such as the presence of a peak at low momentum. The γ − polarization, instead, shouldbe correctly described by a BE distribution with a dispersion relation ω = k + m , where m = − kξ/τ ,22hile δφ , if thermalized, is described by a massless BE distribution.However a very important ingredient modifies the above discussion: when thermalization happens,gauge fields typically develop a thermal mass m T ∝ gT , where g is the gauge coupling, as a resultof having a non zero density of particles correcting the propagator at 1-loop. Note, importantly, that φ instead does not get a thermal mass, since its interactions only renormalize its kinetic term [29]. Weleave to the future a more careful study of the effect of thermal masses, but we anticipate here someinteresting features.Let us parameterize the thermal mass as m T ≡ ¯ gT . (6.1)In the case of SM fermions ¯ g = g (cid:80) i q i /
6, where q i is the charge of a given fermion. At energies ofaround 10 GeV this gives ¯ g (cid:39) .
3. If the thermal mass is included in the equation of motion for thegauge fields it will compete with the instability due to the axial coupling, giving A (cid:48)(cid:48)± + ω T ( k ) A ± = 0 , ω T ( k ) = (cid:18) k ± kξτ + m T H τ (cid:19) . (6.2)When m T ≥ ξH the mass completely shields the instability band. Therefore, gauge fields cannot beproduced and we might expect an equilibrium temperature T eq = ξH ¯ g . (6.3)In fact, starting from a non-thermal case, when we reach thermalization the plasma has an initial temper-ature, ¯ T , typically higher than T eq . Thus, at that point the particle production stops and the plasma ofparticles is expected simply to redshift as radiation, as a − , which is equivalent to say that temperaturedecreases as 1 /a . When the temperature drops below T eq an instability band reappears, most likely asa narrow band centered at | kτ | ≈ ξH , which is the minimum of ω T ( k ) . At this point we expect thesystem either to reach an almost stationary state slightly below T eq or perhaps an oscillatory behavioraround such temperature.When the temperature drops, thermalization could be less efficient. We can estimate the condition forthe plasma to remain in thermodynamic equilibrium by comparing σ eq · n eq with H , using the equilibriumnumber density n eq (cid:39) T and the typical cross section σ eq (cid:39) ¯ α /T , where ¯ α ≡ g / (4 π ) (cid:80) i q i issummed over different species. In the SM case at 10 GeV, ¯ α (cid:39) .
05. As a result we get: TH (cid:29) α (cid:39) . (6.4)To estimate if the plasma keeps thermalized down to T eq we should fulfill the condition T eq H (cid:29) α ⇒ ξ (cid:29) ¯ g ¯ α (cid:39) . (6.5)For a given g and a given set of fermions, this condition can indeed be met for sufficiently large ξ whilewhen considering only SM fermions we get a more precise condition on ξ . Note however that, even if the This certainly happens in the well-known case of having interactions with charged SM particles. It could also happenfor the case of only φ -mediated interactions, but this would require performing the one-loop thermal correction, which wepostpone to a subsequent publication [29]. T > T > T eq we still expect the plasma simply to redshift down to T eq ,because the instability remains screened by the thermal mass. Equilibrium temperature
In the rest of this section we will discuss the case where the system reaches a stationary temperature T eq . In that case we expect a narrow band of the instability around | kτ | ≈ ξH to be constantly present,generating a peak in the γ + distribution, whereas we expect a thermal distribution away from suchnarrow peak, where ω T ( k ) >
0. Observing this dynamics should be numerically feasible but we postponefurther investigation of this important issue to future work. We now briefly discuss possible interestingobservational consequences of such cases.We have to distinguish between the case with SM interactions and the case with only φ -mediatedinteractions. In the first case, as we saw before, we can reach a thermal state if ξ (cid:38) .
9. Then we needto check if φ thermalization is efficient by using, as above, the thermal condition σ eq · n eq (cid:29) H , with σ eq ≈ T /f , so that: T (cid:29) ( f H ) / , (6.6)which should be true at T = T eq , thus imposing ξ (cid:29) ¯ g (cid:18) fH (cid:19) / . (6.7)If the above condition is met we can use the results from the previous section to make predictions for thespectrum of perturbations.If φ is thermal from eq. (5.5) we have that: P therm ζ = ξ ¯ g H ∗ π ˙ φ ∗ = H π f ¯ gξ , (6.8)and the spectral index becomes simply n s − ≡ d ln P therm ζ d ln k = − (cid:15) H + 2 η + ˙ ξHξ = − (cid:15) H + η . (6.9)We remind the reader that we are still working under the assumptions listed at the beginning of section 5,which are valid in the non-backreacting regime and imply that ζ is conserved on superhorizon scales.In the case of tensor modes the couplings are gravitational, suppressed by 1 /M p , and so the crosssection would be given by σ eq ≈ T /M p and so in order to reach thermalization one would need T (cid:29) ( HM p ) / . This requires a temperature larger than the inflationary energy, ρ / inf = 3 H M p , and so,for that reason, it is not possible to thermalize the tensor modes during inflation at this equilibriumtemperature. Therefore, assuming tensor modes to be in the vacuum the tensor to scalar ratio would begiven by r ≡ P h P therm ζ = 16 (cid:15) H T ≈ (cid:15) ¯ gξ . (6.10)As we mention before, thermalization of φ leads to a suppression of the tensor to scalar ratio which, at24he equilibrium temperature, would be ¯ g/ (2 ξ ). If we look at fig. 5, this amounts to at least an O (10 − )suppression.We can also try to estimate whether there is relevant backreaction on the scalar field equation ofmotion, as follows. The energy extraction from φ is given by (cid:104) F ˜ F (cid:105) ˙ φ/f which we can estimate by usingconservation of energy density in the gauge fields ρ γ :˙ ρ γ + 4 Hρ γ = ˙ φ f (cid:104) F ˜ F (cid:105) . (6.11)In a stationary situation ˙ ρ γ ≈ (cid:104) F ˜ F (cid:105) / (4 f ) (cid:38) H ˙ φ . These two conditionsthen require ξ (cid:38) g fH . (6.12)Finally, using the estimation of non-Gaussianity in the case where φ is thermal, eq. (A.26), we can alsofind what would be a conservative viable region for T eq = ξ/ ¯ g by requiring the non-Gaussian parameter f NL < O (10). This would impose ξ (cid:46) −
50 which corresponds to a small viable window in fig. (5).In fig. 5 we plot the different constraints and phenomenological windows as in fig. 4, but now includingthe effects of a thermal mass and assuming a stationary regime. The different regimes are described byeqs. (6.5), (6.7), (6.12), (6.8) and (5.13). Note that, contrary to fig. 4, the region where φ is thermalizedis not completely inside the backreaction region and so there is some parameter space where φ canthermalize without backreaction.In the case of oscillatory temperatures around T eq the above description remains probably roughlycorrect although in that case we should expect larger deviations from the thermal predictions as wellas superimposed oscillations in the relevant quantities at horizon crossing. Clearly, if such oscillationsare large the model is observationally ruled out. However, interestingly, if oscillations are within theobservational bound but non-negligible, this would result in oscillatory power spectra, which could beanother striking evidence of this mechanism. In this work we have analyzed the system composed of an inflaton ( φ ) coupled to a gauge field ( A µ )through a CP-odd term with strength 1 /f . One helicity of the gauge field, γ + , is known to develop aninstability with exponential particle production at a rate controlled by the parameter ξ = ˙ φ/ (2 f H ), formodes with momenta in the range (8 ξ ) − H (cid:46) k/a (cid:46) ξH . We have included in this scenario the effectof scatterings through a set of Boltzmann-like equations. The scattering rates are Bose-enhanced by thelarge occupation numbers of γ + and can drive the system to thermal equilibrium, thus leading to a setupthat we dubbed thermalized axion inflation . We have analyzed two situations: (1) only scalar-gauge fieldinteractions, as in the case of a “dark photon”; (2) identifying the gauge field with a Standard Modelgauge boson (abelian or non-abelian). We stress that (2) is probably the best motivated case, since thissystem anyway needs to be coupled to the SM to provide reheating, and it turns out to be the one inwhich we have good control of the theory, since thermalization happens before backreaction and beforeperturbativity constraints.In the first case the most relevant scatterings are γ + γ + → γ + γ + . By solving the associated Boltzmann-25 �� � � � � � � �� �� � � � � � � � � �� �� � � �� � � ��� �� � � ( � � � � � � � � � � � �� �� �� � � ) �� - �������������� ϕ - �������������� � � � � � � � � � � ��� �� � � ζ = � � ��� - � ξ f / H Figure 5: Same conditions as fig. (4) but in the presence of a thermal mass and assuming an equilibriumtemperature of T eq = ξ/ ¯ g . In this plot we fix ¯ g = 0 .
5. We also show as a dashed line the condition thata thermal spectrum fits the observational value for P ζ .like set of equations of the system we find that N γ + , N γ − and N φ evolve to a Bose-Einstein distributionwhen ξ (cid:38) .
44 ln( f /H ) + 3 .
4. Note, however, that due to the axial coupling the gauge fields have amodified dispersion relation, imaginary in the case of γ + , which is expected to distort the distribution forlow momenta k/a < ξH . In a phenomenologically realistic model one has to require the power spectrumof curvature perturbations to be compatible with observations, and that imposes a lower bound on f /H ,which then implies thermalization at ξ (cid:38)
6. At such values, however, we would have to assume thatthe occupation numbers are not significantly altered by loop effects [9], which are important already at ξ (cid:38) .
5. Moreover, for such a large ξ the gauge fields would also backreact on φ and on δφ , so a dedicatedstudy is needed. We will return to this point in a future work.In the second case, where we consider SM interactions with the gauge field, thermalization is easierto achieve because the cross sections are not suppressed by powers of H/f . This case is also morepredictive because ξ is the only free parameter. For example, if the gauge field is the one associatedto U (1) hypercharge and considering only particle anti-particle productions, we find that thermalizationhappens already at ξ (cid:38) .
9. A similar condition on ξ applies if instead we consider the gauge field to bea gluon and include self-interactions. This is of great interest because in this case one can have a thermalbath of particles well under control, before perturbativity constraints become relevant and also beforebackreaction of gauge fields on φ becomes important.In both cases thermalization has profound implications for the phenomenology of axion inflation. Toelaborate on this we have further pointed out that, after thermalization is reached, the system shouldevolve in a very different regime, due to the presence of thermal masses m T = ¯ gT for the gauge bosons,which tends to screen the instability of γ + . As a result we expect either an oscillatory behavior or astationary solution at a temperature given by T eq = ξH/ ¯ g . Note that this regime is now only linearlydependent on ξ and not exponential anymore and so all constraints on ξ should become much weaker26nd should be properly readdressed. We described such behavior only qualitatively in this work but weplan to address these points in full detail in the future.Assuming the system reaches a stationary solution, or one where the temperature changes adiabati-cally, we derived very interesting consequences for cosmological observables such as the scalar and tensorperturbations, under the assumption that φ is thermalized, which can happen in the region shown infig. 5. In this regime one can also find a region of parameters where φ can thermalize without backreact-ing on the background. Note also that in this case the power spectrum of curvature perturbations can beobtained simply by setting its value at horizon crossing to be that of a thermal state, assuming it will beconserved on superhorizon scales. This leads to a result for P ζ that is enhanced compared to the vacuumcase by the ratio of particle numbers between the thermal state, T ∗ /H ∗ , and the vacuum, 1 /
2. We havealso computed the tilt of the power spectrum to be n s − (cid:15) H − η , which differs from the standardslow-roll formula because of the time dependence of T eq ∗ /H ∗ .Regarding the tensor modes, because of their Planck suppressed couplings it does not seem possiblefor them to thermalize, at least at the equilibrium temperature. For that reason we expect the tensors toremain in the vacuum and so, in the case in which φ is thermalized, a tensor to scalar ratio suppressedby H ∗ / (2 T ∗ ) with respect to the standard formula. This is very interesting because it can allow fora reconciliation between models of inflation where (cid:15) and η are comparable, e.g. polynomial large fieldmodels, with observations.Finally, although we did not present here a detailed study of non-Gaussianities, we provided strongarguments, and one particular estimation, which justify why the previous constraints do not apply to ourcase and why, generically, we expect them to be less restrictive in the thermal regime. Thermalizationredistributes in fact the occupation numbers: it depletes the highly populated horizon-sized modes andit populates higher k modes. As we have argued, loop corrections to cosmological correlators will then begenerically smaller than in the non-thermal case. In addition, odd correlators of ζ are further suppressed,since thermodynamic equilibrium tends to restore parity, driving (cid:104) F ˜ F (cid:105) to small values. We estimated thenon-Gaussian parameter f NL to be proportional to P vac ζ T /H , which corresponds to a large parametricsuppression of e − πξ compared to the non-thermal case. This seems promising since it could allow forinteresting regimes with large ξ , such as the φ -thermalized regime and the backreacting regime [29], tobe in agreement with observations.To sum up we provided a working model in which during inflation a thermal bath with a possibly largetemperature is present, leading to new implications in model building and to new observational features.One of the crucial ingredients of this setup is the fact that the axial coupling respects a continuous shiftsymmetry and so cannot induce thermal mass corrections. This is true in the U(1) case, while in non-Abelian theories it can be broken by non-perturbative effects. If the thermalized backreacting regime canbe achieved in full control it could even remove the standard need of an inflaton flat potential, by thehelp of the dissipative friction. Moreover reheating is already incorporated in this scenario and obtainedjust as a transition when the thermal bath starts dominating over the potential energy of the inflaton. Acknowledgments:
We thank K. Tywoniuk, J. Garriga, C. Germani, F. Mescia, E. Verdaguer, K. To-bioka and M. Sloth for comment and useful discussions. This work is supported by the grants ECFPA2010-20807-C02-02, AGAUR 2009-SGR-168, ERC Starting Grant HoloLHC-306605 and by the Span-ish MINECO under MDM-2014-0369 of ICCUB (Unidad de Excelencia Maria de Maeztu). This estimation holds if φ is not thermalized, while in the opposite case we showed that there can be, at most, oneextra power of T/H . Appendix
A.1 Fourier conventions
We quantize scalars ( δφ or ζ ) and gauge fields as ζ ( x, τ ) = (cid:90) d k (2 π ) e − ikτ (cid:104) ζ k a k + ζ ∗ k a †− k (cid:105) ,(cid:126)A ( x, t ) = (cid:90) d k (2 π ) e − ikτ (cid:88) σ (cid:126)e k,σ (cid:104) A k,σ b k,σ + A ∗ k,σ b †− k,σ (cid:105) , (A.1)where we have used the Coulomb gauge ( A = 0) and where the creation and annihilation operatorssatisfy (cid:104) a k , a †− q (cid:105) = (2 π ) δ (3) ( k + q ) , (A.2) (cid:104) b k,σ , b †− q,σ (cid:48) (cid:105) = (2 π ) δ σ,σ (cid:48) δ (3) ( k + q ) . (A.3)The polarization vectors satisfy the following identities (cid:126)k × (cid:126)e ± ( (cid:126)k ) = ∓ (cid:126)k(cid:126)e ± ( (cid:126)k ) , (cid:126)e σ ( (cid:126)k ) · (cid:126)e σ (cid:48) ( (cid:126)k ) ∗ = δ σσ (cid:48) , e ± ( (cid:126)k ) = e ∓ ( − (cid:126)k ) = e ∓ ( (cid:126)k ) ∗ . (A.4) A.2 Decay estimation from 1-loop correction
In section 3.2 we presented the 1-loop correction to the two point function which we then used as anestimation of the decay rate. In this appendix we present some of the intermediate steps which took usto eq. (3.14). Starting from eq. (3.13) by taking the time derivative of the two point function we arriveat ddτ (cid:104) u k u k (cid:105) loop ( τ ) = 2 Ha (cid:104) δφ k δφ k (cid:105) loop ( τ ) + a (cid:90) τ −∞ dτ (cid:48)(cid:48) (2 π ) (2 f ) (cid:90) d q (cid:12)(cid:12)(cid:12) (cid:126)e + ( (cid:126)q ) · (cid:126)e − ( | (cid:126)k − (cid:126)q | ) (cid:12)(cid:12)(cid:12) ×× (cid:68)(cid:104) A (cid:48) q ( τ ) A | (cid:126)k − (cid:126)q | ( τ ) | (cid:126)k − (cid:126)q | δφ k ( τ ) , (cid:104) A (cid:48)| (cid:126)k − (cid:126)q | ( τ (cid:48)(cid:48) ) A q ( τ (cid:48)(cid:48) ) qδφ k ( τ (cid:48)(cid:48) ) , δφ k ( τ ) δφ k ( τ ) (cid:105)(cid:105)(cid:69) ++ perm. (A.5)The time integrations are dominated by the latest time simply because, before thermalization, that iswhen the gauge field particle number has its largest value. This is only true until horizon crossing,afterwards the integral goes quickly to zero. Therefore, we simplify the previous expression by replacingthe time integration by its integrand evaluated at τ times the integration interval which we approximateto be ∆ τ = 2 ξ/ min( q, | (cid:126)k − (cid:126)q | ) corresponding to the time at which the largest mode enters the resonantband and so where the process becomes non-zero. Thus, we arrive at ddτ (cid:104) u k u k (cid:105) loop ( τ ) = 2 Ha (cid:104) δφ k δφ k (cid:105) loop ( τ ) + 2 ξa (2 π ) (2 f ) (cid:90) d q min( q, | (cid:126)k − (cid:126)q | ) (cid:12)(cid:12)(cid:12) (cid:126)e + ( q ) · (cid:126)e − ( | (cid:126)k − (cid:126)q | ) (cid:12)(cid:12)(cid:12) ×× (cid:68)(cid:104) A (cid:48) q ( τ ) A | (cid:126)k − (cid:126)q | ( τ ) | (cid:126)k − (cid:126)q | δφ k ( τ ) , (cid:104) A (cid:48)| (cid:126)k − (cid:126)q | ( τ ) A q ( τ ) qδφ k ( τ ) , δφ k ( τ ) δφ k ( τ ) (cid:105)(cid:105)(cid:69) ++ perm. (A.6)28ur goal is to rewrite this expression as powers of 2-point functions which we then associate with particlenumbers. However, by isotropy the particle number should not depend on the angle, so it is convenientto perform the angular integral. At this level we do not have access to the generic mode functions inthe presence of scatterings, so we parameterize this unknown by b which was found in similar 1-loopcomputations to be O (10 − ). The drawback is, however, that we now need to specify the angle between k and q at which the angular integral would peak in order to compute | (cid:126)k − (cid:126)q | . We assume | (cid:126)k − (cid:126)q | (cid:39) k + q .Under these simplifications we arrive at ddτ (cid:104) u k u k (cid:105) loop ( τ ) = 2 Ha (cid:104) δφ k δφ k (cid:105) loop + 2 ξa bδ (0)(2 π ) (2 f ) (cid:90) dqq | (cid:126)k − (cid:126)q | min( q, | (cid:126)k − (cid:126)q | ) ×× (cid:68)(cid:104) A (cid:48) q A k + q | (cid:126)k − (cid:126)q | δφ k , (cid:2) A (cid:48) k + q A q qδφ k , δφ k δφ k (cid:3)(cid:105)(cid:69) + perm. , (A.7)where we omitted the time argument in the right hand side because all quantities are evaluated at τ . Bya similar reasoning the first term in the right hand side of the previous equation gives4 Hbξ a δ (0)(2 π ) (2 f ) (cid:90) dq q | (cid:126)k − (cid:126)q | min( q, | (cid:126)k − (cid:126)q | ) (cid:68)(cid:104) A (cid:48) q A k + q | (cid:126)k − (cid:126)q | δφ k , (cid:2) A (cid:48) k + q A q qδφ k , δφ k δφ k (cid:3)(cid:105)(cid:69) . (A.8)We only deal with subhorizon modes, thus, the first term is smaller or at most equal to the second termso we neglect it. Now, in order to write the expression as a particle number we approximate | A (cid:48) k | (cid:39) kA k and make use of the approximate relation 1 / N u ( k ) ≈ k | u k | . The commutator structure implies theexpectation value to give an expression of the formIm [ (cid:104) δφ k δφ k (cid:105) ] Im [ (cid:104) δφ k δφ k (cid:105) (cid:104) A q A q (cid:105) (cid:104) A p A p (cid:105) ] . (A.9)The imaginary part of the two point function can be roughly seen as the vacuum contribution, (2 π ) , whilethe real part would be proportional to (2 π ) (2 N ). This identification is accurate for thermal propagatorsand in fact by making this identity we would find that the previous equation has indeed a structure ofparticle numbers similar to those appearing in a standard decay expression which would be proportionalto N u ( k )(1 + N γ + ( q ))(1 + N γ + ( k + q )) − N γ + ( q ) N γ + ( k + q )(1 + N u ( k )) . However, the 1-loop correction contains more processes than just the decay and for that reason we donot expect to arrive at that precise form. Nevertheless we can isolate the terms which do have the formof a decay and so we finally end up with dN u ( k ) dτ = − ξ × b (2 π ) (2 f ) a k ×× (cid:90) dq q ( k + q )min( q, k + q ) N u ( k )(1 + N γ + ( q ))(1 + N γ + ( k + q )) − N γ + ( q ) N γ + ( k + q )(1 + N u ( k )) . (A.10)The extra factor of 4 comes from identifying Re[ (cid:104) XX (cid:105) ] ∝ N X .29 .3 Matrix elements for scatterings In this appendix we provide the matrix elements for the scatterings depicted in fig. 1. We start by the γ + γ + → γ + γ + scattering which only has non-zero contributions from the s-channel M γ + γ + γ + γ + = 14 a f (cid:15) abcd k a k b e c ( k ) e d ( k ) (cid:15) wxyz k w k x e ∗ y ( k ) e ∗ z ( k )2 | k µ − k µ | + 4 perm. , (A.11)where (cid:15) abcd is the 4-dimensional Levi-Civita tensor and k µ are the 4-momenta of the particles involvedin the scatterings. The processes involving γ + and γ − can be readily obtained from the previous resultby making use of the identities in eq. (A.4) and summing over all the channels u, s, t . On the other handthe gauge field-axion scattering, γ + φ → γ + φ , has instead s and t-channel contributions and is given by M γ + φγ + φ = 14 a f (cid:15) abcd k a k b e c ( k ) (cid:15) wxyz k w k x e ∗ y ( k ) η dz | k µ − k µ | + k ↔ k , (A.12)where η dz is the Minkowski metric. Finally gauge field-axion conversion, γ + γ + → φφ , has matrix elements M γ + γ + φφ = 14 a f (cid:15) abcd k a k b e c ( k ) (cid:15) wxyz k w k x e ∗ y ( k ) η dz | k µ − k µ | + 2 perm. (A.13) A.4 Particle number and Bogolyubov coefficients
In curved space the particle number is typically defined as N k = | β k | , where [31] β k = f ( k ) f (cid:48) ( k ) − f (cid:48) ( k ) f ( k ) W k , (A.14)is the Bogolyubov coefficient which relates the mode functions between two different basis of states and W k = ( f (cid:48) ( k ) f ∗ ( k ) − f (cid:48)∗ ( k ) f ( k )) / (2 i ) is the Wronskian. If we want to know the particle number comparedto the vacuum state we fix f = e − ikτ / √ k and so N k = | β k | = | f (cid:48) ( k ) | + k | f ( k ) | k , (A.15)in agreement with our definition in eq. (3.2). A.5 Numerical distribution of particle number
In fig. 6 we show the numerical distribution of particle number for the different species as a function ofmomenta and at the end of the simulation, − k min τ = 1, for ξ = 3 . f = H . In the same plot weshow Bose-Einstein distributions fro massless particles with and without a chemical potential. We cansee that, apart from the longest mode, a Bose-Einstein distribution with T = 45 H is in rough agreementwith the distribution. The inclusion of a chemical potential µ = 10 − . T gives a perfect fit, includingalso the first mode. A.6 Estimation of loop corrections
In the main text we have argued that in the thermal regime all constraints, such as non-Gaussianities, non-perturbative constraints and backreaction, are strongly modified by the fact that the energy moves from30 γ - N γ + N ϕ Bose - Einstein ( T =
45 H ) Bose - Einstein ( μ = - T, T =
45 H ) - k η f Figure 6: Particle number as a function of momenta at a time τ f for γ + , γ − and φ . The dashed linesare fits using thermal distributions with and without chemical potentialthe horizon scale to the UV. Moreover, odd correlators are also further suppressed by parity arguments.In this appendix we give a rough estimation of the 1-loop correction to the 3 point function. We study thecase where the photons are thermalized but the inflaton perturbations are not. We also derive an upperbound for the case with thermalized inflaton perturbations. This estimation is also easily generalizableto higher order correlators. Let us start by the 3-point function. The 1-loop correction is proportional to (cid:104) ζ k ζ q ζ p (cid:105) ( τ ) = i (cid:90) τ −∞ dτ (cid:90) τ −∞ dτ (cid:90) τ −∞ dτ (cid:104) [ H int ( τ ) [ H int ( τ ) , [ H int ( τ ) , ζ k ( τ ) ζ q ( τ ) ζ p ( τ )]]] (cid:105) , (A.16)where H int is the interaction Hamiltonian given by [6, 8] H int = ξ (cid:90) d x ζ F µν ˜ F µν . (A.17)We are interested in the parametric dependence of the correlators with the temperature in order tocompare it to previous results. For that reason we will ignore factors of 2 π and other O (1) numbers.We will assume the result to be dominated by the thermal part of the gauge field propagators. We alsoneglect the commutator structure, which should be an overestimation, since there could be cancelationsdue to its structure.Let us first analyze the vertex structure. Generically, the 1-loop correction to the n-point functionwill be proportional to (cid:68)(cid:16) F µν ˜ F µν (cid:17) n (cid:69) . (A.18)As we mentioned in the text, in the thermal regime we expect (cid:68) F ˜ F (cid:69) to be suppressed because paritytends to be restored. In that case, the correlators involving the gauge fields will be dominated by (cid:104) F F (cid:105) F µν ˜ F µν = F µν F µν we estimate (cid:68)(cid:16) F µν ˜ F µν (cid:17) n (cid:69) (cid:39) (cid:104) F F (cid:105) n , if n is even (cid:68) F ˜ F (cid:69) (cid:104) F F (cid:105) n − , if n is odd (A.19)up to permutation factors which we neglected here. The term in (cid:104) F F (cid:105) is dimensionally an energy densityand, using eq. (3.2), we approximate in the thermal regime (cid:104)
F F (cid:105) (cid:39) kN ( ω + ( k )), barring cancellations.On the other hand, the term in (cid:68) F ˜ F (cid:69) is proportional to the parity asymmetry (cid:68) F µν ( k ) ˜ F µν ( k ) (cid:69) = kd/dτ ( | A + | − | A − | ) (cid:39) d/dτ ( N ( ω + ( k )) − N ( ω − ( k )) where ω + , − ( k ) = ( k ± kξ/τ ) / . Assuming athermal distribution we find ,in limit of large temperature, and for − kτ (cid:29) ξddτ [ N ( ω + ( k )) − N ( ω − ( k ))] (cid:39) ξTHk τ , (A.20)Note, however, that in the in-in computation the vertices are evaluated at different times and so theycould connect different physical momenta. Nevertheless, because the integrals are regulated in the UVby the particle number and in the IR by the momentum dependence of the integrand, we expect theintegrals to peak when the gauge fields running in the loop have physical momenta of the order of thetemperature s phys (cid:39) T .Under these approximations and after Fourier transforming all the interaction Hamiltonians andintegrating over the delta functions, the 1-loop correction simplifies to (cid:104) ζ k ζ q ζ p (cid:105) ( τ ) (cid:39) − ξ (cid:16) P vac ζ (cid:17) k δ (3) (cid:16) (cid:126)k + (cid:126)q + (cid:126)p (cid:17) (cid:90) d s (cid:90) τ −∞ dτ τ sN ( ω + ( s )) | τ (cid:90) τ −∞ dτ τ sN ( ω + ( s )) | τ (cid:90) τ −∞ dτ τ ddτ [ N ( ω + ( s )) − N ( ω − ( s ))] | τ , (A.21)where we assumed the integrals to peak in the equilateral configuration, k (cid:39) q (cid:39) p , similarly to previousresults in the non-thermal case [6], although there might be some changes in the thermal regime. Inthe last expression we have also replaced ζ ’s in the time integrals by their subhorizon expression ζ k = H τ / ( ˙ φ √ k ), while the external ones by ζ k = H / ( ˙ φ √ k ), because they are evaluated at horizon crossing.Finally, to keep the treatment simple we further assume the integral to peak at a scale s (cid:29) k, q, p . Wecould also include corrections with s (cid:39) k, q, p but we expect not to have higher powers of the temperatureso we neglect them for this estimation.The time integrals should peak when the s/a = − sHτ = T . At that scale the particle number is N + , − ( T ) (cid:39) (cid:104) ζ k ζ q ζ p (cid:105) ( τ ) = − ξ (cid:16) P vac ζ (cid:17) k δ (3) (cid:16) (cid:126)k + (cid:126)q + (cid:126)p (cid:17) (cid:90) d s (cid:18) ξs (cid:19) (cid:34)(cid:18) T sH (cid:19) s (cid:35) , (A.22)Therefore, the final result is parametrically given by (cid:104) ζ k ζ q ζ p (cid:105) ( τ ) = c ξ (cid:16) P vac ζ (cid:17) k O (cid:18) T H (cid:19) δ (3) (cid:16) (cid:126)k + (cid:126)q + (cid:126)p (cid:17) , (A.23)where c is a small number that contains inverse powers of (2 π ) and we neglected a ln( T /H ) term. This32esult means that the non-Gaussian parameter f NL (cid:39) (cid:10) ζ (cid:11) / (cid:10) ζ (cid:11) would be f NL (cid:39) c ξ P vac ζ O (cid:18) T H (cid:19) . (A.24)One should now compare this with the result in the non-thermal regime, where f NL (cid:39) − e πξ P vac ζ /ξ [6], keeping in mind that the constant c could also be a small number. Let us look at a case of instantaneousthermalization at a temperature ¯ T , where the energy in the plasma is the same as the energy beforethermalization ρ γ (cid:39) − H e πξ /ξ . Therefore, ¯ T (cid:39) ρ / (cid:39) . He πξ/ /ξ / . This means that, even atthe very high temperature ¯ T , non-Gaussianity is proportional to e πξ and so it is suppressed by e − πξ compared to the non-thermal case. Using the observational bound f NL < O (10), we get ξ (cid:46) − − (cid:46) c (cid:46) − ). If we put instead the temperature T eq = ξ/ ¯ g , derived in section 6, due tothe competition between the thermal mass and the instability, we get ξ (cid:46) −
80, using the same rangeof c and using ¯ g = 0 . ζ propagator by itsthermal counterpart, which amounts to replace 1 / / N u . A simple bound can be derived bynoting that N u ( k ) ≤ T /H , for modes inside the horizon. Therefore, (cid:104) ζ k ζ q ζ p (cid:105) thermal1-loop (cid:46) d ξ (cid:16) P vac ζ (cid:17) k O (cid:18) T H (cid:19) δ (3) (cid:16) (cid:126)k + (cid:126)q + (cid:126)p (cid:17) × T H . (A.25)where d is another small constant. This implies that f thermal NL = (cid:10) ζ (cid:11) therm / ( (cid:10) ζ (cid:11) therm ) would then begiven by f thermal NL (cid:46) d ξ P vac ζ O (cid:18) T H (cid:19) . (A.26)This result is only an upper bound, while a proper calculation could turn out to give something smaller.Although this derivation is rough and a more proper study should be done, it shows parametricallythat the loop corrections are strongly suppressed compared to the vacuum case. Apart from the technicalcaveats hidden in the approximations there is another caveat to this computation. As we argued in thetext there are deviations for the thermal bath at low momenta and in particular we might expect a peakat those scales (in the case of an equilibrium temperature T eq , the peak would be located at ξH , as weargued in section 6). In this computation we neglected the effect of such a peak on the non-Gaussianity. A.7 Negligible superhorizon sourcing of curvature perturbation
In this section we show that the superhorizon sourcing of curvature perturbation due to the presence ofthe gauge fields is negligible. The amount of superhorizon sourcing is given by˙ ζ ( k ) = − H δp nad ρ + p (A.27)where δp nad is the non-adiabatic pressure perturbation, ρ is the total energy density and p the totalpressure. The non-adiabatic pressure in the case where the energy in the gauge fields is assumed to be33ubdominant is given by [32] δp nad = δp γ − ˙ p ˙ ρ δρ γ (A.28)where γ denotes the pressure and energy associated with the gauge fields. To leading order in slow-roll,and using p γ = 1 / ρ γ , then δp nad = 4 / ρ γ ( k ). On the other hand, ρ + p = 23 (cid:15)ρ φ + 43 ρ γ . (A.29)Therefore, assuming ρ γ (cid:28) (cid:15)ρ φ one gets ˙ ζ (cid:39) − H δρ γ (cid:15)ρ φ . (A.30)At horizon crossing we know that the gauge fields can imprint a sizable effect in ζ . Here, however,we are only interested in computing the sourcing on superhorizon scales. Therefore, we integrate theprevious equation from a time τ i where the modes are already superhorizon, − kτ i <
1. The superhorizoncorrection to the 2-point of ζ in Fourier space is then given by (cid:104) ζ k ( τ ) ζ p ( τ ) (cid:105) = (cid:104) ζ k ( τ ) ζ p ( τ ) (cid:105) ∗ + 4 (cid:90) ττ i dτ dτ τ τ (cid:15) ( τ ) ρ φ ( τ ) 1 (cid:15) ( τ ) ρ φ ( τ ) ( (cid:104) δρ γ ( k, τ ) δρ γ ( p, τ ) (cid:105) + c.c.) (A.31)where the star denotes the quantity evaluated at horizon crossing. For simplicity, and because that doesnot affect the leading order result, we assume (cid:15) and ρ φ to be roughly constant. Then, the computationboils down to compute (cid:104) δρ γ ( k ) δρ γ ( p ) (cid:105) where δρ γ ( k, τ ) = 12(2 π ) a (cid:90) d q (cid:2) A (cid:48) q A (cid:48) k − q + q | k − q | A q A k − q (cid:3) (cid:126)e ( (cid:126)q ) · (cid:126)e ( (cid:126)k − (cid:126)q ) . (A.32)and we consider here only the + polarization. When computing the 2-point function of δρ γ severalpermutations will appear but here we just want to show that the result is suppressed and so, for thatreason, we just use the magnetic part of δρ γ . Moreover, both in the thermal and non-thermal regime,the bulk of the energy is in modes with momenta − qτ (cid:38) H . Therefore, because we are evaluating thetime integrals at times where − kτ (cid:28) − qτ (cid:29) − kτ and so q − k (cid:39) q . Therefore, after performing the contractions we are left with (cid:104) δρ γ ( k, τ ) δρ γ ( p, τ ) (cid:105) ≈ τ τ H δ (3) ( (cid:126)p + (cid:126)k ) (cid:90) d q q [ A q ( τ ) A q ( τ ) ∗ ] . (A.33)We now specify the computation to the non-thermal case where the analytical solutions are easier tohandle. In the thermal case the reasoning is similar and would lead to a similar result.In the non-thermal case the integrals peak at horizon crossing which means when τ (cid:39) τ (cid:39) − /q .Therefore, by approximating the times to be the same we arrive at (cid:104) ζ k ( τ ) ζ p ( τ ) (cid:105) ≈ (cid:104) ζ k ( τ ) ζ p ( τ ) (cid:105) ∗ + 2 H (cid:15) ρ φ (cid:0) | A q | q (cid:1) ∗ (cid:90) ττ i dτ τ δ (3) ( (cid:126)p + (cid:126)k ) . (A.34)Therefore, while the first term is proportional to k − to give the standard scale invariance, the secondterms is proportional to τ where − kτ (cid:28)
1. Therefore we conclude that the superhorizon sourcing of34urvature perturbation is negligible.
References [1] A. Berera,
Warm inflation , Phys. Rev. Lett. (1995) 3218–3221, [ astro-ph/9509049 ].[2] A. Berera, I. G. Moss, and R. O. Ramos, Warm Inflation and its Microphysical Basis , Rept. Prog.Phys. (2009) 026901, [ arXiv:0808.1855 ].[3] M. Morikawa and M. Sasaki, Entropy Production in the Inflationary Universe , Prog. Theor. Phys. (1984) 782.[4] M.-a. Sakagami and A. Hosoya, Fate of Order Parameter in the Inflationary Universe , Phys. Lett. (1985) 342–346.[5] M. M. Anber and L. Sorbo,
Naturally inflating on steep potentials through electromagneticdissipation , Phys. Rev.
D81 (2010) 043534, [ arXiv:0908.4089 ].[6] N. Barnaby, R. Namba, and M. Peloso,
Phenomenology of a Pseudo-Scalar Inflaton: NaturallyLarge Nongaussianity , JCAP (2011) 009, [ arXiv:1102.4333 ].[7] A. Linde, S. Mooij, and E. Pajer,
Gauge field production in supergravity inflation: Localnon-Gaussianity and primordial black holes , Phys. Rev.
D87 (2013), no. 10 103506,[ arXiv:1212.1693 ].[8] R. Z. Ferreira and M. S. Sloth,
Universal Constraints on Axions from Inflation , JHEP (2014)139, [ arXiv:1409.5799 ].[9] R. Z. Ferreira, J. Ganc, J. Norea, and M. S. Sloth, On the validity of the perturbative description ofaxions during inflation , JCAP (2016), no. 04 039, [ arXiv:1512.06116 ]. [Erratum:JCAP1610,no.10,E01(2016)].[10] M. Peloso, L. Sorbo, and C. Unal,
Rolling axions during inflation: perturbativity and signatures , JCAP (2016), no. 09 001, [ arXiv:1606.00459 ].[11] A. Notari and K. Tywoniuk,
Dissipative Axial Inflation , JCAP (2016), no. 12 038,[ arXiv:1608.06223 ].[12] N. Bartolo, S. Matarrese, M. Peloso, and M. Shiraishi,
Parity-violating CMB correlators withnon-decaying statistical anisotropy , JCAP (2015), no. 07 039, [ arXiv:1505.02193 ].[13] C.-M. Lin and K.-W. Ng,
Primordial Black Holes from Passive Density Fluctuations , Phys. Lett.
B718 (2013) 1181–1185, [ arXiv:1206.1685 ].[14] E. Bugaev and P. Klimai,
Axion inflation with gauge field production and primordial black holes , Phys. Rev.
D90 (2014), no. 10 103501, [ arXiv:1312.7435 ].[15] L. Sorbo,
Parity violation in the Cosmic Microwave Background from a pseudoscalar inflaton , JCAP (2011) 003, [ arXiv:1101.1525 ]. 3516] K. Freese, J. A. Frieman, and A. V. Olinto,
Natural inflation with pseudo - Nambu-Goldstonebosons , Phys. Rev. Lett. (1990) 3233–3236.[17] A. Berera, Warm inflation solution to the eta problem , hep-ph/0401139 .[PoSAHEP2003,069(2003)].[18] S.-L. Cheng, W. Lee, and K.-W. Ng, Numerical study of pseudoscalar inflation with an axion-gaugefield coupling , Phys. Rev.
D93 (2016), no. 6 063510, [ arXiv:1508.00251 ].[19] P. Adshead, J. T. Giblin, T. R. Scully, and E. I. Sfakianakis,
Gauge-preheating and the end ofaxion inflation , JCAP (2015), no. 12 034, [ arXiv:1502.06506 ].[20] I. I. Tkachev,
Coherent scalar field oscillations forming compact astrophysical objects , Sov. Astron.Lett. (1986) 305–308. [Pisma Astron. Zh.12,726(1986)].[21] M. M. Anber and L. Sorbo, N-flationary magnetic fields , JCAP (2006) 018,[ astro-ph/0606534 ].[22] J. S. Schwinger,
On gauge invariance and vacuum polarization , Phys. Rev. (1951) 664–679.[23] T. Hayashinaka, T. Fujita, and J. Yokoyama, Fermionic Schwinger effect and induced current in deSitter space , JCAP (2016), no. 07 010, [ arXiv:1603.04165 ].[24] W. Tangarife, K. Tobioka, L. Ubaldi, and T. Volansky,
Dynamics of Relaxed Inflation , arXiv:1706.03072 .[25] P. Adshead, J. T. Giblin, T. R. Scully, and E. I. Sfakianakis, Magnetogenesis from axion inflation , JCAP (2016) 039, [ arXiv:1606.08474 ].[26] V. F. Mukhanov, H. A. Feldman, and R. H. Brandenberger,
Theory of cosmological perturbations.Part 1. Classical perturbations. Part 2. Quantum theory of perturbations. Part 3. Extensions , Phys.Rept. (1992) 203–333.[27] S. Bartrum, M. Bastero-Gil, A. Berera, R. Cerezo, R. O. Ramos, and J. G. Rosa,
The importanceof being warm (during inflation) , Phys. Lett.
B732 (2014) 116–121, [ arXiv:1307.5868 ].[28] A. D. Linde, V. Mukhanov, and M. Sasaki,
Post-inflationary behavior of adiabatic perturbationsand tensor-to-scalar ratio , JCAP (2005) 002, [ astro-ph/0509015 ].[29] R. Z. Ferreira and A. Notari,
In preparation , .[30]
Planck
Collaboration, P. A. R. Ade et al. , Planck 2015 results. XVII. Constraints on primordialnon-Gaussianity , Astron. Astrophys. (2016) A17, [ arXiv:1502.01592 ].[31] V. Mukhanov and S. Winitzki,
Introduction to quantum effects in gravity . Cambridge UniversityPress, 2007.[32] S. Nurmi and M. S. Sloth,
Constraints on Gauge Field Production during Inflation , JCAP (2014) 012, [ arXiv:1312.4946arXiv:1312.4946