The Cosmological Heavy Ion Collider: Fast Thermalization after Cosmic Inflation

Heavy-ion colliders have revealed the process of"fast thermalization". This experimental break-through has led to new theoretical tools to study the thermalization process at both weak and strong coupling. We apply this to the reheating epoch of inflationary cosmology, and the formation of a cosmological quark gluon plasma (QGP). We compute the thermalization time of the QGP at reheating, and find it is determined by the energy scale of inflation and the shear viscosity to entropy ratio $\eta/s$; or equivalently, the tensor-to-scalar ratio and the strong coupling constant at the epoch of thermalization. Thermalization is achieved near-instantaneously in low-scale inflation and in strongly coupled systems, and takes less than an e-fold of expansion for weakly-coupled systems or after high-scale inflation. We then consider the potential for observing this process: we demonstrate that the cosmic microwave background is largely insensitive, and the shift in $n_s$ and $r$ is well below the sensitivity of CMB S4 and Simons Observatory. We also find a stochastic background of gravitational waves at frequencies accessible by interferometers, albeit at a level unobservable by even next generation experiments.


I. INTRODUCTION
Thermalization in quantum field theory has a long history (see e.g. [1]), and, in particular, thermalization in non-Abelian gauge theory has been a driving research question in Quantum Chromo-Dynamics (QCD) for decades. There are two places in the history of the universe when a non-Abelian gauge theory at high density has thermalized: at the present day at the Relativistic Heavy Ion Collider (RHIC) [2,3], and in the early universe, at the beginning of the radiation dominated era of standard cosmology.
A primary question that has arisen in the QCD community has been the mechanism behind the rapid approach to thermal equilibrium, so-called fast thermalization, that has been observed at RHIC [2,3]. This has lead to radically different approaches, both at weak coupling [4][5][6][7][8], where a particle or hydrodynamics description can be applied, and at strong coupling [9][10][11][12], where thermalization of the boundary gauge theory corresponds to black hole formation in the bulk spacetime [13].
Parallel to this, a central question in early universe cosmology, which has risen in prominence with the precision of cosmic microwave background experiments [14], is the duration of the period of reheating that necessarily follows cosmic inflation [15][16][17][18][19]; for reviews see e.g. [20][21][22][23]. The uncertainty in the duration of reheating leads to an uncertainty in two of the hallmark predictions of inflationary cosmology, namely the scalar spectral index n s and the tensor-to-scalar ratio r [24][25][26] of the primordial perturbations. Narrowing and quantifying these uncertainties has been a primary aim of the reheating literature [27][28][29][30][31][32]; and has spurred on the detailed on numerical calculations [33][34][35][36][37][38][39][40]. This issue aside, the self-consistency of the inflationary universe scenario requires the universe reach a state of near thermal equilibrium at a temperature of or above ∼MeV [41,42], else a signature success of standard cos-mology, namely big bang nucleosynthesis, would be forgone. This imposes a constraint that the universe thermalize before the energy density of the universe redshifts below MeV 4 , providing a baseline test of self-consistency, which one may intuit to be satisfied, but ultimately must be derived or numerically demonstrated.
It is conventional in the cosmology community to define the transition from reheating to standard cosmology as the moment when the universe is dominated by relativistic particles, or alternatively, when the equation of state of the universe is sufficiently close to w = 1/3, see e.g. [30]. However, by ignoring the thermalization process, this abandons the rich underlying physics, and potentially rich phenomenology, that can be accessed by applying the machinery developed for thermalization at heavy ion colliders.
In this work we will study the thermalization process during the reheating epoch after cosmic inflation, which has has been neglected in all but a few previous works [43][44][45][46][47][48][49] 1 . We apply results from the thermalization literature, namely the relation between thermalization time, temperature, and shear-viscosity-to-entropy ratio η/s, and the theoretical predictions for η/s at strong and weak coupling. From this we will demonstrate that the post-inflationary universe, like RHIC, undergoes fast thermalization. We will compute the imprint on observables, namely on n s , r, and the production of a stochastic background of gravitational waves; each far below experimental sensitivities. This sets the stage for further applications of modern QCD techniques to the very early universe.

