Gravitational wave background from non-Abelian reheating after axion-like inflation

A pseudoscalar inflaton $\varphi$, coupled to the topological charge density $F\tilde{F}$ of a non-Abelian sector, can decay to gauge bosons ($\varphi\to g g $), which may thermalize rapidly. The friction felt by $\varphi$ is then increased by non-Abelian"strong sphalerons", leading to a self-amplifying process that can efficiently heat up the medium. We determine a lower bound for the gravitational wave production rate from such a process, originating via hydrodynamic fluctuations and particle collisions, in terms of a minimal number of parameters. Only a moderate fraction of energy density is converted to gravitational waves, suggesting that non-Abelian models may avoid the overproduction observed in some Abelian cases.


Introduction
The dynamics of the early universe is dictated, on one hand, by gravitational physics, accounting for a period of inflationary expansion that produced the seeds for structure formation. On the other hand, the microscopic properties of matter are governed by particle physics. The two descriptions connect to each other during reheating, in which the vacuum energy density, sustaining inflation, is converted to thermal radiation, carried by particles. The process establishes the highest temperature, T max , that can meaningfully be talked about.
It would be interesting to understand theoretically the dynamics of reheating, and to develop empirical tests for it. As we are studying early moments, penetrating probes need to be considered, with gravitational waves as an obvious candidate [1][2][3][4]. As tensor perturbations, gravitational waves are also produced by inflation itself [5]. However, the spectra that originate from inflation and reheating have different shapes, with the latter peaking at higher frequencies, all the way to the microwave range MHz...GHz. This could be within the reach of future observation, even if novel avenues need to be explored (cf. ref. [6] for a review).
The purpose of this paper is to study reheating within a framework that is simple enough to be tractable without heavy numerics, yet rich enough to offer for variants with different observational consequences.
Specifically, we consider a pseudoscalar inflaton field ϕ [7][8][9], whose interactions are governed by the Lagrangian where F c µν ≡ ∂ µ A c ν − ∂ ν A c µ + gf cde A d µ A e ν is the Yang-Mills field strength, N c is the number of colours, c ∈ {1, ..., N 2 c − 1}, g 2 ≡ 4πα is the gauge coupling, and f a is the axion decay constant. In the context of inflation, the allure of this model stems from its incorporation of interactions with Standard Model-like gauge fields, but in the special way that they do not spoil the desired flatness of the potential. That said, eq. (1.1) does involve a nonrenormalizable operator, so to keep the description self-consistent, it can only be applied for low energy and momentum scales of ϕ, i.e. ǫ ϕ ≪ max(4πf a , f a /α). The inflationary predictions of eq. (1.1) depend on the shape of the potential V (ϕ). A much-used example mimics an instanton-induced periodic structure, where on the semiclassical level f b = f a (this may change after renormalization). 1 However, there are other possibilities, for example as given in eqs. (4)-(8) of ref. [10]. As we are concerned with a late stage, we may expand the potential around its global minimum (ϕ min = 0). This then leads to a universal shape that is also familiar from chaotic inflation [11], According to eqs. (1.1) and (1.3), our framework depends on the parameters f a , m, α. It would be easy to enrich the setup by introducing more parameters, for instance via the shape of V (ϕ). Furthermore, if the plasma contains fermions (ψ k , k = 1, ..., N f ), then it is natural to couple ϕ also to the corresponding pseudoscalar operators, ∂ µ [ψ k γ µ γ 5 ψ k ] and im kψk γ 5 ψ k , where m k is a fermion mass. The coefficient of each operator is, a priori, independent. If the fermions couple to the gauge fields A a µ , then there is a certain redundancy in the couplings, as dictated by the axial anomaly equation. At the same time, fermions affect the dynamics of the gauge sector, with anomalous processes generating an effective chemical potential for left and right-chiral modes, which influences the friction felt by ϕ [12]. However, in the following, we restrict ourselves to the minimal setup of eqs. (1.1) and (1.3).
Our presentation is organized as follows. We start with a brief review of approaches to reheating in sec. 2, motivating the idea that in the non-Abelian case the density matrix of the system can be parametrized with a rapidly increasing temperature-like variable. The main part is sec. 3: after reviewing how gravitational waves are produced from such a system (sec. 3.1), we compute the contributions at long (sec. 3.2) and short wavelengths (sec. 3.3), and illustrate the results numerically (sec. 3.4). The findings are summarized in sec. 4.

