Thermalized Axion Inflation

TL;DR

This paper studies axion-like inflation models with gauge field couplings, showing that gauge fields thermalize during inflation, which significantly alters the predicted cosmological perturbations and relaxes previous constraints on model parameters.

Contribution

It demonstrates that gauge fields can thermalize during inflation, leading to a thermal bath that modifies inflationary predictions and constraints, especially in models with Standard Model interactions.

Findings

Gauge fields thermalize at .9 for SM interactions.

Thermalization alters the spectrum of perturbations.

Tensor-to-scalar ratio is suppressed by H/(2T).

Abstract

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

Figures

Equations

$$S=\int d^{4}x\sqrt{-g}\left[\frac{M_p^2}{2}R+\frac{1}{2}{\partial }_{\mu }\varphi {\partial }^{\mu }\varphi -V(\varphi )-\frac{1}{4}{F}_{\mu \nu }{F}^{\mu \nu }-\frac{\varphi }{4f}{F}_{\mu \nu }{\tilde{F}}^{\mu \nu }\right].$$

$$A_{\pm }^{\prime \prime }+\left(k^2\pm \frac{2k\xi }{\tau }\right)A_{\pm }=0,\xi \equiv \frac{\dot{\varphi }}{2fH},$$

$$A_{+}(k,\tau )=\frac{1}{\sqrt{2k}}{e}^{\pi \xi /2}{W}_{-i\xi ,1/2}(2ik\tau ),$$

$$N_{\text{BE}}(k)={\left(e^{-\frac{\omega (k)}{aT}}-1\right)}^{-1},$$

$$1/2+N_X(k)\equiv \frac{k^2|X_k{|}^2+|X_k^{\prime }{|}^2}{2\omega (k)}.$$

$$u_k^{\prime \prime }+\left(k^2-\frac{2}{\tau }\right)u_k=0,$$

$$N_{{\gamma }_{+}}^{\prime }(k)$$

$$N_{u_{}}^{\prime }(k)$$

$$N_{{\gamma }_{-}}^{\prime }(k)$$

$$u=\frac{e^{ik\tau }}{\sqrt{2k}}\left(1+\frac{i}{k\tau }\right),g_u=\frac{i\left(k^2{\tau }^2+ik\tau -1\right)}{\tau (k\tau +i)},$$

$$N_{{\gamma }_{+}}^{\prime }(k)$$

$$N_{u_{}}^{\prime }(k)$$

$$N_{{\gamma }_{-}}^{\prime }(k)$$

$$S_k=\frac{1}{\omega (k)}\int \prod _{i=2}^4\left(\frac{d^3{\vec{k}}_i}{(2\pi {)}^3(2{\omega }_i)}\right){\left|M_i\right|}^2(2\pi {)}^4{\delta }^{(4)}\left(k^{\mu }+k_2^{\mu }-k_3^{\mu }-k_4^{\mu }\right)B(k,k_2,k_3,k_4),$$

$$B(k_1,k_2,k_3,k_4)$$

$${⟨u_k{u}_k⟩}_{\text{loop}}(\tau )$$

$$\frac{d{N}_u(k)}{d\tau }\simeq -\frac{8\xi b}{(2\pi {)}^3(2f{)}^2{a}^{2}k}\times$$

$$\times \int \frac{dq{q}^3|\vec{k}-\vec{q}|}{\text{min}(q,|\vec{k}-\vec{q}|)}{N}_u(k)(1+N_{{\gamma }_{+}}(q))(1+N_{{\gamma }_{+}}(|\vec{k}-\vec{q}|))-N_{{\gamma }_{+}}(q)N_{{\gamma }_{+}}(|\vec{k}-\vec{q}|)(1+N_u(k)),$$

$$S^{++}=\frac{1}{\omega (k_1)}\int d{k}_{2}d{k}_3\frac{c_S}{f^4}\left[N_{{\gamma }_{+,3}}{N}_{{\gamma }_{+,4}}(1+N_{{\gamma }_{+,1}})(1+N_{{\gamma }_{+,2}})-N_{{\gamma }_{+,1}}{N}_{{\gamma }_{+,2}}(1+N_{{\gamma }_{+,3}})(1+N_{{\gamma }_{+,4}})\right],$$

$$N_{{\gamma }_{+,1}}{N}_{{\gamma }_{+,2}}+N_{{\gamma }_{+,1}}{N}_{{\gamma }_{+,2}}{N}_{{\gamma }_{+,3}}+N_{{\gamma }_{+,1}}{N}_{{\gamma }_{+,4}}{N}_{{\gamma }_{+,2}}-$$

$$-N_{{\gamma }_{+,3}},N_{{\gamma }_{+,4}}{N}_{{\gamma }_{+,2}}-N_{{\gamma }_{+,1}}{N}_{{\gamma }_{+,3}}{N}_{{\gamma }_{+,4}}-N_{{\gamma }_{+,3}}{N}_{{\gamma }_{+,4}}.$$

$$S^{++}\approx \frac{{\omega }^5}{{\beta }_S{f}^4}{N}_{{\gamma }_{+}}^3.$$

$$N_{{\gamma }_{+}}\gg \sqrt{{\beta }_S\frac{H{f}^4}{{\omega }^5}}.$$

$$\xi \gtrsim 0.45\ln \left(\frac{f}{H}\right)+2.7,$$

$$D^{+\varphi }$$

$$D\approx \frac{{\omega }^3}{{\beta }_D{f}^2}{N}_{{\gamma }_{+}}^2,$$

$$N_{{\gamma }_{+}}(\omega )\gg {\beta }_D\frac{H{f}^2}{{\omega }^3}⟹\xi \gtrsim 0.45\ln \left(\frac{f}{H}\right)+4.3,$$

$$N_{{\gamma }_{+}}\gg \sqrt{{\beta }_{-}\frac{H{f}^4}{{\omega }^5}}⟹\xi \gtrsim 0.45\ln \left(\frac{f}{H}\right)+2.2,$$

$$N_{{\gamma }_{+}}\gg \sqrt{{\beta }_{\varphi }\frac{H{f}^4}{{\omega }^5}}⟹\xi \gtrsim 0.45\ln \left(\frac{f}{H}\right)+2.4,$$

$${\frac{d\sigma }{d\Omega }}_{\gamma \gamma \to e^{-}{e}^{+}}=\frac{2\pi {\alpha }^2}{(4\omega {)}^2}\left(\frac{1+{\cos }^2\theta }{{\sin }^2\theta }\right),$$

$${\sigma }_{\gamma \gamma \to \gamma \gamma }^{\varphi }=\frac{|M{|}^2}{64{\pi }^2(2\omega {)}^2}\simeq \frac{1}{4}{\left(\frac{\omega }{f}\right)}^4\frac{1}{64{\pi }^2(2\omega {)}^2}.$$

$$\frac{{\sigma }_{\gamma \gamma \to e^{-}{e}^{+}}}{{\sigma }_{\gamma \gamma \to \gamma \gamma }^{\varphi }}\simeq \frac{\frac{2\pi {\alpha }^2}{(4\omega {)}^2}}{\frac{1}{4}{\left(\frac{\omega }{f}\right)}^4\frac{1}{64{\pi }^2(2\omega {)}^2}}=128{\pi }^3{\alpha }^2\frac{f^4}{{\omega }^4},$$

$$N_{{\gamma }_{+}}\gg \frac{8}{\pi {\alpha }^2},$$

Full text

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

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

1 Introduction

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

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

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

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

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

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

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

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

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

2 The axial coupling

Let us consider an axion-like scalar field $\phi$ with a potential $V(\phi)$, coupled to a gauge field $A_{\mu}$, with field strength $F_{\mu\nu}$, via an axial coupling, as described by the action

[TABLE]

Here $\tilde{F}^{\mu\nu}=\epsilon^{\mu\nu\alpha\beta}F_{\alpha\beta}/(2\sqrt{-g})$ is the dual of $F_{\mu\nu}$. Note that the above coupling induces a periodic potential, of period $f$, only in the non-abelian case, while in the abelian case $V(\phi)$ can be completely unrelated to $f$. Actually, in the context of this paper the only relevant ingredient is the axial coupling to gauge fields, and so our considerations apply also to the case of a non-periodic potential.