II. THERMALIZATION AT THE COSMOLOGICAL HEAVY ION COLLIDER
There are many theoretical approaches to the quark gluon plasma, with different regimes of applicability, depending on the coupling and occupation number of the plasma. At weak coupling and perturbative occupation numbers n 1/α s , one can derive an effective kinetic theory governed by a Boltzmann equation [7,53,54]. At weak coupling but non-perturbative occupation numbers n k ∼ 1/α s , one can take a different approach and model the system as weakly coupled fluctuations of a condensate [55], a so-called Color-Glass Condensate. In this same regime, lattice gauge theory studies suggest that there may occur a process of Turbulent Thermalization [56]. For couplings that are large but still smaller than the non-perturbative value 1/α s , the system is effectively classical and perturbative calculations can be made in classical Yang-Mills theory [57].
A general result of these approaches is that the thermalization time and temperature of the quark gluon plasma are both related to shear viscosity to entropy ratio. The three quantities are intertwined as [58,59], where η/s is the shear viscosity to entropy ratio. This is in turn determined by the interaction strength of the quantum field theory. At strong coupling, the only analytical tool available is holography, namely the Anti-de-Sitter/Conformal Field Theory correspondence (AdS/CFT) [60]. In this context the thermalization of a non-Abelian plasma can be described as black hole formation in an anti-de-Sitter spacetime [13]. A remarkable result from this field is a lower bound on the shear-viscosity-to-entropy ratio [12,61,62], This is computed in the gravity dual as the zerofrequency limit of the graviton absorption cross-section of black three-branes. This lower bound can be violated by higher-derivative terms in the gravitational action [63], or by the inclusion of bulk viscosity. For the purposes of this work, we consider the above to a lower bound.
In what follows, we will extend (1) to an FRW universe. This will allow us to calculate the number of e-folds of the thermalization process after inflation. Importantly, we note that (1) is independent of the expansion of space: solving self-consistently for the redshifting of energy and the thermalization time alters the scaling with η/s of both t th and T th , leaving (1) unchanged.
The expansion of space is related to the thermalization temperature via the Friedmann equation. That is, the thermalization temperature is related to the Hubble parameter at the end of thermalization as Meanwhile, the time evolution of the Hubble constant, for an FRW universe with equation of state w > − 1 3 , is given by where the latter equality applies for a constant w, ∆t ≡ t − t i , andŵ is defined aŝ Thus the Hubble constant at the end of thermalization can also be written as Combining equations (3) and (6) with (1) uniquely fixes t th and T th in terms of η/s, w, and the initial value of the Hubble constant H i . The Hubble constant at the culmination of thermalization is, From this we see that the equation of state w acts as a simple rescaling of η/s. One can also see that H th → H i at strong coupling, i.e. the limit η/s → 0, while at weak coupling, η/s → ∞, H th tends to zero: This is a simple consequence of the redshifting of H over the course of thermalization. This in turn determines the thermalization time t th via (6), and the thermalization temperature T th via (3). The latter is given by, We can relate this to number of e-folds of thermalization N th , using dN = Hdt, and assuming w a constant, which, inverted for ∆t, gives, Plugging this in equation (1) with ∆t = t th gives, with T th given in terms of η/s via equation (9), andŵ in range 1 ≤ŵ ≤ 2, as defined in equation (5). These results allow for N th and T th to be specified uniquely via η/s and the energy scale of inflation H i . This can be directly related to the scalar-to-tensor ratio [64], Hence the thermalization process is itself sensitive to the tensor-to-scalar ratio.