Mechanism of sphaleron-induced reheating
Recently, a large body of literature has appeared in which the gauge fields in the operator χ of eq. (1.1) have Abelian nature, and contribute then to the gravitational wave background (an incomplete list can be found in refs. ). Physically, Abelian fields could represent either a "dark photon", or the Standard Model hypercharge gauge field. In the following, we restrict ourselves instead to non-Abelian gauge fields, which is arguably more natural if one thinks of embedding the Lagrangian in some Grand Unified framework. 2 The methods employed in the above studies follow, roughly speaking, two lines. One relies on the solution of mode equations (cf., e.g., ref. [8]), the other on classical field theory simulations (cf., e.g., refs. [45][46][47]). While simulations should account for the full non-perturbative dynamics of momentum modes with large occupation numbers, they are not sensitive to phenomena where the occupation number is of order unity. In particular, the quantum mechanical vacuum decays that dominate at early stages, or the thermalization towards the Bose distribution that takes place at the end, are not captured.
The framework that we adopt is a different one, modelling the complicated dynamics of the non-Abelian sector through the assumption that the system effectively thermalizes. In practice, this is closely related to the paradigm of warm inflation [48], well established in the axion inflation context [49][50][51][52][53][54]. Our new implementation guarantees that the initial vacuum-like decays are correctly incorporated as well [55].
The physical thinking behind this framework is that non-Abelian gauge fields tend to thermalize more rapidly than other types of interactions, as both particle number and momenta are changed at each cubic vertex. In fact, the heat bath has been argued to represent an attractor solution even during the inflationary stage [56,57]. General aspects of non-Abelian thermalization have been reviewed in ref. [58].
For the case of eq. (1.1), the friction felt by ϕ takes a special form. The reason is that, in the limit of low frequencies, the real-time topological susceptibility is rendered non-zero through the non-perturbative dynamics mediated by "strong sphalerons" [59]. Incorporating this physics, an interesting warm inflation scenario has been found [60][61][62].
In concrete terms, adopting the signature (+−−−), and splitting overall energy-momentum conservation into field and radiation parts [63], the equations of motion take the form where T is a local temperature; u µ is a local plasma four-velocity; Υ is a friction coefficient; and e r and p r are respectively the energy density and pressure of radiation.
The key issue is to fix Υ, as it originates from the operator in eq. (1.1). This problem was addressed in ref. [55], where it was shown that the value is determined by the response function of the medium, evaluated at an appropriate frequency scale, ω. Thereby both vacuum physics (for ω ≫ πT ) and sphaleron physics (for ω < ∼ α 2 T ) are incorporated.
For the present investigation, we can for the most part treat Υ as an independent parameter. Its dynamical form in terms of f a , m, α and T can be inserted at the end (cf. sec. 3.4). We also keep in mind that, as demonstrated in ref. [55] in the context of the "weak regime" of warm inflation, the system can reheat up to temperatures T max ∼ f a /α after inflation, where the influence of the operator χ in eq. (1.1) becomes of order unity.