We split the scalar field into a spatially homogeneous background value, which drives inflation, and a perturbation, $\phi=\phi(\tau)+\delta\phi(\vec{x},\tau)$ and we approximate the metric with a de-Sitter form in conformal coordinates $ds^{2}=a^{2}(\tau)[d\tau^{2}-d\vec{x}^{2}]$, where $a=(-H\tau)^{-1}$ and $\tau$ goes from $-\infty$ (past) to $0^{-}$ (future). Here $H\equiv a^{\prime}/a^{2}$ is the Hubble constant during inflation and a prime denotes the derivative taken with respect to conformal time $\tau$. In such a background, the gauge field of comoving momentum $k$ satisfies the following equation of motion in Coulomb gauge (see [20] and references therein):

[TABLE]

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

We consider the simple case $\xi\simeq{\rm constant}$, which is a good approximation if $\phi$ is slowly-rolling down its potential. It is also a good approximation if we are in the regime of strong backreaction of the gauge fields on $\phi$ [5, 11]. In both cases $\xi=\sqrt{\epsilon/2}M_{p}/f$, where $\epsilon\equiv\dot{\phi}^{2}/(2M^{2}_{p}H^{2})$ is the first slow-roll parameter. The equation of motion then has analytical solutions that can be written for example by a Whittaker function:

[TABLE]

which, in rough terms, connects the flat space oscillatory regime deep inside the horizon666Note however that even at $\tau\rightarrow-\infty$ the mode functions always have a logarithmic time-dependent phase. with an exponential growth for $-k\tau<2\xi$, until the solution approaches a constant $A_{k}=e^{\pi\xi}/(2\sqrt{\pi k\xi})$ well outside the horizon, when $-k\tau\lesssim(8\xi)^{-1}$.

3 Thermalized Axion Inflation

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