III. THERMALIZATION TIME AT STRONG AND WEAK COUPLING
We are now in a position to make quantitative statements as to the thermalization time of the universe following inflation. We note that the energy scale of the thermalization process is not a priori the inflationary energy scale, but instead is the energy of the decay products predominantly produced at (p)reheating. This is sensitive to the (p)reheating mechanism; whether it is broad or narrow resonance, or perturbative decays, and could be in IR modes or UV modes relative to H.
To compute the thermalization time, one needs the shear-viscosity-to-entropy ratio η/s. This is calculable via different methods, and at weak coupling we have [5,58] η s | αs 1 = 34.784 where λ is the t 'Hooft coupling, λ = g 2 N c for N c colors of gluon, which is related to the strong coupling constant α s ≡ g 2 4π by λ = 4πN c α s . As a fiducial example, for λ = 0.1 this evaluates to η/s 10 3 . On the other hand, at strong coupling, η/s is bounded: This bound arises from fundamental considerations in the AdS/CFT correspondence. From this one can immediately deduce the AdS/CFT prediction for the thermalization time at reheating. It follows from (1), as This is much less than 1 even for high scale inflation.
Thus we see that thermalization at strong coupling is always cosmologically fast, lasting for much less than one e-fold of expansion. Moreover, the thermalization temperature is, To leading order in H i /m pl this is the standard result of instant reheating. It follows that reheating at strong coupling is in no danger from a low thermalization temperature.
On the other hand, to make predictions at weak coupling one must understand in detail the renormalization group flow of the strong coupling constant. In standard QCD, from the β-function for SU (3) with N f = 6, the energy-scale dependence of α s is given by, with M Z = 91.2 GeV. At very high energies, e.g. µ = m pl , this is only as small as 0.0198. However, the highenergy behavior of α s is sensitive to new particles that emerge at high energies. The inclusion of new gauge bosons accelerates asympotic freedom, as the β-function scales with N c . Additionally, at energies near the Planck one expects quantum gravity effects to become important, which may radically alter the evolution of the effective gauge coupling. With these possibilities in mind, here we consider the number of e-folds of thermalization at strong and weak coupling. See Figure 1. For simplicity here we identify the energy of scatterings µ with the energy scale of inflation, and we fix w = 1/3. From this plot one can appreciate that the thermalization time is sensitive to the energy scale of inflation, and thermalization generically takes longer at weak coupling and/or in large field inflation, taking less than an e-fold of expansion. On the other hand, for both low-scale inflation and moderate coupling, thermalization takes much less than an e-fold.