Gravitational wave production from reheating
As already alluded to, an important physics phenomenon associated with the reheating process is the production of a gravitational wave background [1][2][3][4]. In the Abelian case, this can be used for constraining the model parameters (cf., e.g., ref. [10] and references therein). In the current section, we study gravitational wave production within the framework of sec. 2.
The empirical significance of a gravitational wave background depends on its wavelength. The physics of inflation corresponds to large wavelengths, or low frequencies, ≪ Hz today. In contrast, gravitational waves produced at reheating peak at high frequencies, which are in fact similar to those originating from a thermal plasma [64][65][66].
To carry out the computation, we organize it in an "adiabatic" approximation, assuming that plasma reactions are fast compared with the Hubble rate (α 2 T > H) and considering physical momenta well within the horizon (k ≫ H). The first part is certainly satisfied during and after reheating [55]. Under these conditions we can operate in a local Minkowskian frame. 3 It is generally believed that non-equilibrium processes produce more gravitational waves than equilibrium ones, as the latter reflect the maximal entropy of the underlying system, so that we should not be able to discern many features. Therefore thermal production as considered here is likely to set a lower bound for the full rate. 3 Recently, the regime of typical inflationary momenta exiting the horizon (k < H) has been considered [67].
The analysis of hydrodynamic fluctuations bears a conceptual similarity to sec. 3.2, however the shear viscosity was taken over from self-interacting scalar field theory [68] rather than from interactions with a heat bath. In principle our formalism could be generalized to apply to such a situation and then be used for addressing the full range of momenta relevant for warm inflation (cf., e.g., ref. [69] for a concrete realization), however the approximation of a local Minkowskian frame needs to be abandoned for the momentum modes that exit the horizon, even if the basic physical processes remain the same.

General features
In a local Minkowskian frame, the production rate of the energy density carried by gravitational radiation can be expressed as [64,65] where f GW is the polarization-average phase-space density of gravitons, and the dot stands for a time derivative. On general grounds [70], the evolution equation for f GW takes the forṁ where n B (k) ≡ 1/(e k/T −1) is the Bose distribution, and m Pl ≈ 1.221×10 19 GeV is the Planck mass.
In practice, f GW ≪ n B (k), so the right-hand side can be approximated as Γ(k) n B (k).
The dynamical information about the processes taking place is encoded in the interaction rate Γ(k), which in turn can be expressed as [65] is the retarded correlator related to the energy-momentum tensor T µν ,