Crucially the scattering probabilities are Bose-enhanced if the number of $\gamma_{+}$ is large, which will happen when $\xi\gtrsim 1$. As a consequence, if the particle number overcomes some threshold, which depends on $\xi$ and $f$, the scatterings will lead to an equilibrium distribution and as a result the modes will be redistributed from $k\lesssim 2\xi H$ to a new scale given by the temperature $T$. Roughly speaking, when the rates $\Gamma_{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, $\gamma_{+}$ and $\gamma_{-}$, as well as $\phi$ or other charged particles (depending on the size of the cross sections) and would lead, a priori, to a Bose-Einstein (BE) distribution defined as

[TABLE]

where $\omega(k)=\sqrt{k^{2}+a^{2}m^{2}}$ is the comoving energy of the particle. Note, however, that $\gamma_{+}$, due to the imaginary dispersion relation modes and the associated sourcing of low momenta modes $k/a<2\xi H$, would be infinitely populated if one could wait an arbitrarily long time. In our inflationary background each mode spends a finite amount of time in the instability region and so it never reaches an infinite occupation number but, as a remnant of such dynamics, we expect relevant deviations from a BE distribution, such as the presence of a peak at low momentum. On the other hand, $\gamma_{-}$ should still be described by a BE distribution with a modified dispersion relation with $\omega^{2}=k^{2}+m^{2}$, where $m^{2}=-2k\xi/\tau$ while the $\phi$ 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 $\rho_{R}\approx 10^{-4}H^{4}e^{2\pi\xi}/\xi^{3}$ [21], due to the continuous excitation of modes, as described by eq. (2.2). A simple estimate of the expected temperature at thermalization is roughly given by ${\bar{T}}\approx\rho_{R}^{1/4}\approx 0.1He^{\pi\xi/2}$. Note that in our system three scales of interest are present: the horizon scale $H$, the instability scale $2\xi H$ (each mode starts getting excited when its momentum satisfies $k/a<2\xi H$) and the temperature at thermalization ${\bar{T}}$. If $\xi\gtrsim{\cal O}(1)$ there is a hierarchy of scales: $H<2\xi H<\bar{T}$. We will consider $\xi$ in the range 1$\sim$10: in fact, if it gets very large backreaction will start and one should consider the dynamics of the background together with the mode evolution [5, 11], which we postpone to future work. In any case, even in presence of backreaction $\xi$ should be related logarithmically to the parameters of the potential and typically be at most ${\cal O}(10)$ [5].

After thermalization is reached, however, the system can evolve in a non-trivial way. In fact, having an interacting plasma typically implies that gauge fields have thermal masses, which should screen the instability and, as a result, quench the particle production, so that the temperature is expected to decrease. Moreover, the energy extraction from the scalar field is proportional to $\dot{\phi}/(4f)\langle F\tilde{F}\rangle=a^{-4}\int kd^{3}d/d\tau(|A_{+}|^{2}-|A_{-}|^{2})$ [5, 11]; such a quantity should be suppressed if $\gamma_{+}$ and $\gamma_{-}$ are in equilibrium, since they tend to compensate each other when averaged over a thermal distribution. One could imagine that the system could reach a stationary configuration at a temperature smaller than $\bar{T}$, or perhaps an oscillatory behavior. We postpone the study of the full evolution after thermalization to future work, although we will comment on some expected features in section 6.

3.1 Boltzmann equation and particle numbers

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

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

For a given field $X$, with mode functions $X_{k}$, we define a comoving particle number as the ratio between energy density $\rho_{k}$ and energy per particle $\omega(k)$:

[TABLE]

Deep inside the horizon, in the Bunch-Davies vacuum, the fields have a clear particle interpretation, since they behave as $X_{k}=e^{ik\tau}/\sqrt{2k}$ at $\tau\rightarrow-\infty$, with vacuum particle number $N_{X}=0$. In the scalar field case the previous definition will be valid for the canonically normalized field $u\equiv a\delta\phi$, which satisfies

[TABLE]

where we neglected, for simplicity, slow-roll corrections. In appendix A.4 we also show that the definition of $N_{X}(k)$ matches the more standard one given in terms of the Bogolyubov coefficient $N_{k}=|\beta_{k}|^{2}$.

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

Using the above definition of effective particle number we are able to rewrite the equations of motion, eq. (3.3) and eq. (2.2), as an equation for the number of right-handed gauge fields, $N_{\gamma_{+}}(k)$, and $u$ particles, $N_{u}(k)$:

[TABLE]

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

[TABLE]

while for $g_{A}$ we can use the exact solution eq. (2.3).

3.2 Collision terms

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

[TABLE]

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

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

• •

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

• •

the particle number entering in the collision terms is given by eq. (3.2) with $\omega(k)\simeq k$, as discussed above;

• •

backreaction effects of gauge fields on $\phi$ are neglected for simplicity.

Scatterings

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

We consider the scatterings shown in fig. 1. Each scattering term $S_{k}$ has the form

[TABLE]

where $k_{i}^{\mu}=(\omega_{i},\vec{k}_{i})$ are the momenta of the external legs, $k_{i}=|\vec{k_{i}}|$, $k^{\mu}=(\omega,\vec{k})$, $M_{n}$ is the matrix element of the process (given in appendix A.3) and $B$ are the phase space factors

[TABLE]

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

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

Decays

In flat space the decay rate of a massive particle of mass $m$ with an axial coupling is simply given by $\Gamma_{d}\propto m_{\phi}^{3}/f^{2}$, 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 $\gamma_{+}$ or, in other words, due to the energy transfer from the background to the fields.

We estimate such a decay rate by looking at the 1-loop correction to the 2-point function $\langle u_{k}u_{k}\rangle$, with gauge fields $\gamma_{+}$ running in the loop. Using the fact that the two point function is related to the mode functions via $\langle u_{k}u_{k^{\prime}}\rangle=|u_{k}|^{2}(2\pi)^{3}\delta^{3}(\vec{k}-\vec{k^{\prime}})$, we can relate the loop correction to an increase in the particle number and thus to an inverse decay rate999Note that using the Bogolyubov coefficients to define the particle number, as explained in sec. A.4, one also finds that the two point function is related to the particle number as $\langle u_{k}u_{k}\rangle\propto 1+2N_{k}+2\text{Re}\left[\beta_{k}\alpha_{k}^{*}\right]$ where $\alpha_{k}$ is the other Bogolyubov coefficient.

To compute the loop correction we use the in-in formalism by considering the interaction Hamiltonian $H_{\text{int}}=\int d^{3}x\sqrt{-g}\,\delta\phi\,F\tilde{F}/(4f)$. After Fourier transforming and integrating over the delta functions according to appendix A.1 we get

[TABLE]

where we omitted a factor of $\delta^{3}(0)$, $\vec{e}_{+,-}$ are the polarization vectors defined in A.4 and we used $q\equiv|\vec{q}|$. Here […] stands for a commutator and $\left<...\right>$ is an expectation value. Then, by taking the time derivative of this expression and by making several approximations, described in appendix A.2, we arrive at

[TABLE]

where $b$ is a number of order ${\cal O}(10^{-3})$ related to the angular integrals. Written in this form we can now identify the right hand side with the decay term $D^{+\phi}$ of eq. (3.9).

3.3 Analytical estimates

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

• •

only subhorizon modes participate in the collisions, which is possible if $\xi\gtrsim 1$;

• •

the collision integrals peak where the $N_{\gamma_{+}}$ is maximal, which means at energies $H\lesssim\omega\lesssim 2\xi H$;

• •

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

For the gauge fields the dominant scatterings are $\gamma_{+}\gamma_{+}\leftrightarrow\gamma_{+}\gamma_{+}$ because they are enhanced by the most powers of $N_{\gamma_{+}}$. After integrating the delta functions over the angles the scattering term has the form (we omit here factors of $a$):

[TABLE]

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

[TABLE]

We assume the integrals to peak at energies $\omega_{i}\simeq\omega$, somewhere between $H$ and $2\xi H$, where the modes are more densely populated, with $N_{\gamma_{+,i}}\equiv N_{\gamma_{+}}\gg 1$. In this limit the above expression goes as $N_{\gamma_{+}}^{3}$, hence, we estimate

[TABLE]

Here $\beta_{S}={\cal O}(10^{4})$ is a numerical factor that comes from the collision integrals and the cross sections. The scatterings become relevant when this term is comparable with the left-hand side of eq. (3.8) which is of order $N_{\gamma_{+}}H$ (or than the source term, which is of the same order). Therefore, the condition in the particle number to have thermalization is

[TABLE]

In the next section we compare this estimate with the numerical results. This can also be translated into a bound on $\xi$, for a given $f/H$, by using the fact that, in absence of thermalization, the particle number depends exponentially on $\xi$. We can estimate the total particle number in the band $1<-k\tau<2\xi$ by numerically integrating $N_{\gamma_{+}}\approx 2k|A_{+}|^{2}$, using the exact mode functions and fitting with an exponential. This yields $N_{\gamma_{+}}\approx 10^{-4}e^{4.5\xi}$. Using this expression in the above condition for thermalization, and setting $\omega\approx H$, then gives

[TABLE]

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

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

[TABLE]

where $c_{D}$ is a dimension 3 coefficient. Similarly to the scatterings, the decay term can be approximated by

[TABLE]

where we assumed $N_{\gamma_{+}}\gg N_{\phi}$. Here $\beta_{D}\approx 10^{5}/\xi$ is the inverse of the pre-factor of eq. 3.14. Comparing $D$ with $N_{\gamma_{+}}H$ gives

[TABLE]

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

Finally we can estimate the threshold of thermalization for $\phi$ and $\gamma_{-}$. In the case of $\gamma_{-}$ the collision terms can be estimated as $S^{+-}\approx\omega^{5}N_{\gamma_{+}}^{2}N_{\gamma_{-}}/(\beta_{-}f^{4})$, with $\beta_{-}={\cal O}(10^{2})$, which should be compared to the left-hand side of the Boltzmann equation, of order $N_{\gamma_{-}}H$. This leads to the condition

[TABLE]

which is actually realized even more easily than eq. (3.18), since $\beta_{-}\ll\beta_{S}$. The case of $\phi$ is also analogous and leads to

[TABLE]

where $\beta_{\phi}\approx 5\cdot 10^{2}$ and so it is also realized almost at the same time as the above eq. (3.23).

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

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

3.4 Standard Model couplings

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

The purpose of this section is to study whether the SM self-interactions and scatterings involving gauge bosons and SM charged fields could help in thermalizing the system. The interesting feature of such interactions is that they are not suppressed by powers of $1/f$ like the previous diagrams. Instead, they are proportional to gauge coupling constants and so, if $f\gg H$, they are the relevant interactions for thermalization. Moreover, this makes the scenario more predictive since there is only one parameter ($\xi$) that sets the thermalization condition for the gauge bosons. However, as we mentioned before, in this case the $\phi$ 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 leave a deeper and more precise study with all such interactions for future work. Nevertheless, we can still estimate the onset of thermalization due to scatterings with SM particles and compare it with the interactions generated from the axial coupling.

Photon-Fermion scatterings

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

[TABLE]

where $\alpha=e^{2}/(4\pi)$, $e$ is the electric charge and $2\omega$ is the center-of-mass energy. The total cross section has a $\log$ divergence for small angles which can be regulated by imposing an IR cutoff, in our case given by the Hubble scale $H$, or, in the case of thermalization, by a thermal mass of order $eT$. Now we can compare this expression with photon-photon scattering mediated by $\phi$, which in the center-of-mass frame is

[TABLE]

Comparing the two cross sections we have

[TABLE]

where we assumed the logarithmic term to give a coefficient of order one and the typical energy of the process to be $\omega$. It is clear that, if $f\gg H$, the cross section for pair production is much larger than the scatterings mediated by $\phi$. In particular if we extrapolate the Standard Model to high energies the relevant coupling is the U(1) hypercharge, whose $\alpha$ changes from $1/40$ at $10^{14}$ GeV to $1/60$ at $10^{2}$ GeV. Therefore, the ratio of cross-sections is always above one for any $f>\omega$, which is anyway required for the validity of effective field theory.

On the other hand, pair production is less Bose-enhanced than photon-photon scattering by one power of $N_{\gamma_{+}}$. To estimate the threshold for thermalization we compare $H$ with $n_{\gamma_{+}}\sigma_{\gamma\gamma\rightarrow e^{-}e^{+}}$ where $n_{\gamma_{+}}\approx N_{\gamma_{+}}H^{3}$ is roughly the particle number density. Therefore, the thermalization condition is

[TABLE]

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

We can also compare these scattering rates with the rate of particle production by the Schwinger effect given by $\Gamma_{\text{S}}\simeq\frac{e^{2}E^{2}}{6\pi^{3}H^{3}}\exp\left({\frac{-\pi m_{e}^{2}}{eE}}\right)$ [22, 23, 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^{2}\simeq N_{\gamma_{+}}H^{4}$ and consider the regime $E\gg m^{2}_{e}$, we find a similar threshold for thermalization in terms of $\xi$.

Gauge-field self interactions

Another case in which thermalization of gauge fields can happen quite easily is if the gauge bosons belong to $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 $\propto g_{s}^{2}$, with $g_{s}$ the gauge group coupling constant. The cross-section of gluon-gluon scattering is121212We neglect factors of order 1 in the total cross-section.

[TABLE]

where $\alpha_{s}=g_{s}^{2}/(4\pi)$. Thus, the condition for thermalization is

[TABLE]

For $SU(3)$, whose coupling constant $\alpha_{s}=g_{s}^{2}/(4\pi)$ runs from $1/40$ at $10^{14}$ GeV to $\simeq 1$ at GeV scale 131313Note that this analysis is only valid for gauge-fields in the perturbative regime where the propagator and the equation of motion remains unmodified by self-interactions., this means that gluon scatterings are the most efficient process for thermalization and the system would thermalize also around $\xi\simeq 2.9$.

4 Numerics

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

• •

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

[TABLE]

where $\Delta k$ is the mode spacing. We typically use 10 discrete modes, multiples of a given $k_{\text{min}}$.

• •

We only consider subhorizon modes, where the flat space result and the particle interpretation are a good approximation. We start our simulation when the largest mode $k_{\text{min}}$ is subhorizon, using $|k_{\text{min}}\tau|=2$, and stop it when such mode goes superhorizon ($|k_{\text{min}}\tau|<1$).

• •

We give as initial condition for $\gamma_{+}$ the solution of the equation of motion eq. (2.2). For $\gamma_{-}$ and $\phi$ one can use the vacuum solution and eq (3.7), respectively, although the final results will basically be insensitive to such choice.

• •

To estimate the temperature of the system at thermalization, $\bar{T}$, one would need to ensure that the chosen window of modes covers the bulk of the energy, which consists of modes of momenta $k/a\lesssim\bar{T}$. However, $\bar{T}$ grows exponentially with $\xi$ and the instability band consists of modes of $k/a$ between $H$ and $\xi H$. Thus, for large $\xi$ we would need a very large window of momenta to extract correctly the temperature. For numerical reasons we do not consider such a large box but only a band of modes between $H$ and ${\cal O}(10)\,H$ approximately. For this reason we cannot evaluate correctly $\bar{T}$ when $\xi$ is large, although the system still approaches a Bose-Einstein distribution in the chosen window.

• •

We neglect any backreaction on the scalar field equation of motion. This would be relevant when $\langle F\tilde{F}\rangle/f$ becomes of order of the derivative of the potential $dV/d\phi$ and would require to follow the treatment of [5, 11] and to include in the system the background equation for $\phi$. We postpone this richer 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), to be always given by the solution of the equations of motion in absence of collisions. This assumption might fail after thermalization, but should be valid until there.

• •

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

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

• •

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

[TABLE]

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

• •

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

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

Let us start by considering only $\gamma_{+}$ gauge fields and their scatterings, with matrix elements given in eq. (A.11). In fig. 2 we show the total particle number and the average difference to a thermal distribution $\Delta N/N$, defined in eq. (4.1), for fixed $f$ and varying $\xi$ and vice-versa, evaluated at the end of the simulation (when the longest mode is horizon size). We assume the system has thermalized when $\Delta N/N\lesssim 0.1$. When $\xi=3.9$, $\Delta N/N$ decreases sharply to roughly $0.1$ for $f\lesssim 20H$, while for $f=H$ that happens for $\xi\gtrsim 2.7$. This is in agreement with the plots on the right which show that the analytical estimate for the threshold of thermalization, using eq. (3.18), intersects the total particle number $N^{\text{total}}_{\gamma_{+}}=\sum_{k}N_{\gamma_{+}}(k)$ at roughly the same values of $f$ and $\xi$. In the upper right plot we also see that the particle number depends exponentially on $\xi$ and is well approximated by $N_{\gamma_{+}}\simeq 3\times 10^{-4}e^{4.8\xi}$, as expected from the analytical estimates in sec. 3.1.

In the second case we consider all scatterings and decays involving $\phi,\gamma_{+}$ and $\gamma_{-}$. The presence of decays and of more degrees of freedom leads to a decrease in the $\gamma_{+}$ particle number which implies less efficient thermalization, i.e, it happens for smaller (larger) value of $f$ ($\xi$). In fig. 3 we show the same plots as in fig. 2. By looking at the average differences to a Bose-Einstein distribution (left plots) one can see that for $\xi=3.9$, lower plots, the two gauge field polarizations thermalize at $f\lesssim 4H$ while for $f=H$ all species thermalize when $\xi\gtrsim 3.9$. Comparing the average differences, on the left, with the particle numbers and estimates, on the right, we see that thermalization requires a value of $\xi$ about 0.6 larger and 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 that the power spectrum of both polarizations of the gauge field and of the inflaton perturbations follow each other while, at the same time, there is a transfer of energy to the UV (fig. 7 of [25]).

Finally, we also evaluate numerically the threshold for thermalization by collecting the values of $\xi$ and $f$ such that $\Delta N/N<0.1$. The threshold is well fitted by

[TABLE]

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

5 Phenomenology

After thermalization is reached, however, the evolution is expected to be non-trivial. Indeed the temperature should decrease and the system could also depart from a thermal spectrum at low momenta. In fact, having a temperature normally generates thermal masses for the gauge fields, which tend to screen the instability and, as a result, the source will be less efficient and the temperature is expected to decrease. Moreover, the energy density $\rho_{\phi}$ extracted from the scalar field is given by $d\rho_{\phi}/dt=\dot{\phi}/(4f)\langle F\tilde{F}\rangle=\dot{\phi}/f\int(d^{3}k)/(2\pi)^{3}kd/d\tau(|A_{+}|^{2}-|A_{-}|^{2})/a^{4}$, averaged over a thermal distribution; such a quantity should be suppressed if $\gamma_{+}$ and $\gamma_{-}$ are in equilibrium, since the scatterings tend to restore parity symmetry. We do not try to estimate such quantity exactly here because, while $N_{\gamma_{-}}$ is easy to compute in equilibrium, the same is not true for $N_{\gamma_{+}}$, due to its complex frequency $\omega_{+}$. We postpone to future work a refined numerical treatment, including such effects, but we anticipate in this section some of the expected features, under the assumption that thermalization is successful. Generically we can anticipate how the scalar and tensor curvature perturbation, $\zeta$ 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 density decreases as $a^{-4}$. Because the sourcing of adiabatic curvature perturbation is proportional to this quantity, we expect a negligible isocurvature sourcing of the scalar curvature perturbation in the uniform density gauge, $\zeta$. In appendix A.7 we elaborate on this point. Therefore, it should be a good estimate to evaluate correlators of $\zeta$ at horizon crossing. We also disregard here the possibility that a backreaction regime, which could change the predictions on perturbations, can be present at any time during inflation.

• •

Loop corrections:

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

• •

Parity symmetry:

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

5.1 Power spectrum and tensor to scalar ratio

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

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

[TABLE]

which means that the power spectrum of the scalar curvature perturbation in the uniform density gauge, $\zeta$, in the superhorizon limit, $-k\tau\rightarrow 0$, is [26]

[TABLE]

Alternatively, one could have evaluated the subhorizon expression for $u$ at horizon crossing, $-k\tau\simeq 1$, and plugged it into the definition of $\zeta$, which is constant on superhorizon:

[TABLE]

where $*$ denotes quantities evaluated at horizon crossing. This last result is actually more general because $H_{*}$ and $\dot{\phi}_{*}$ are evaluated at horizon crossing and so it also holds in the quasi de Sitter case.

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

If the inflaton perturbations are instead thermal, the expected outcome is much simpler by following the same strategy of evaluating $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^{\prime}_{k}|^{2}\simeq k|u_{k}|^{2}$, we have

[TABLE]

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

[TABLE]

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

Indeed the spectral index is given by:

[TABLE]

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

[TABLE]

The reason why we introduced $\epsilon_{H}$ is to keep the analysis fully general: in the backreaction dominated case 141414Note, however, that in the backreaction case one might also worry that $\zeta$ could have a non-negligible contribution from gauge field fluctuations at horizon crossing, which is not taken into account by the above eq. (5.5)., in fact, this parameter differs from the previous definition, $\epsilon=\dot{\phi^{2}}/(2H^{2}M_{P}^{2})$.

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

[TABLE]

This case is phenomenologically very interesting because it would help polynomial large field models ($V(\phi)\propto\phi^{n}$, with $n>1$), but also others models of inflation where $\epsilon$ and $\eta$ are comparable, to agree with observations. In fact the observational bound $r\lesssim 0.1$ usually puts a quite stringent bound on $\epsilon$, while in our case if $T_{*}/H_{*}$ is sufficiently large the bound is relaxed and those models are not anymore in tension with observations.

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

[TABLE]

In this case the tensor-to-scalar ratio would remain unchanged, $r=16\epsilon$, while the tensor tilt would be non-trivial:

[TABLE]

This regime would actually be very interesting. In fact, normally the tensor tilt is necessarily red ($n_{T}<0$), reflecting the fact that the Hubble constant decreases during inflation (and so $\epsilon_{H}>0$); here, instead, we could find a blue spectrum depending on the behavior of $T_{*}/H_{*}$. It is however difficult to have thermalization of tensors, due to their weaker coupling, as we will see in section 6.

5.2 Constraints on $f$ and $\xi$

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

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

[TABLE]

where $l$ is some loop factor. For $l\simeq 10^{2}$ and imposing $P_{\zeta}^{\text{vac}}=2.2\times 10^{-9}$, the previous bound constrains $\xi\lesssim 3.5$ although specific 1-loop computations suggest a weaker bound $\xi\lesssim 4.4$ [9, 10].

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

[TABLE]

In the case without collisions we have $\langle F\tilde{F}\rangle\simeq 10^{-4}H^{4}e^{2\pi\xi}/\xi^{4}$ [5], which means that backreaction is absent if $f/H\gg 4\times 10^{-3}e^{\pi\xi}/\xi^{5/2}$.

However, we stress here that these constraints do not directly apply if the system is thermalized simply because the occupation numbers completely change and new constraints should be derived. In fact, 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 whether the spectrum of perturbations is thermal or not, the vacuum spectrum cannot exceed the observational bound, $P_{\zeta}^{\text{vac}}\leq P_{\zeta}^{\text{obs}}\equiv 2.2\times 10^{-9}$. The reason is that vacuum fluctuations are a lower bound and adiabatic perturbations cannot be erased once they cross the horizon [28]. If we now recall the definition of $\xi$, this imposes a lower bound on $f$,

[TABLE]

In fig. 4 we plot the different thresholds for thermalization derived in sections 3.3 and 4. Jointly with the thermalized regions we overlap the perturbative constraint eq. (5.11) (derived in absence of thermalization) for $l=10^{2}$, the backreaction threshold (also in absence of thermalization) eq. (5.12) and the observational constraint eq. (5.13). In absence of SM interactions the phenomenological allowed region for $\phi$ to be thermalized is such that one would need to take into account higher order loop corrections to the propagators as well as backreaction, at least before the onset of thermalization. Therefore, a dedicated study is required. However, we stress that in the presence of SM interactions the system can thermalize before backreaction although, in this case, it is also crucial to take into account thermal masses to establish when $\phi$ thermalizes and when backreaction is relevant. We address this issue in the next section.

Finally, we argue why, generically, the constraints on $\xi$ are very much alleviated when the system thermalizes. The reason is two-folded. First, the fact that energy moves from the horizon scale to the UV, with a consequent decrease in the particle number, changes the peak of the integrals from the horizon scale to the temperature scale. This will affect all correlators. Moreover, odd correlators are further suppressed by the following reason. The interaction Hamiltonian associated with the axial coupling is given by [6, 8]

[TABLE]

Therefore, all odd correlators will necessarily be proportional to, at least, one power of $\left<F\tilde{F}\right>\propto d/d\tau(N_{+}-N_{-})$. In a thermal bath $N_{+}$ and $N_{-}$ tend to equilibrate through scatterings and so any parity asymmetry tends to be suppressed. This asymmetry would be completely erased if the two helicities had exactly the same dispersion relation. However this is not the case, since the axial coupling changes the dispersion relations of $\omega_{+}$ and $\omega_{-}$ in different ways, and so a small asymmetry should remain there, proportional to $\xi H$, even in thermal equilibrium.

Computing loop corrections in the thermal regime requires a dedicated study [29]. Nonetheless, in appendix A.6 we provide a rough derivation of the parametric dependence of the non-Gaussian parameter $f_{NL}\simeq\left<\zeta^{3}\right>/\left<\zeta^{2}\right>^{2}$ with the temperature, considering the case where gauge fields are thermalized but $\phi$ is not. We find

[TABLE]

where $c$ is a small number containing inverse powers of $(2\pi)$ and the result of the angular integration and so needs to be derived in a more accurate computation. If we consider the case of instantaneous thermalization of the plasma, then, the energy density in the gauge fields is simply given by the initial energy $\rho_{\gamma}\simeq 10^{-4}H^{4}e^{2\pi\xi}/\xi^{3}$ [21]. Therefore, the plasma would have a temperature $\bar{T}\simeq\rho^{1/4}\simeq 0.1He^{\pi\xi/2}/\xi^{3/4}$ which means that $f_{NL}$ would be proportional to $e^{2\pi\xi}$. Comparing this result with the non-thermal case where $f_{NL}\simeq 10^{-7}e^{6\pi\xi}P_{\zeta}^{\text{vac}}/\xi^{8}$ [6] there is a parametrical suppression of $e^{-4\pi\xi}$ which translates into a weaker constraint for $\xi$. In particular, assuming $10^{-4}>c>10^{-7}$ and using the constraint $f_{NL}<{\cal O}(10)$ [30] one would find $\xi\lesssim 6-7$. However such a high temperature $\bar{T}$ is actually not what we finally expect, as we discuss in the next section. In reality, we expect a much smaller temperature to be reached, thus relaxing even more the constraints. In appendix A.6 we also argue that, if $\phi$ is also thermalized, $f_{NL}$ can be at most multiplied by an additional power of $T/H$

[TABLE]

where $d$ is another small coefficient.

6 Effect of thermal masses

So far we have studied the conditions for thermalization to occur due to the initially large occupation number of gauge bosons. However, it is more subtle to understand what happens after thermalization is established.

As we discussed in previous sections, due to the presence of an imaginary dispersion relation for $\gamma_{+}$, it is unclear how is the shape of the distribution at low momenta, $k/a<2\xi H$. In fact such modes are exponentially produced, although in a finite window of time, and so one expects relevant deviations from a BE distribution such as the presence of a peak at low momentum. The $\gamma_{-}$ polarization, instead, should be correctly described by a BE distribution with a dispersion relation $\omega^{2}=k^{2}+m^{2}$, where $m^{2}=-2k\xi/\tau$, while $\delta\phi$, if thermalized, is described by a massless BE distribution.

However a very important ingredient modifies the above discussion: when thermalization happens, gauge fields typically161616This certainly happens in the well-known case of having interactions with charged SM particles. It could also happen for the case of only $\phi$-mediated interactions, but this would require performing the one-loop thermal correction, which we postpone to a subsequent publication [29]. develop a thermal mass $m_{T}\propto gT$, where $g$ is the gauge coupling, as a result of having a non zero density of particles correcting the propagator at 1-loop. Note, importantly, that $\phi$ instead does not get a thermal mass, since its interactions only renormalize its kinetic term [29]. We leave to the future a more careful study of the effect of thermal masses, but we anticipate here some interesting features.

Let us parameterize the thermal mass as

[TABLE]

In the case of SM fermions $\bar{g}^{2}=g^{2}\sum_{i}q^{2}_{i}/6$, where $q_{i}$ is the charge of a given fermion. At energies of around $10^{14}$ GeV this gives $\bar{g}^{2}\simeq 0.3$. If the thermal mass is included in the equation of motion for the gauge fields it will compete with the instability due to the axial coupling, giving

[TABLE]

When $m_{T}\geq\xi H$ the mass completely shields the instability band. Therefore, gauge fields cannot be produced and we might expect an equilibrium temperature

[TABLE]

In fact, starting from a non-thermal case, when we reach thermalization the plasma has an initial temperature, $\bar{T}$, typically higher than $T_{eq}$. Thus, at that point the particle production stops and the plasma of particles is expected simply to redshift as radiation, as $a^{-4}$, which is equivalent to say that temperature decreases as $1/a$. When the temperature drops below $T_{eq}$ an instability band reappears, most likely as a narrow band centered at $|k\tau|\approx\xi H$, which is the minimum of $\omega_{T}(k)^{2}$. At this point we expect the system either to reach an almost stationary state slightly below $T_{eq}$ or perhaps an oscillatory behavior around such temperature.

When the temperature drops, thermalization could be less efficient. We can estimate the condition for the plasma to remain in thermodynamic equilibrium by comparing $\sigma_{eq}\cdot n_{eq}$ with $H$, using the equilibrium number density $n_{eq}\simeq T^{3}$ and the typical cross section $\sigma_{eq}\simeq\bar{\alpha}^{2}/T^{2}$, where $\bar{\alpha}^{2}\equiv g^{4}/(4\pi)^{2}\sum_{i}q^{4}_{i}$ is summed over different species. In the SM case at $10^{14}$ GeV, $\bar{\alpha}^{2}\simeq 0.05$. As a result we get:

[TABLE]

To estimate if the plasma keeps thermalized down to $T_{eq}$ we should fulfill the condition

[TABLE]

For a given $g$ and a given set of fermions, this condition can indeed be met for sufficiently large $\xi$ while when considering only SM fermions we get a more precise condition on $\xi$. Note however that, even if the condition is not met in the range $\bar{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\tau|\approx\xi H$ to be constantly present, generating a peak in the $\gamma_{+}$ distribution, whereas we expect a thermal distribution away from such narrow peak, where $\omega_{T}(k)^{2}>0$. Observing this dynamics should be numerically feasible but we postpone further investigation of this important issue to future work. We now briefly discuss possible interesting observational consequences of such cases.

We have to distinguish between the case with SM interactions and the case with only $\phi$-mediated interactions. In the first case, as we saw before, we can reach a thermal state if $\xi\gtrsim 2.9$. Then we need to check if $\phi$ thermalization is efficient by using, as above, the thermal condition $\sigma_{eq}\cdot n_{eq}\gg H$, with $\sigma_{eq}\approx T^{2}/f^{4}$, so that:

[TABLE]

which should be true at $T=T_{eq}$, thus imposing

[TABLE]

If the above condition is met we can use the results from the previous section to make predictions for the spectrum of perturbations.

If $\phi$ is thermal from eq. (5.5) we have that:

[TABLE]

and the spectral index becomes simply

[TABLE]

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 $\zeta$ is conserved on superhorizon scales.

In the case of tensor modes the couplings are gravitational, suppressed by $1/M_{p}$, and so the cross section would be given by $\sigma_{eq}\approx T^{2}/M_{p}^{4}$ and so in order to reach thermalization one would need $T\gg(HM_{p}^{4})^{1/5}$. This requires a temperature larger than the inflationary energy, $\rho_{inf}^{1/4}=3H^{2}M_{p}^{2}$, and so, for that reason, it is not possible to thermalize the tensor modes during inflation at this equilibrium temperature. Therefore, assuming tensor modes to be in the vacuum the tensor to scalar ratio would be given by

[TABLE]

As we mention before, thermalization of $\phi$ leads to a suppression of the tensor to scalar ratio which, at the equilibrium temperature, would be $\bar{g}/(2\xi)$. If we look at fig. 5, this amounts to at least an ${\cal O}(10^{-2})$ suppression.

We can also try to estimate whether there is relevant backreaction on the scalar field equation of motion, as follows. The energy extraction from $\phi$ is given by $\langle F\tilde{F}\rangle\dot{\phi}/f$ which we can estimate by using conservation of energy density in the gauge fields $\rho_{\gamma}$:

[TABLE]

In a stationary situation $\dot{\rho}_{\gamma}\approx 0$ while backreaction is achieved if $\langle F\tilde{F}\rangle/(4f)\gtrsim 3H\dot{\phi}$. These two conditions then require

[TABLE]

Finally, using the estimation of non-Gaussianity in the case where $\phi$ is thermal, eq. (A.26), we can also find what would be a conservative viable region for $T_{eq}=\xi/\bar{g}$ by requiring the non-Gaussian parameter $f_{NL}<{\cal O}(10)$. This would impose $\xi\lesssim 20-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 including the effects of a thermal mass and assuming a stationary regime. The different regimes are described by eqs. (6.5), (6.7), (6.12), (6.8) and (5.13). Note that, contrary to fig. 4, the region where $\phi$ is thermalized is not completely inside the backreaction region and so there is some parameter space where $\phi$ can thermalize without backreaction.

In the case of oscillatory temperatures around $T_{eq}$ the above description remains probably roughly correct although in that case we should expect larger deviations from the thermal predictions as well as superimposed oscillations in the relevant quantities at horizon crossing. Clearly, if such oscillations are large the model is observationally ruled out. However, interestingly, if oscillations are within the observational bound but non-negligible, this would result in oscillatory power spectra, which could be another striking evidence of this mechanism.

7 Conclusions

In this work we have analyzed the system composed of an inflaton ($\phi$) coupled to a gauge field ($A_{\mu}$) through a CP-odd term with strength $1/f$. One helicity of the gauge field, $\gamma_{+}$, is known to develop an instability with exponential particle production at a rate controlled by the parameter $\xi=\dot{\phi}/(2fH)$, for modes with momenta in the range $(8\xi)^{-1}H\lesssim k/a\lesssim 2\xi H$. We have included in this scenario the effect of scatterings through a set of Boltzmann-like equations. The scattering rates are Bose-enhanced by the large occupation numbers of $\gamma_{+}$ and can drive the system to thermal equilibrium, thus leading to a setup that we dubbed thermalized axion inflation. We have analyzed two situations: (1) only scalar-gauge field interactions, as in the case of a “dark photon”; (2) identifying the gauge field with a Standard Model gauge boson (abelian or non-abelian). We stress that (2) is probably the best motivated case, since this system anyway needs to be coupled to the SM to provide reheating, and it turns out to be the one in which we have good control of the theory, since thermalization happens before backreaction and before perturbativity constraints.

In the first case the most relevant scatterings are $\gamma_{+}\gamma_{+}\rightarrow\gamma_{+}\gamma_{+}$. By solving the associated Boltzmann-like set of equations of the system we find that $N_{\gamma_{+}}$, $N_{\gamma_{-}}$ and $N_{\phi}$ evolve to a Bose-Einstein distribution when $\xi\gtrsim 0.44\ln(f/H)+3.4$. Note, however, that due to the axial coupling the gauge fields have a modified dispersion relation, imaginary in the case of $\gamma_{+}$, which is expected to distort the distribution for low momenta $k/a<2\xi H$. In a phenomenologically realistic model one has to require the power spectrum of curvature perturbations to be compatible with observations, and that imposes a lower bound on $f/H$, which then implies thermalization at $\xi\gtrsim 6$. At such values, however, we would have to assume that the occupation numbers are not significantly altered by loop effects [9], which are important already at $\xi\gtrsim 3.5$. Moreover, for such a large $\xi$ the gauge fields would also backreact on $\phi$ and on $\delta\phi$, so a dedicated study 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 easier to achieve because the cross sections are not suppressed by powers of $H/f$. This case is also more predictive because $\xi$ is the only free parameter. For example, if the gauge field is the one associated to $U(1)$ hypercharge and considering only particle anti-particle productions, we find that thermalization happens already at $\xi\gtrsim 2.9$. A similar condition on $\xi$ applies if instead we consider the gauge field to be a gluon and include self-interactions. This is of great interest because in this case one can have a thermal bath of particles well under control, before perturbativity constraints become relevant and also before backreaction of gauge fields on $\phi$ becomes important.

In both cases thermalization has profound implications for the phenomenology of axion inflation. To elaborate on this we have further pointed out that, after thermalization is reached, the system should evolve in a very different regime, due to the presence of thermal masses $m_{T}=\bar{g}T$ for the gauge bosons, which tends to screen the instability of $\gamma_{+}$. As a result we expect either an oscillatory behavior or a stationary solution at a temperature given by $T_{eq}=\xi H/\bar{g}$. Note that this regime is now only linearly dependent on $\xi$ and not exponential anymore and so all constraints on $\xi$ should become much weaker and should be properly readdressed. We described such behavior only qualitatively in this work but we plan to address these points in full detail in the future.

Assuming the system reaches a stationary solution, or one where the temperature changes adiabatically, we derived very interesting consequences for cosmological observables such as the scalar and tensor perturbations, under the assumption that $\phi$ is thermalized, which can happen in the region shown in fig. 5. In this regime one can also find a region of parameters where $\phi$ can thermalize without backreacting on the background. Note also that in this case the power spectrum of curvature perturbations can be obtained simply by setting its value at horizon crossing to be that of a thermal state, assuming it will be conserved on superhorizon scales. This leads to a result for $P_{\zeta}$ that is enhanced compared to the vacuum case by the ratio of particle numbers between the thermal state, $T_{*}/H_{*}$, and the vacuum, $1/2$. We have also computed the tilt of the power spectrum to be $n_{s}-1=4\epsilon_{H}-\eta$, which differs from the standard slow-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 possible for them to thermalize, at least at the equilibrium temperature. For that reason we expect the tensors to remain in the vacuum and so, in the case in which $\phi$ is thermalized, a tensor to scalar ratio suppressed by $H_{*}/(2T_{*})$ with respect to the standard formula. This is very interesting because it can allow for a reconciliation between models of inflation where $\epsilon$ and $\eta$ are comparable, e.g. polynomial large field models, with observations.

Finally, although we did not present here a detailed study of non-Gaussianities, we provided strong arguments, and one particular estimation, which justify why the previous constraints do not apply to our case and why, generically, we expect them to be less restrictive in the thermal regime. Thermalization redistributes in fact the occupation numbers: it depletes the highly populated horizon-sized modes and it populates higher $k$ modes. As we have argued, loop corrections to cosmological correlators will then be generically smaller than in the non-thermal case. In addition, odd correlators of $\zeta$ are further suppressed, since thermodynamic equilibrium tends to restore parity, driving $\langle F\tilde{F}\rangle$ to small values. We estimated the non-Gaussian parameter $f_{NL}$ to be proportional to171717This estimation holds if $\phi$ is not thermalized, while in the opposite case we showed that there can be, at most, one extra power of $T/H$. $P_{\zeta}^{\text{vac}}\,T^{4}/H^{4}$, which corresponds to a large parametric suppression of $e^{-4\pi\xi}$ compared to the non-thermal case. This seems promising since it could allow for interesting regimes with large $\xi$, such as the $\phi$-thermalized regime and the backreacting regime [29], to be in agreement with observations.

To sum up we provided a working model in which during inflation a thermal bath with a possibly large temperature 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 shift symmetry 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 can be achieved in full control it could even remove the standard need of an inflaton flat potential, by the help of the dissipative friction. Moreover reheating is already incorporated in this scenario and obtained just 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. Tobioka and M. Sloth for comment and useful discussions. This work is supported by the grants EC FPA2010-20807-C02-02, AGAUR 2009-SGR-168, ERC Starting Grant HoloLHC-306605 and by the Spanish MINECO under MDM-2014-0369 of ICCUB (Unidad de Excelencia “Maria de Maeztu”).

Appendix A Appendix

A.1 Fourier conventions

We quantize scalars ($\delta\phi$ or $\zeta$) and gauge fields as

[TABLE]

where we have used the Coulomb gauge ($A_{0}=0$) and where the creation and annihilation operators satisfy

[TABLE]

The polarization vectors satisfy the following identities

[TABLE]

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 an estimation of the decay rate. In this appendix we present some of the intermediate steps which took us to eq. (3.14). Starting from eq. (3.2) by taking the time derivative of the two point function we arrive at

[TABLE]

The time integrations are dominated by the latest time simply because, before thermalization, that is when 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 replacing the time integration by its integrand evaluated at $\tau$ times the integration interval which we approximate to be $\Delta\tau=2\xi/\text{min}(q,|\vec{k}-\vec{q}|)$ corresponding to the time at which the largest mode enters the resonant band and so where the process becomes non-zero. Thus, we arrive at

[TABLE]

Our goal is to rewrite this expression as powers of 2-point functions which we then associate with particle numbers. However, by isotropy the particle number should not depend on the angle, so it is convenient to perform the angular integral. At this level we do not have access to the generic mode functions in the presence of scatterings, so we parameterize this unknown by $b$ which was found in similar 1-loop computations to be ${\cal O}(10^{-3})$. 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 $|\vec{k}-\vec{q}|$. We assume $|\vec{k}-\vec{q}|\simeq k+q$. Under these simplifications we arrive at

[TABLE]

where we omitted the time argument in the right hand side because all quantities are evaluated at $\tau$. By a similar reasoning the first term in the right hand side of the previous equation gives

[TABLE]

We only deal with subhorizon modes, thus, the first term is smaller or at most equal to the second term so we neglect it. Now, in order to write the expression as a particle number we approximate $|A^{\prime}_{k}|\simeq kA_{k}$ and make use of the approximate relation $1/2+N_{u}(k)\approx k|u_{k}|^{2}$. The commutator structure implies the expectation value to give an expression of the form

[TABLE]

The imaginary part of the two point function can be roughly seen as the vacuum contribution, $(2\pi)^{3}$, while the real part would be proportional to $(2\pi)^{3}(2N)$. This identification is accurate for thermal propagators and in fact by making this identity we would find that the previous equation has indeed a structure of particle numbers similar to those appearing in a standard decay expression which would be proportional to

[TABLE]

However, the 1-loop correction contains more processes than just the decay and for that reason we do not expect to arrive at that precise form. Nevertheless we can isolate the terms which do have the form of a decay and so we finally end up with

[TABLE]

The extra factor of $4$ comes from identifying $\text{Re}[\left<XX\right>]\propto 2N_{X}$.

A.3 Matrix elements for scatterings

In this appendix we provide the matrix elements for the scatterings depicted in fig. 1. We start by the $\gamma_{+}\gamma_{+}\rightarrow\gamma_{+}\gamma_{+}$ scattering which only has non-zero contributions from the s-channel

[TABLE]

where $\epsilon^{abcd}$ is the 4-dimensional Levi-Civita tensor and $k^{\mu}$ are the 4-momenta of the particles involved in the scatterings. The processes involving $\gamma_{+}$ and $\gamma_{-}$ can be readily obtained from the previous result by making use of the identities in eq. (A.4) and summing over all the channels $u,s,t$. On the other hand the gauge field-axion scattering, $\gamma_{+}\phi\rightarrow\gamma_{+}\phi$, has instead s and t-channel contributions and is given by

[TABLE]

where $\eta_{dz}$ is the Minkowski metric. Finally gauge field-axion conversion, $\gamma_{+}\gamma_{+}\rightarrow\phi\phi$, has matrix elements

[TABLE]

A.4 Particle number and Bogolyubov coefficients

In curved space the particle number is typically defined as $N_{k}=|\beta_{k}|^{2}$, where [31]

[TABLE]

is the Bogolyubov coefficient which relates the mode functions between two different basis of states and $W_{k}=(f^{\prime}(k)f^{*}(k)-f^{\prime*}(k)f(k))/(2i)$ is the Wronskian. If we want to know the particle number compared to the vacuum state we fix $f_{1}=e^{-ik\tau}/\sqrt{2k}$ and so

[TABLE]

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 of momenta and at the end of the simulation, $-k_{\text{min}}\tau=1$, for $\xi=3.9$ and $f=H$. In the same plot we show Bose-Einstein distributions fro massless particles with and without a chemical potential. We can see that, apart from the longest mode, a Bose-Einstein distribution with $T=45H$ is in rough agreement with the distribution. The inclusion of a chemical potential $\mu=10^{-1.66}T$ gives a perfect fit, including also 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 from the 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 the case where the photons are thermalized but the inflaton perturbations are not. We also derive an upper bound for the case with thermalized inflaton perturbations. This estimation is also easily generalizable to higher order correlators. Let us start by the 3-point function. The 1-loop correction is proportional to

[TABLE]

where $H_{\text{int}}$ is the interaction Hamiltonian given by [6, 8]

[TABLE]

We are interested in the parametric dependence of the correlators with the temperature in order to compare it to previous results. For that reason we will ignore factors of $2\pi$ and other ${\cal O}(1)$ numbers. We will assume the result to be dominated by the thermal part of the gauge field propagators. We also neglect the commutator structure, which should be an overestimation, since there could be cancelations due to its structure.

Let us first analyze the vertex structure. Generically, the 1-loop correction to the n-point function will be proportional to

[TABLE]

As we mentioned in the text, in the thermal regime we expect $\left<F\tilde{F}\right>$ to be suppressed because parity tends to be restored. In that case, the correlators involving the gauge fields will be dominated by $\left<FF\right>$ terms. By noting that $\tilde{F}_{\mu\nu}\tilde{F}^{\mu\nu}=F_{\mu\nu}F^{\mu\nu}$ we estimate

[TABLE]

up to permutation factors which we neglected here. The term in $\left<FF\right>$ is dimensionally an energy density and, using eq. (3.2), we approximate in the thermal regime $\left<FF\right>\simeq 4kN(\omega_{+}(k))$, barring cancellations. On the other hand, the term in $\left<F\tilde{F}\right>$ is proportional to the parity asymmetry $\left<F_{\mu\nu}(k)\tilde{F}^{\mu\nu}(k)\right>=kd/d\tau(|A_{+}|^{2}-|A_{-}|^{2})\simeq d/d\tau(N(\omega_{+}(k))-N(\omega_{-}(k))$ where $\omega_{+,-}(k)=(k^{2}\pm 2k\xi/\tau)^{1/2}$. Assuming a thermal distribution we find ,in limit of large temperature, and for $-k\tau\gg\xi$

[TABLE]

Note, however, that in the in-in computation the vertices are evaluated at different times and so they could connect different physical momenta. Nevertheless, because the integrals are regulated in the UV by the particle number and in the IR by the momentum dependence of the integrand, we expect the integrals to peak when the gauge fields running in the loop have physical momenta of the order of the temperature $s_{\text{phys}}\simeq T$.

Under these approximations and after Fourier transforming all the interaction Hamiltonians and integrating over the delta functions, the 1-loop correction simplifies to

[TABLE]

where we assumed the integrals to peak in the equilateral configuration, $k\simeq q\simeq p$, similarly to previous results in the non-thermal case [6], although there might be some changes in the thermal regime. In the last expression we have also replaced $\zeta$’s in the time integrals by their subhorizon expression $\zeta_{k}=H^{2}\tau/(\dot{\phi}\sqrt{k})$, while the external ones by $\zeta_{k}=H^{2}/(\dot{\phi}\sqrt{k^{3}})$, because they are evaluated at horizon crossing. Finally, to keep the treatment simple we further assume the integral to peak at a scale $s\gg k,q,p$. We could also include corrections with $s\simeq k,q,p$ but we expect not to have higher powers of the temperature so we neglect them for this estimation.

The time integrals should peak when the $s/a=-sH\tau=T$. At that scale the particle number is $N_{+,-}(T)\simeq 1$ and so we get

[TABLE]

Therefore, the final result is parametrically given by

[TABLE]

where $c$ is a small number that contains inverse powers of $(2\pi)$ and we neglected a $\ln(T/H)$ term. This result means that the non-Gaussian parameter $f_{NL}\simeq\left<\zeta^{3}\right>/\left<\zeta^{2}\right>^{2}$ would be

[TABLE]

One should now compare this with the result in the non-thermal regime, where $f_{NL}\simeq 10^{-7}e^{6\pi\xi}P_{\zeta}^{\text{vac}}/\xi^{8}$ [6], keeping in mind that the constant $c$ could also be a small number. Let us look at a case of instantaneous thermalization at a temperature $\bar{T}$, where the energy in the plasma is the same as the energy before thermalization $\rho_{\gamma}\simeq 10^{-4}H^{4}e^{2\pi\xi}/\xi^{3}$. Therefore, $\bar{T}\simeq\rho^{1/4}\simeq 0.1He^{\pi\xi/2}/\xi^{3/4}$. This means that, even at the very high temperature $\bar{T}$, non-Gaussianity is proportional to $e^{2\pi\xi}$ and so it is suppressed by $e^{-4\pi\xi}$ compared to the non-thermal case. Using the observational bound $f_{NL}<{\cal O}(10)$, we get $\xi\lesssim 6-7$ (assuming $10^{-7}\lesssim c\lesssim 10^{-4}$). If we put instead the temperature $T_{eq}=\xi/\bar{g}$, derived in section 6, due to the competition between the thermal mass and the instability, we get $\xi\lesssim 40-80$, using the same range of $c$ and using $\bar{g}=0.5$.

In the case where inflaton perturbations are thermalized one should replace the $\zeta$ propagator by its thermal counterpart, which amounts to replace $1/2$ by $1/2+N_{u}$. A simple bound can be derived by noting that $N_{u}(k)\leq T/H$, for modes inside the horizon. Therefore,

[TABLE]

where $d$ is another small constant. This implies that $f^{\text{thermal}}_{NL}=\left<\zeta^{3}\right>^{\text{therm}}/(\left<\zeta^{2}\right>^{\text{therm}})^{2}$ would then be given by

[TABLE]

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 parametrically that the loop corrections are strongly suppressed compared to the vacuum case. Apart from the technical caveats hidden in the approximations there is another caveat to this computation. As we argued in the text there are deviations for the thermal bath at low momenta and in particular we might expect a peak at those scales (in the case of an equilibrium temperature $T_{eq}$, the peak would be located at $\xi H$, as we argued 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 of the gauge fields is negligible. The amount of superhorizon sourcing is given by

[TABLE]

where $\delta p_{\rm nad}$ is the non-adiabatic pressure perturbation, $\rho$ is the total energy density and $p$ the total pressure. The non-adiabatic pressure in the case where the energy in the gauge fields is assumed to be subdominant is given by [32]

[TABLE]

where $\gamma$ denotes the pressure and energy associated with the gauge fields. To leading order in slow-roll, and using $p_{\gamma}=1/3\rho_{\gamma}$, then $\delta p_{\rm nad}=4/3\rho_{\gamma}(k)$. On the other hand,

[TABLE]

Therefore, assuming $\rho_{\gamma}\ll\epsilon\rho_{\phi}$ one gets

[TABLE]

At horizon crossing we know that the gauge fields can imprint a sizable effect in $\zeta$. Here, however, we are only interested in computing the sourcing on superhorizon scales. Therefore, we integrate the previous equation from a time $\tau_{i}$ where the modes are already superhorizon, $-k\tau_{i}<1$. The superhorizon correction to the 2-point of $\zeta$ in Fourier space is then given by

[TABLE]

where the star denotes the quantity evaluated at horizon crossing. For simplicity, and because that does not affect the leading order result, we assume $\epsilon$ and $\rho_{\phi}$ to be roughly constant. Then, the computation boils down to compute $\left<\delta\rho_{\gamma}(k)\delta\rho_{\gamma}(p)\right>$ where

[TABLE]

and we consider here only the $+$ polarization. When computing the 2-point function of $\delta\rho_{\gamma}$ several permutations will appear but here we just want to show that the result is suppressed and so, for that reason, we just use the magnetic part of $\delta\rho_{\gamma}$. Moreover, both in the thermal and non-thermal regime, the bulk of the energy is in modes with momenta $-q\tau\gtrsim H$. Therefore, because we are evaluating the time integrals at times where $-k\tau\ll 1$ then, the integrals will peak at momenta $-q\tau\gg-k\tau$ and so $q-k\simeq q$. Therefore, after performing the contractions we are left with

[TABLE]

We now specify the computation to the non-thermal case where the analytical solutions are easier to handle. 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 $\tau_{1}\simeq\tau_{2}\simeq-1/q$. Therefore, by approximating the times to be the same we arrive at

[TABLE]

Therefore, while the first term is proportional to $k^{-3}$ to give the standard scale invariance, the second terms is proportional to $\tau^{3}$ where $-k\tau\ll 1$. Therefore we conclude that the superhorizon sourcing of curvature perturbation is negligible.