IV. TOWARDS A DETECTION
The analysis thus far is an important consistency check on the inflationary universe scenario. However, ideally we would have an observable associated with the detailed physics that occurs during the process of thermalization.
The canonical observational probe of reheating is a shift in n s and r [24-27, 29, 30, 65]. This occurs solely due to the modified expansion history of the universe, and not the detailed microphysics. A period of expansion lasting ∆N e-folds with w ∆ is completely degenerate with a shift in the number of e-folds before the end of inflation that the CMB pivot scale exited the horizon, by an amount [24,26], This completely encodes the reheating dependence of inflationary predictions. In particular, the Hubble constant at the moment a fluctuation of wavenumber k exits the horizon H k , from which one can derive n s and r, is given by [26], where N k is the number of e-folds before the pivot scale exits the horizon, and for simplicity the numerical factor of 61.6 accounts for the spin-degrees of freedom (fixed in the above to be 100) and the CMB pivot scale, taken to be .05 Mpc −1 . This determines n s and r in a model-dependent and -indepedent fashion respectively. For polynomial potentials, V = m 4−p φ p , one finds [26] As a numerical example, consider m 2 φ 2 inflation, which is not a good fit to CMB data, but as a toy model is suitable for our purposes. The modification to the spectral index from a non-thermal phase is [26], Indeed we see that the modification to the spectral index vanishes if w = 1/3, since in this case, the expansion history of the universe from standard cosmology. At energies far above the QCD scale, the quark gluon plasma is approximately conformal. It can be approximated as 'pure glue', in which case it is easy to check that w = 1/3 at all times. Namely, the Maxwell stress tensor, is traceless: T µ µ = 0. This in turn implies that ρ − 3p = 0 and hence w = 1/3. Thus the thermalization process at sufficiently high energies will have no impact on the spectral index of primordial perturbations.
At lower energies, the situation is more subtle. We again emphasize that the energy scale of the thermalization process, µ, is not the inflationary energy scale, but instead the energy of the decay products predominantly produced at (p)reheating. This energy scale could be far below the inflationary scale, e.g. super-horizon modes produced via broad parametric resonance, or far above the energy scale, produced by perturbative decays of the inflaton.
The equation of state of the quark gluon plasma can be computed using lattice QCD, and at temperatures around 10 2 -10 3 MeV can differ substantially from 1/3, see e.g. [66]. For the present analysis it is sufficient to consider 0 < w < 1/3; the former bound arising from the limit of a collection of non-relativistic particles and the latter the high-temperature limit. In this case, using the previous result that N th 1, we find, This shifts n s by at most the level of 1 part in 10 4 ; e.g. for N th = 1 and w = 0.1, the shift in n s is δn s = 3.9 × 10 −5 . This is two orders of magnitude below the sensitivity of CMB S4 [67]. However, this is not the whole story: The non-Abelian gauge fields themselves will source gravitational waves. Gravitational waves are in general sourced by anisotropic stress, as encoded in the equation of motion, where Π ij is spatial off-diagonal components of the stresstensor, or more formally, the transverse traceless projection of the spatial stress tensor. The spatial stress tensor is given in full generality as, One can easily see that it is η which will play a key role in generating gravitational waves. The sourcing of gravitational waves from a non-Abelian plasma was studied in detail in [51]. The sourcing of gravitational wave spectra is, where the TT ij are polarization tensors. Evaluating the sum over polarizations, and choosing k = k e 3 , leads to Pl X e ik(t−z) T 12 (0) T 12 (X ) . (27) Utilizing the stress tensor, one can simplify this on large scales, which gives the power in gravitational waves on large scales as, We can relate this to η/s using s = g * T 3 in a radiation dominated universe, with g * the number of spin degrees of freedom, as well as H 2 T 4 /3m pl 2 , to find, The net production of gravitational waves at thermalization is then given by an integral over the thermalization time period.
As an estimate, guided by the results of the previous section, we consider the limit of fast thermalization, N th 1. We change the time coordinate to the number of e-folds, and approximate thermalization as an instantaneous process. The produced gravitational waves are, The signal observed today is related to the above by a simple redshifting by a(t) −4 : The redshift factor can itself be re-expressed in terms of the thermalization temperature, by relating to known quantities at matter-radiation equality; i.e. a th /a 0 = (a th /a eq )(a eq /a 0 ) where a eq /a 0 = 1/(1 + z eq ) 1/3000 and a th /a eq ∼ T eq /T th , with T eq ∼ 3 eV and T th ∼ H th m pl . The resulting present day spectrum of gravitational waves is given by, From this one can appreciate that the redshifting of gravitational waves has flipped the dependence on H th : highscale thermalization will lead to a lower amount of gravitational waves in the present universe. The gravitational wave spectrum is heavily blue tilted, and thus suppressed on large scales. This makes the effect unimportant for CMB B-modes, but opens the possibility for observing it with interferometer experiments. The quantity of interest is, see e.g. [64], where ρ c is the critical energy density of the universe today. We translate from k-mode to frequency space via f = k/2π, following the unit conversion k/Mpc −1 = 6.5× 10 14 f /Hz. From this we find, Ω GW (f ) = 5.6 × 10 −108 g * η s GeV The extremely small prefactor prevents this signal from reaching any appreciable level. This is simply a reflection of the hierarchy between scales in the problem, e.g. the ratio m pl /Mpc −1 10 56 .
The current LIGO constraint [68] on a stochastic gravitational wave background is given by, Even in the most optimistic case, saturating the BBN bound H th ∼ MeV, and assuming (η/s) ∼ 10 3 , g * ∼ 10 2 , the gravitational wave signal from thermalization is well below the observational bound.

V. DISCUSSION
In this work we have applied and generalized results from the thermalization of non-Abelian gauge theories, originally developed to understand the formation of quark gluon plasma at heavy ion colliders, to the epoch of reheating after cosmic inflation. We derived a relation between the number of e-folds of the thermalization process, the shear-viscosity-to-entropy ratio η/s, and the scalarto-tensor ratio r. We find that thermalization generically completes within much less than a single e-fold of expansion, and the observable effects, namely shifts in n s and r and a stochastic background of gravitational waves, are well below observational sensitivities.
It is imperative to study further possibilities for observing the microphysics of the formation of the primordial quark gluon plasma. For example, it is known that thermalization of the quark gluon plasma can be delayed by a long-lived phase nearby a non-thermal fixed point [69]; it would be interesting to consider the cosmological realization of this, and the impact on n s and r. It is also possible that there could be a phase transition during the thermalization process itself, which would further source gravitational waves. We leave these interesting possibilities for future work.
Finally, we have not endeavored to consider the detailed interplay of thermalization and preheating. Holographic methods have been applied to preheating in [70,71]. It would be interesting to consider the gravitational wave signal once these effects are included.