4)
... T denotes a thermal average, and Ä is a transverse and traceless projector, In practice, it is sometimes convenient to choose the special frame in which k = k e z , whereby rotational symmetry implies that 4 The method to compute Im G R xy;xy depends on the momentum range considered. For very small momenta, we find ourselves in the so-called hydrodynamic domain. Then Im G R xy;xy 4 An algebraic way to verify the prefactor is given in footnote 2 of ref. [65]. Physically, there are two transverse-traceless polarizations, and if we choose the one in non-diagonal components as a representative, viz. (T xy + T yx )/ √ 2, then the equality T xy = T yx leads to the additional factor ( evaluates to ηω, where η is the shear viscosity [64]. The task therefore is to determine the contribution of the inflaton field ϕ to the shear viscosity, a topic that we address in sec. 3.2.
In contrast, for typical thermal momenta, parametrically k ∼ πT , elementary gauge bosons and axions can be resolved as quasi-particles. Then we are faced with a Boltzmann type of a computation, which is presented in sec. 3.3.

Hydrodynamic domain
We start by computing the correlator in eq. (3.7) in the hydrodynamic domain, i.e. at very small frequencies and momenta. 5 In this regime, elementary particles cannot be resolved, and the degrees of freedom relevant for describing the gauge plasma are hydrodynamic fluctuations. The contribution of a weakly coupled scalar field to eq. (3.7) has been determined in ref. [71] in this domain, and here we improve upon that computation, by ameliorating its ultraviolet (UV) sensitivity.
In the hydrodynamic domain, the traceless part of the energy-momentum tensor reads where T r µν is the contribution of radiation. Total energy-momentum is conserved, with the coefficient Υ extracting energy from ϕ, according to eq. (2.1), and transmitting it to the plasma, according to eq. (2.2). For us the information needed about this dynamics is the retarded propagator of ϕ, determined in ref. [71]. In the local rest frame, it takes the form In order to now compute eq. (3.7), we make use of the real-time formalism of thermal field theory, in the so-called Keldysh (r/a) basis (cf., e.g., ref. [72]). Then the propagator becomes a matrix, where all components are determined by eq. (3.9), (3.11) As for the vertices, they are obtained by substituting ϕ 1 = ϕ r + ϕ a /2 and ϕ 2 = ϕ r − ϕ a /2 in the Lagrangian L(ϕ 1 ) − L(ϕ 2 ). Applying this first to how a metric perturbation h couples 5 The scales at which hydrodynamics applies can be estimated as follows. Consider a sound wave, carrying the energy ω = c s k, where the speed of sound is c s ≃ 1/ √ 3. It is damped by viscous effects, with the rate Γ s ∼ ηk 2 /(e + p), where e + p = T s is the enthalpy and s the entropy density. Shear viscosity is parametrically of order η ∼ T 3 /α 2 . Hydrodynamics is applicable for Γ s ≪ ω, which converts to ω, k ≪ α 2 T . to the energy-momentum tensor T , we obtain (3.12) Repeating the same with the specific structure from eq. (3.8) yields (3.13) We can now consider the retarded correlator of the energy-momentum tensor. Given that the propagator G aa vanishes, the contribution of ϕ originates from − iG R xy;xy k=k e z = T r T a = (ϕ r,x ϕ r,y )(ϕ a,x ϕ r,y + ϕ r,x ϕ a,y ) . (3.14) Carrying out the Wick contractions, going over to momentum space, inserting eqs. (3.10) and (3.11), taking the imaginary part, and symmetrizing in P 1 ↔ P 2 , yields Im G R xy;xy (ω, k) k=k e z = P 1 ,P 2δ × ρ ,x,x (P 1 )ρ ,y,y (P 2 ) + ρ ,y,y (P 1 )ρ ,x,x (P 2 ) + 2ρ ,x,y (P 1 )ρ ,y,x (P 2 ) , (3.15) where Pδ (P) ≡ 1.
As a next step, we may integrate over ω 1 . This can be done with the residue theorem. We first note that the Bose distribution n B (ω 1 ) has poles at ω 1 = iω n ≡ i2πnT , n ∈ . The pole at ω 1 = 0 is lifted by the expression in the square brackets, and the same is true for the pole at ω − ω 1 = 0, from n B (ω − ω 1 ). If we close the contour in the upper half-plane, the remaining contributions are from ω 1 ∈ {iω n , ω + iω n , ±ǫ p + iΥ/2, ±ǫ pk + ω + iΥ/2}, with n ≥ 1.
Even if the residues are readily determined and the corresponding expression could be integrated numerically, it is helpful to put it in a more transparent form, by considering ω, k, Υ ≪ ǫ p ∼ πT . (3.17) The magnitude of ǫ p originates from looking at the domain where the p-integrand becomes suppressed, because of factors ∼ 1/(ǫ p ± iω n ) m or ∼ n B (ǫ p ), with m ∈ AE + . The inequality part of eq. (3.17) is certainly true, given that in the hydrodynamic domain, ω, k < ∼ α 2 T , and that for the maximal temperature, viz. T max ∼ f a /α, Υ max ∼ α 5 T 3 max /f 2 a ∼ α 3 T max . Now, the contributions originating from ω 1 ∈ {iω n , ω + iω n } are strongly suppressed in the limit of eq. (3.17), being parametrically ∼ ωΥ 2 T . The other residues yield much larger contributions, parametrically ∼ ωT 4 /Υ. In the domain of eq. (3.17) we can furthermore approximateǫ Then the leading contributions combine into Im G R xy;xy (ω, k) .

(3.19)
This illustrates a crossover from the regime ω, k ≫ Υ, where the result is suppressed by Υ, to that at ω, k ≪ Υ, where the result is enhanced by 1/Υ. The physical reason for why the gravitational wave production rate at very low frequencies (or the shear viscosity), is inversely proportional to the coupling [68], is that the most weakly interacting particle species display the strongest hydrodynamic fluctuations. The angular integral in eq. (3.19) can be carried out, by going over to spherical coordinates. Inserting 2π 0 dφ cos 2 φ sin 2 φ = π/4 and denoting which at light-cone ω = k has the limiting values 6 The full expression reads F(ω, vk, Υ) = 2(9ω 2 −5v 2 k 2 −3Υ 2 )  .7) can be expressed as For the largest wavelengths, where F can be approximated by the first line of eq. (3.21), this is illustrated numerically in fig. 3(middle). When we go from k ≪ Υ to k ≫ Υ, eqs. (3.21) and (3.22) indicate that the growth of de GW /(dt d ln k) moderates from ∼ k 3 into ∼ k. But the function is still growing, and in fact most of the energy density carried by gravitational waves lies at larger momenta, k ∼ πT . We now turn to how this dominant contribution can be determined.

Boltzmann domain
At larger momenta, elementary particle excitations can be resolved, and we need to consider the microscopic form of the energy-momentum tensor. Omitting trace parts, which drop out when projected with eq. (3.5), the fields appearing in eq. (1.1) give the contribution to the traceless part. These components couple to the propagating part of the graviton field (h). We are interested in the contribution to graviton production that involves one appearance of the vertex in eq. (1.1), as the processes without this vertex were already considered in ref. [65]. Various processes are depicted in fig. 2. (The 2 → 1 channel ϕϕ → h is not kinematically allowed with on-shell particle states.) A way to represent and evaluate the rates of the reactions in fig. 2 has been presented in ref. [74]. In the following, we adopt its methods and notation. The procedure starts by considering the processes in fig. 2(b), which are not kinematically allowed, but have a simple would-be algebraic structure, as the non-equilibrium particle and the plasma particles are on different sides of the reaction. This contribution is represented as where scat 1 → 3 is a phase-space average, 7 and g 1 , g 3 label two (identical) gauge bosons.
The dynamical information concerning the production process enters through the function Θ in eq. (3.24), which may be referred to as "matrix element squared". More precisely, if we couple T µν to polarization vectorsĥ µν of a would-be source (L ⊃ĥ µν T µν ), and replace the sum over the polarizations through eq. (3.5), viz.
where s 1 , s 3 label the helicities of the final-state gauge bosons, and 1 2 accounts for the gauge bosons being identical (with this factor we can integrate over the full phase space without the danger of overcounting). 8 The tensor Ä αβ;µν , defined in eq. (3.5), contains a quadratic appearance of the projector Ã T . After contractions with the metric tensor, linear and zeroth order terms in Ã T can appear as well (cf. eq. (3.6)). Redundancies can be eliminated by making use of the fact that Ã T is with analogous relations obtained through the relabellings 1 ↔ 3 and 2 ↔ 3.
A further ingredient of the computation is the Levi-Civita tensor ǫ µνρσ from eq. (1.1). In 7 In ref. [74], the non-equilibrium particle was defined to be the initial state, i.e. time was running in the opposite direction, which explains the reference to a 1 → 3 process. In the particle production language of fig. 2(b), it is more intuitive to depict the non-equilibrium particle as a final state, yielding a 3 → 1 reaction. 8 Alternatively, Θ can be extracted from the definition in eq. (3.24), i.e. by computing the retarded 2-point correlator of the energy-momentum tensor and taking its imaginary part, as illustrated in fig. 2(a). Either way, it is important to crosscheck the gauge independence of the result.
The numerical integration of eq. (3.34) can be carried out by modifying the algorithm provided in ref. [74]. We note that the production peaks at temperatures T ≫ m, and thus appearances of m can be omitted in practice. In this limit there are no poles in eq. (3.33), whereby the virtual corrections that were discussed in ref. [74] play no role.

Numerical estimates
The purpose of this section is to summarize the parametric forms of the results that were obtained in secs. 3.2 and 3.3, and to illustrate the corresponding prefactors numerically.
If we set k ∼ πT , where the production rate proportional to k 3 n B (k) peaks (cf. eq. (3.34)), then de GW dt d ln k ∼ αT 7 /m 2 Pl for the Standard Model contribution to gravitational wave production [64]. The result in eq. (3.33) contains the prefactor α 2 /f 2 a , implying that axion-like inflation leads to the additional contribution de GW dt d ln k ∼ α 2 T 9 /(m 2 Pl f 2 a ). The numerical coefficients associated with these parametric behaviours are illustrated in figs. 3(left) and (right), respectively, where the running of α has been taken into account, assuming N c = 3 and a QCD-like initial value at low energies (numerically, α ∼ 0.015 in the temperature range shown). In addition we have plotted the estimate in the extreme hydrodynamic domain, from eq. (3.22), in fig. 3(middle), even though we do not think that this result has significance for our main conclusions. The hydrodynamic prediction scales as de GW dt d ln k ∼ f 2 a k 3 T 2 /(α 5 m 2 Pl ), assuming Υ ≃ κα 5 T 3 /f 2 a , with a numerical coefficient κ ≃ 100 [75]. Given the similar shapes but different normalizations in figs. 3(left) and (right), we may expect the axion contribution to gravitational wave production to exceed the Standard Model one at T > ∼ 10 3 f a . However, the numerical solutions in ref. [55] only reached T max ∼ 200f a . It thus appears that the axion contribution does not exceed the Standard Model one. Furthermore, for the Standard Model contribution, T max = 2 × 10 17 GeV increases the massless degrees of freedom only by ∆N eff ≈ 10 −3 [65], which is very demanding to observe [66]. Therefore, reheating through a coupling between an axion-like inflaton and non-Abelian gauge fields is not excluded at present, and represents a viable scenario for the foreseeable future.

Summary and conclusions
Depending on its magnitude, the gravitational wave background produced by a reheating process [1][2][3][4] can lead to one of two possible consequences. If the background is substantial, this would be exciting as a motivation for possible future experiments (cf., e.g., ref. [6]). If it is moderate, we can be confident that the model is not already excluded, as could happen in the case of an axion-like inflaton coupled to Abelian gauge fields (cf., e.g., ref. [10]). For non-Abelian reheating after axion-like inflation, the magnitude of the gravitational wave production rate depends on the parameters α, f a , m and T max (cf. eq. (1.1)). The dependence on m is small enough to be negligible in practice, provided that m ≪ πT max , as is the case towards the end of the reheating period [55]. Within the setup of eq. (1.1), and in the domain where most of the energy density lies, the dependence on f a is given by the power law 1/f 2 a , and the dependence on T max by dimensional analysis. Therefore the main task has been to sort out the dependence on α, and to determine the associated prefactor.
Our numerical results are illustrated in fig. 3. Parametrically, the axion contribution exceeds the Standard Model for T > f a / √ α. Numerically, this has turned into T > 10 3 f a , which is unlikely to be reached according to ref. [55]. A main reason for the numerical suppression is the small factor c χ in eq. (1.1), which appears quadratically in the production rate. If the production rate does not exceed the Standard Model one, it is not strongly constrained in the temperature range that is associated with normal inflationary scenarios, T max ≪ 10 17 GeV, given that in this range the Standard Model contribution increases the energy density as parametrized by massless degrees of freedom only by ∆N eff ≪ 10 −3 [65]. It may be wondered why the non-Abelian case differs so notably from the Abelian one, where an efficient tachyonic instability has been claimed to convert a significant fraction of energy density to gravitational waves. The reason is that backreaction effects lead arguably to rapid thermalization. In a thermal system, as we have assumed to be the case, tensor modes are excited only through interactions, whereby their production is suppressed by α 2 .
To sharpen our conclusions, it would be nice to fix T max in terms of f a , m and α such that inflationary predictions are in line with observation. This requires going beyond the universal eq. (1.3), by defining V (ϕ) away from the minimum. Ultimately, it would also be great to probe non-equilibrium effects, and to employ a UV complete description, as reheating easily takes us to a domain where the non-renormalizable operator in eq. (1.1) is having a substantial influence. We hope to return to some of these issues in the future.