Gravitational wave signatures of post-fragmentation reheating

After cosmic inflation, coherent oscillations of the inflaton field about a monomial potential $V(\phi)\sim \phi^k$ result in an expansion phase characterized by a stiff equation-of-state $w\simeq(k-2)/(k+2)$. Sourced by the oscillating inflaton condensate, parametric (self)resonant effects can induce the exponential growth of inhomogeneities eventually backreacting and leading to the fragmentation of the condensate. In this work, we investigate realizations of inflation giving rise to such dynamics, assuming an inflaton weakly coupled to its decay products. As a result, the transition to a radiation-dominated universe, i.e. reheating, occurs after fragmentation. We estimate the consequences on the production of gravitational waves by computing the contribution induced by the stiff equation-of-state era in addition to the signal generated by the fragmentation process for $k=4,6,8,10$. We find that the signal generated during the fragmentation process gives a larger contribution than the one induced by the stiff equation-of-state era in given frequency ranges for all values of $k$. Our results are independent of the reheating temperature provided that reheating is achieved posterior to fragmentation. Our work shows that the dynamics of such weakly-coupled inflaton scenario can actually result in characteristic gravitational wave spectra with frequencies from Hz to GHz, in the reach of future gravitational wave observatories, in addition to the complementarity between upcoming detectors in discriminating (post)inflation scenarios. We advocate the need of developing high-frequency gravitational wave detectors to gain insight into the dynamics of inflation and reheating.


Introduction
The inflationary paradigm, which predicts the existence of an early epoch of accelerated expansion of the universe, provides one of the most compelling solutions to the initial conditions problems of the standard Big Bang cosmological model. 1 In its arguably simplest realization, the quasiexponential growth of the universe is driven by the slow roll of a spatially homogeneous scalar field (the inflaton) over a sufficiently flat potential.The excitation of the quantum fluctuations of this field and the spacetime metric would source an almost scale-invariant power spectrum of adiabatic scalar fluctuations, which is in remarkable agreement with Cosmic Microwave Background (CMB) and Large Scale Structure observations.Nevertheless, what is considered the smoking gun signal for inflation, a primordial spectrum of parity odd (B-mode) polarization fluctuations in the CMB, sourced by tensor modes excited during inflation, has so far escaped the capabilities of existing CMB observatories.Currently, the bound on the ratio of the tensor and scalar spectra at the pivot scale k * = 0.05 Mpc −1 is r < 0.036, at the 95% CL [5].The non-observation of such primordial tensor modes seems to favor models with concave, asymptotically flat potentials [6], among which the Starobinsky [7], Higgs [8], and α-attractor models [9][10][11][12][13] are theoretically well motivated examples.
The superhorizon stretching of vacuum tensor modes during the slow roll of the inflaton is however not the only source of gravitational waves predicted by inflation.The epoch of accelerated expansion must eventually come to an end, after which the energy density of the inflaton is to be transferred to the Standard Model (and possibly dark matter) degrees of freedom, during what is known as the reheating epoch.In the slow roll picture, this "graceful exit" mechanism is built in, where the coherent oscillations of the inflaton field about the minimum of its potential source particle production via the coupling of the inflaton to light degrees of freedom [14][15][16][17][18][19][20][21].Although for small couplings this dissipation process can be modeled perturbatively, for large couplings collective bosonic enhancement (or fermionic blocking) effects require the use of non-perturbative techniques.For bosonic fields, this preheating leads to the exponential growth of the fluctuations of the decay products, due to the non-adiabatic change of their effective masses, manifested as parametric resonance [15,[22][23][24][25][26].If this growth is sustained, the coherent inflaton condensate will eventually be fragmented, due to the large gradients sourced by mode-mode couplings (backreaction) [18,[27][28][29][30][31].Importantly, these large amplitude gradients can in turn efficiently source a spectrum of high frequency gravitational waves [32][33][34][35][36][37][38][39][40][41][42][43][44].On top of this efficient mechanism, gravitational waves can also be sourced directly by the oscillating condensate [45], or from graviton bremsstrahlung during reheating, from inflaton decay or inflaton annihilation [46][47][48][49][50]. Additionally, if certain conditions are met during the preheating epoch, long-lived solitonic-like configurations of the inflaton field called oscillons can be formed through inflaton fragmentation during reheating and can also yield a significant gravitational wave signal [51][52][53].Hence, the primordial stochastic gravitational wave background is a potential powerful probe of post-inflationary dynamics [44].
In this work we revisit the production of gravitational waves due to the resonant growth of inflaton fluctuations due to its self-resonance.Namely we consider models in which the potential of the inflaton ϕ can be approximated as V (ϕ) ∝ ϕ k (with k even) near its minimum.For k > 2, the self-interaction of the inflaton will eventually fragment the homogeneous inflaton condensate after the end of inflation, unless reheating can be completed promptly.After fragmentation, the universe is dominated by the relativistic free inflatons, which must subsequently decay into the Standard Model primordial plasma.This implies that the equation of state of the universe transitions as w → −1/3 → k−2 k+2 → 1/3 from the end of inflation to the end of backreaction.Here we will work under this late reheating assumption, taking into account the evolution of the equation of state parameter and the post-fragmentation completion of reheating, recently studied in [54][55][56], to determine the present-day amplitude of the stochastic gravitational wave background sourced by the inflaton self-resonance.Moreover, although the super-horizon fluctuations sourced during inflation are generically impervious to the nonlinear preheating dynamics [57][58][59], modes that reenter the horizon prior to the end of backreaction will experience a relative enhancement due to the stiffer-than-radiation equation of state parameter [60][61][62][63][64][65][66], which we also quantify here.By means of numerical computation and lattice techniques, we keep track of the time-dependence of the equationof-state parameter from a stage deep into the inflation era until the completion of the fragmentation process.We use both analytical and lattice techniques to estimate the total resulting gravitational wave signal that we confront to prospects from potential future gravitational wave observatories.
This paper is organized as follows: In Section 2 we review the transition from the inflation slow roll phase (Sec.2.1) to the coherent oscillation of the homogeneous field (Sec.2.2), the growth of inhomogeneities sourced by parametric resonances (Sec.2.3), the numerical simulation of backreaction (Secs.2.4 and 2.5), until the eventual completion of reheating via inflaton decays (Sec.2.6).Section 3 contains the main results of this work, and is devoted to the computation of the primordial gravitational wave signal.Sec.3.1 presents the resulting energy densities and spectra originated from inflaton self-fragmentation, for a selection of inflaton potentials.In Sec.3.2 we address the stochastic spectrum generated from vacuum fluctuation growth during inflation, accounting for the effect of the time-dependent equation of state during reheating.Secs.3.3 and 3.4 present a brief discussion of additional contributions to the tensor spectrum, and the comparison with observational prospects, respectively.Our summary and conclusions are presented in Section 4.
2 Post-fragmentation reheating: from inflation to a radiation-dominated phase

Inflation
We consider a phase of cosmic inflation driven by a single inflaton field denoted by ϕ, minimally coupled to gravity, with a corresponding action of the form where GeV is the reduced Planck mass, g is the metric determinant and R the Ricci scalar.L int represents interaction terms between the inflaton and additional species, responsible for reheating.For the inflation model we consider the α-attractor T-models described by the potential [12] V (ϕ) = λM 4 which are compatible with the current CMB constraints [5,6], and which appear in the context of no-scale supergravity [21].The exponent k controls the behavior of the potential close to the minimum Along with the exponent k, the inflaton potential can be fully determined by fixing λ, which controls the potential at large field values.One can reproduce the central value of the amplitude of the scalar power spectrum, A S = e 3.044 /10 10 ≃ 2.1 × 10 −9 as measured by Planck [67], for the choice of parameter where N * is the number of e-folds between the horizon-crossing of the CMB fiducial scale k * = 0.05 Mpc −1 (corresponding to k * = a * H * ), and the end of inflation.It can be related to the field values evaluated at these times, ϕ * and ϕ end respectively, as In the following, for definiteness, we chose N * = 55. 2ccelerated expansion is maintained until the end of inflation when ä = 0, where dots represent differentiation with respect to cosmic time and a is the scale factor, or equivalently when the derivative of the field is φend = − V (ϕ end ).The field value ϕ end can be approximated by [20,68] This corresponds to an energy density at the end of inflation where H ≡ ȧ/a is the Hubble parameter and k = 4.This energy density only mildly depends on k, with ρ end ≃ (4.9 × 10 15 GeV) 4 for k = 10.

Coherent inflaton oscillations
After the end of inflation the inflaton oscillates about the minimum of its potential.The dynamics can be described by the potential (2.3) which is a good approximation to Eq. (2.2) already after one oscillation.Variation of the action (2.1) with respect to the homogeneous inflaton field and metric yields the equation of motion for the inflaton field and the Friedmann equation Multiplying (2.8) by ϕ, and averaging over oscillations, it is straightforward to obtain the following virial-like relation for the mean kinetic energy, with V ϕ ≡ ∂V /∂ϕ.This allows us to express the oscillation-averaged energy and pressure densities as [20] (2.12) The (oscillation-averaged) Equation-Of-State (EOS) parameter w ≡ P ϕ /ρ ϕ can thus be uniquely determined in terms of the exponent k through This corresponds to a redshift of the energy density of ϕ of the form ρ ϕ ≃ ρ end (a/a end ) − 6k k+2 .During its oscillations about the minimum, the inflaton field can be expanded as where ϕ 0 (t) ∝ ρ 1/k ϕ ∝ a −6/(k+2) is a decaying envelope and P(t) is periodic function encapsulating the fast inflaton oscillations.To facilitate a k-universal treatment, it is convenient to define the rescaled periodic function P ≡ (ϕ end /M P )P and to introduce a new k-dependent time variable related to cosmic time t via dτ α ≡ (a/a end ) −3(k−2)/k+2) √ λM P dt such that the EOM (2.8) can be recast as Solutions to this equation can be found in implicit form as where 2 F 1 are hypergeometric functions.The oscillating functions are represented in Fig. 1 for k = 2, 4, 6, 8, 10 as a function of the α-time over one period of oscillation T α where the initial time τ α = 0 as been chosen when P is zero.Starting from a sine solution for k = 2, by increasing k higher harmonics of the periodic function P contribute more significantly, driving the solution towards a "triangular" periodic function.

Growth of inhomogeneities
Inhomogeneities in the post-inflation era are characterized by space and time dependent inflaton fluctuations δϕ(t, x). 3 By initially neglecting couplings of the inflaton to other degrees of freedom (corresponding to the term L int in (2.1)), variation of the action leads to the equation of motion for the inflaton perturbation One can rewrite the equation of motion in terms of the canonically normalized inflaton fluctuation X = a δϕ expanded in Fourier space where p is the comoving momentum.The creation and annihilation operators âp and â † p satisfy the canonical commutation relations In addition, the Wronskian constraint X p X * ′ p − X * p X ′ p = i is imposed to fulfill the canonical commutation relations between X and its canonically conjugated momentum.Here ′ ≡ d/ dτ represents differentiation with respect to conformal time, dτ = dt/a.Substitution of δϕ = X/a into (2.18)gives an equation describing the evolution of mode functions where the effective frequency reads To keep track of the growth of inhomogeneities during the post-inflation epoch, one solves the mode equation by assuming at some initial time τ 0 the positive-frequency Bunch-Davies vacuum when the mode is deep inside the Hubble radius (i.e.p > a(τ 0 )H(τ 0 )).The mode equation (2.20) features parametric instabilities driven by the oscillating term in the frequency (2.21) which can be made more explicit in terms of the periodic function P were we used the second Friedmann equation a ′′ /a = (1 − 3w)H/2 with H ≡ aH the conformal Hubble parameter.Importantly, even for a weakly coupled theory with λ ≪ 1, necessary for CMB compatibility c.f. (2.4), parametric resonances ensure a significant growth of inhomogeneities if the oscillating term is active for a substantial amount of time.According to Floquet's theorem, the solution to the mode equation can be expressed as [75] X p (τ ) = e µpτ g 1 (τ ) + e −µpτ g 2 (τ ) , where g 1,2 are periodic functions and µ p are complex numbers called the Floquet exponent.Eq. (2.20) can be solved by identifying instability bands corresponding to Re(µ p ) > 0 signifying exponentially growing modes.We refer the reader to [54,76] for further details on the Floquet analysis.In general, instability bands are time-dependent through the scale factor dependence of the second terms in brakets in Eq. (2.23).However, for k = 4 all scale factor dependence cancels out as a consequence of conformal invariance of the quartic potential, as discussed in Ref. [54].In our case, exponential growth will occur once modes cross the first narrow instability band [76].One can estimate the typical scales subject to parametric resonances p/(a end √ λM P ) ≃ O(1), 1.1×10 −2 , 1.6× 10 −5 , 7.6 × 10 −9 respectively for k = 4, 6, 8, 10 [55].The growth of inhomogeneities eventually becomes large enough to trigger significant non-linearities and mode-mode couplings whose proper treatment requires numerical integration.

Simulating fragmentation
In order to account for non-linearities induced by the growth of inhomogeneities, we perform lattice simulations with the use of CosmoLattice [77,78].We first solve the background inflaton equation of motion during inflation and determine the value of various background parameters at the end of inflation.We use these values as initial conditions with CosmoLattice to simulate the post-inflation dynamics.The inflaton field is treated as a classical field on the lattice and initial conditions for inhomegeneities are randomly assigned to reproduce the statistical properties of the Bunch-Davies vacuum.The time variable used for the simulation is the α−time τ α (see (2.15)) and the lattice parameters are chosen appropriately to fully capture parametric resonances.The total energy density of the inflaton is computed on the lattice from the spatially-averaged energy-momentum tensor of ϕ(x, τ ), denoted with an over-bar, Following [54,55], one can define a condensate (i.e.homogeneous) component of the energy density as the sum of kinetic and potential energies of the spatially-averaged inflaton field [30] Inhomogeneities are characterized by the comoving occupation number of ϕ, which in terms of mode functions and their time-derivatives is given by the adiabatic invariant [26] By integrating this quantity, one can infer the number density and the energy density of the inflaton fluctuations expressed as the UV-regular forms [26,30] ) (2.29) Notably, the later expression corresponds to ρ δϕ = ρ ϕ − ρ φ [54,56].Each contribution to the energydensity resulting from our simulations is represented in Fig. 2 as a function of the scale factor for k = 4, 6, 8, 10.The corresponding oscillation-averaged EOS parameter is in turn shown in the bottom panels.For all cases, initially the EOS is constant and follows the expected scaling relation of Eq. (2.13) as the universe is dominated by the coherently oscillating condensate.The growth of inhomogeneities can be observed from these plots as a sharp increase in ρ δϕ until backreaction is reached and the inflaton condensate is fragmented.We estimate that fragmentation is reached at a frag /a end ≃ 180, 4.5 × 10 4 , 6 × 10 6 , 7 × 10 8 for k = 4, 6, 8, 10 consistently with Ref. [54,55,76,79].Following this event, the condensate component redshifts away and becomes subdominant as the universe become filled with a collection of relativistic inflaton quanta dominating the energy budget.The EOS relaxes progressively to 1/3 over typical time scales of O(2−3)(a frag /a end ), slightly dependent on the value of k.

Transition to an inflaton-quanta dominated phase
As it can be appreciated in Fig. 2, the transition from the condensate EOS to a radiation dominated universe is not instantaneous.For this reason, it is convenient to introduce the scale factor a RD > a frag at which the universe is almost exclusively composed of a population of fragmented relativistic inflaton particles.That is, k−2 k+2 > w > 1/3 for a frag < a < a RD , and w ≃ 1/3 for a > a RD .In practice we take a RD as the scale factor evaluated at the end of our simulations.The expansion dynamics of the universe during the post-inflationary stage can be parameterized by the time-dependence of its equation-of-state w(a).Straightforward integration of the continuity equation ρ + 3Hρ(1 + w) = 0 allows to express the total energy-density ρ at a given time as (2.30) We now introduce a redshift factor accounting for time-dependent deviations with respect to a strict radiation-like (i.e.w = 1/3) dominated universe [43] ε(a) ≡ a a end defined in terms of the instantaneous (logarithmically) averaged equation-of-state w(a) ≡ 1 log(a/a end ) a a end w(ã) d log ã . (2.32) By construction, values of ε(a) are larger (smaller) than unity for an EOS w > 1/3 (w < 1/3).The redshift factor ε(a) also tends to a constant value after fragmentation, for a ≥ a RD ≫ a frag , as the EOS asymptotes 1/3.The numerical evaluation of this factor from our lattice simulations is represented in Fig. 3, showing how ε(a) exhibits an asymptotic value reached at late time as the oscillation-averaged EOS of Fig. 2 asymptotes 1/3.In what follows, abusing somewhat notation, we will refer to the asymptotic value of the redshift factor by simply dropping the argument, ε ≡ ε(a ≥ a RD ).In terms of it we can express the scale factor during the radiation-like dominated era as which holds in particular at reheating a = a reh .This allows us to express the scale factor at the end of inflation as where we have used entropy conservation, ρ rad,0 a 0 (g reh /g 0 ) 1/4 (g s,0 /g s,reh ) 1/3 .Here ρ rad,0 denotes the radiation energy density at the present day, and g ρ and g s are the energy density and entropy degrees of freedom, respectively, which have been evaluated at reheating and present day.
A simple analytical estimate of the parameter ε can be obtained by assuming instantaneous fragmentation and approximating the EOS parameter posterior to inflation by a step function with Θ being the Heaviside function.This immediately yields  ).The small deviation from ε = 1 for k = 4 comes from the fact that the inflation potential slightly departs, in the first post-inflation instants, from an exactly quartic potential.Similar effects are also relevant for larger k but less significant as ε ≫ 1.Our numerical treatment allow us to fully capture these effects.

Reheating
After inflation ends, the energy density of the inflaton would have been transferred to relativistic degrees of freedom via the reheating mechanism.Reheating can be achieved by coupling the inflaton to other species that could be fermions (f ) or bosons (b) and that could belong to (or could be connected to) the Standard-Model thermal plasma.For simplicity, in this work we consider the energy transfer from the inflaton to a relativistic species without specifying their interaction properties with respect to the Standard-Model.Such energy transfer could be induced via renormalizable interaction terms in the Lagrangian such as where y, σ (µ) are dimensionless (dimensionful) parameters.At the early stages of reheating, prior to fragmentation, the particle production process is dominated by the transition from the homogeneous oscillating condensate ϕ(t) to the final-state free particles.Disregarding collective enhancement/blocking effects (preheating), which can be done for sufficiently small Lagrangian couplings, the interactions (2.37) induce a source term for the energy density of the radiation component (R) of the form [20,55] ρR in addition to modifying the Friedmann equation (2.39) The corresponding energy-transfer rate from the inflaton condensate into bosonic or fermionic final states can be expressed as with the time-dependent inflaton mass given by m ϕ (t) . The effective couplings appearing in Eq. (2.40) can be expressed for fermion decay products as (2.41) Similarly, bosonic effective couplings can be expressed as [20] and (2.43) These effective couplings depend on the details of the oscillations throughout the exponent k.In addition, they involve the efficiency parameters α y (k, R), α µ (k, R) and α σ (k, R) which account for the kinematic effects resulting from effective masses m eff of the decay products induced by the time-dependent inflaton field vev, R ∝ (m eff /m ϕ ) 2 ϕ→ϕ 0 .Details can be found in Ref. [20].If the decay of the inflaton field is not completed prior to fragmentation, then the completion of the reheating process relies on the decay of the free inflaton particles (fluctuations), which now dominate the energy budget of the universe, into Standard-Model fields.In this case, the continuity equation (2.38) is modified to [54][55][56] ρR where n δϕ denotes the free inflaton number density, and Γ δϕ is the free inflaton particle decay rate.
To obtain an exact solution, the evaluation of the source term of (2.44) requires the numerical determination of the inflaton occupation numbers from the lattice code [54,56].Nevertheless, this source term can be approximated upon the estimation of the non-relativistic particle energy density in terms of the corresponding Lorentz boost factor, δϕ . (2.45) The full expressions for the particle production rates can be found in [55].
The reheating time a reh can be defined by the condition ρ R (a reh ) = ρ ϕ (a reh ) where ρ ϕ accounts for both contributions from the condensate and the fragmented quanta.The conditions for reheating to occur after fragmentation, i.e. a reh > a frag , were derived in Ref. [55] and are summarized here.
For fermionic final states, requiring that reheating is achieved after fragmentation is equivalent to the condition a frag a end where we introduced the quantity parametrizing the power of m ϕ appearing in the production rates.For 8 + k − 6kl < 0, this condition becomes In the case where inflaton decay products are scalars, the post-fragmentation reheating condition imply a frag a end for the dimensionful coupling and a frag a end for the dimensionless coupling.An important point to emphasize is that since the equation of state w = 1/3 is maintained throughout the transition between the quanta dominated universe and the thermal plasma era, the precise moment of reheating does not play any role on relevant physical observables.Particularly, as long as the previously stated conditions are satisfied, the evolution of the scale factor would remain independent on the inflaton couplings to the radiation bath and the gravitational wave signal computed in the following section would remain unaffected. 4ummary.The conditions (2.48)-(2.50)ensure that reheating is achieved after the fragmentation process is completed and that the equation of state transits smoothly from w = (k − 2)/(k + 2) (where the conformal Hubble rate scales as H −1 ∼ a 2(k−1)/(k+2) ) during the coherently-oscillating phase to w = 1/3 (with H −1 ∼ a) during subsequent phases of the universe dominated by fragmented inflaton quanta or ultra-relativistic particles belonging to the thermal plasma.This completes the transition from cosmic inflation (with H −1 ∼ a −1 ) to the standard phase of radiation domination.The cosmological history of our setup in summarized in Fig. 4, where the evolution of the comoving Hubble horizon H −1 is presented as a function of the scale factor throughout the various phases for each value of k = 4, 6, 8, 10 considered in this work.

Gravitational waves
In this section we detail the computation of the gravitational wave signal generated throughout the entire cosmological history, accounting for the transition from cosmic inflation to a universe dominated by a coherently oscillating condensate followed by the domination of inflaton quanta finally a thermal bath of relativistic particles.Limiting ourselves to tensor metric fluctuations, the perturbed Friedmann-Robertson-Walker metric can be parameterized by Gravitational waves are described by the metric perturbation h ij which represents two independent tensor degrees of freedom satisfying the transverseness and tracelessness (TT) conditions ∂ i h ij = 0 and h i i = 0.In cosmological perturbation theory, tensor metric perturbations are decoupled from scalar perturbations at first order.Scalar perturbations only appear as a source of tensor metric perturbations at second order through the transverse-traceless component of the anisotropic stress Π TT ij = ∂ i ϕ∂ j ϕ TT .Accounting for this source term, the corresponding equation of motion for tensor modes expressed in Fourier space reads where p is the (comoving) Fourier scale, p = |p| and primes denote differentiation with respect to conformal time.As Eq. (3.2) is linear in h ij , its solutions can be decomposed in terms of the sum of homogeneous and inhomogeneous solutions.During singe-field slow-roll inflation the source term is essentially negligible making the inhomogeneous solution irrelevant.For modes wavelengths smaller than the Hubble radius at all times (p > H), the homogeneous solution to this equation describes quantum fluctuations in the Bunch Davies vacuum and do not contribute to the Gravitational Wave (GW) signal.Modes experiencing a super-horizon phase (p < H) during inflation would be excited, freezing-out outside the horizon, and contributing to the homogeneous solution if they remain outside the horizon, or sourcing the B-mode polarization spectrum of the CMB if they re-enter the horizon during the recombination epoch.The inhomogeneous solution is only significant after the end of inflation, once non-linearities induce sizable mode-mode couplings resulting in a non-vanishing anisotropic stress, efficiently sourcing gravitational waves during the fragmentation process.This section is dedicated to estimating contributions coming from the inhomegenous solution (i.e.induced by the fragmentation process) and the homogenous solution (i.e.tensor modes excited during inflation and re-entering during a stiff EOS era) and reconstructing the total gravitational wave signal.

Gravitational waves from fragmentation
As mentioned previously both contributions to the GW signal can be estimated separately, in this subsection we estimate the GW signal generated solely during the fragmentation process.At the onset of fragmentation, non-linearities become large enough to contribute significantly to gravitational wave production by acting as a source term through the anisotropic stress in the right-hand side of the equation of motion (3.2).Projected onto the transverse-traceless component in Fourier space, the anisotropic stress can be written explicitly as [33] Π TT ij (p, τ ) = P iℓ (p)P jm (p) − with P ij = δ ij − pi pj , pi = p i /p and Π ℓm (p, τ ) = d 3 q (2π) 3/2 q ℓ q m ϕ(q, τ ) ϕ(p − q, τ ) . (3.4) We simulate numerically the production of gravitational waves with CosmoLattice which implements a discretized TT projector.5

Comoving scales and frequencies
The gravitational wave frequency f at the present epoch can be related to the comoving Fourier scale p via [33] g s,0 g s,reh where we introduced the dimensionless comoving momentum κ ≡ p/(a end √ λM P ).An evaluation of this expression gives the following typical frequencies where normalization for κ has been appropriately chosen as the typical scales prone to parametric resonance.A range of scales close to the typical κ is represented as colored bands in Fig. 4. We denote the largest Fourier scale to leave the horizon during inflation by p end ≡ H end , and the scale which re-enters the horizon at fragmentation by p frag ≡ H frag .Both scales represent respectively the scale tangent to the minimum of H −1 and the scale at which the vertical line corresponding to a frag and H −1 meet in Fig. 4, for each value of the exponent k.Under the approximation of instantaneous fragmentation, the corresponding frequencies can be related via where f end can be expressed in term of the energy density at the end of inflation through a end via Eq.(2.34).We find (3.8)

Energy density and spectra
The (normalized) GW energy density per logarithmic frequency interval can be expressed as where ρ GW is the total GW energy density and ρ c is the (time-dependent) critical energy density defined via ρ c ≡ 3H 2 M 2 P .From the lattice we can infer the full spectral and time dependence of the normalized GW energy-density in Eq. (3.9).The total GW energy density relative to the critical density at a given time can be computed in a straightforward way The right-hand-side is also estimated at a given time by CosmoLattice as an average over the lattice volume V (see Ref. [78] for details) and reads with hij ≡ ah ij .This quantity is represented as a function of the scale factor in Fig. 5 on the left panels for k = 4, 6, 8, 10.One can notice a sharp increase of the relative gravitational wave energy density for a ≃ a frag when non-linearities become significant.The spectral dependence of the GW energy density is represented in Fig. 5 on the right panels in terms of the present day frequency.The colors of the various curves on the right panels code the evaluation time which can be inferred from the left panels.From the right panels of Fig. 5 one can see scales prone to parametric resonances progressively redshift towards the infrared with time for k ̸ = 4,6 and being imprinted in the primordial tensor spectrum.While the resonance structure for k = 4 can be observed in the final spectrum, which can be estimated from the Boltzmann approximation [54], the time dependence of the resonance bands for k = 6, 8, 10, and the subsequent non-linearities in the scalar perturbations, would eventually erase any remnant of such resonances by generating a relatively broad spectrum.This effect was already noticed at the level of the inflaton-fluctuation occupation number in Ref [55].
The relative energy density after fragmentation, deep into the inflaton-quanta dominated era at a = a RD ,7 can be deduced from the asymptotic values reached on the left panels of Fig. 5 and correspond to which is the largest for k = 4 and decreases as k increases.This quantity can be related to the present GW energy density The present day normalized GW energy density can thus be expressed as where we took g reh = g s,reh = 106.75,g = 3.909 and g 0 = 3.363.The total energy density at the present day is The energy density corresponding to the gravitational wave signal can potentially increase the effective number of relativistic species ∆N eff ≡ N eff − N SM eff with respect to the Standard-Model expected value N SM eff = 3.046.This contribution can be estimated as yielding ∆N eff at most of order O(10 −5 ) for k = 4 and smaller by one or two orders of magnitude for larger values of k.Such small values of ∆N eff are far below current and upcoming sensitivity such as COrE [81] or Euclid [82] that are expected to reach ∆N eff < 0.013 at the 2σ level.

Tensor modes excitation from inflation
As discussed above, during slow-roll inflation the source term in the equation of motion for the tensor modes, i.e. the right-hand side of Eq. (3.2), can be neglected and the equation of motion reduces to its homogeneous form The rescaled metric perturbation hij ≡ ah ij can be decomposed as a sum over polarization states where ϵ λ ij (p) are the normalized spin-2 polarization tensors satisfying ij ϵ λ ij (ϵ λ ′ ij ) = 2δ λλ ′ .Each polarization component of the rescaled tensor mode satisfies then the following oscillator equation with the corresponding frequency given by ω 2 p = p 2 − a ′′ /a.By assuming Bunch-Davies initial conditions for the mode functions deep inside the horizon (|τ one can solve the mode-function equation and show that tensor perturbations become frozen after horizon crossing p ≃ H ≡ aH yielding h λ ∼ p −3/2 at the end of inflation and a corresponding power spectrum evaluated at horizon crossing.The result, as it is well known, is an approximately scale-invariant power spectrum which remains frozen for superhorizon scales after the end of inflation regardless of the EOS of the universe. Short wavelengths.Once tensor perturbations re-enter the horizon, the mode function starts oscillating with a decaying envelope h λ ∼ p −3/2 (a p /a) where a p is the scale factor at horizon reentry p ≡ H(a p ).This induces a spectral distortion of the frozen spectrum [60,65,83] through the time dependence of the conformal Hubble rate H ∼ a −2(k−1)/(k+2) implying a p ∼ p −(k+2)(2k−2) and the corresponding energy-density spectrum [84] Ω GW h 2 ≃ 1 24 Ω rad h 2 g reh g 0 g s,0 g s,reh where T p (a) = (a p /a) 2 ∼ p −(k+2)(k−1) is the transfer function accounting for the time evolution of the tensor spectrum from its frozen value.This implies an energy-density spectrum scaling as 1) .For scales re-entering during radiation domination (p < p frag ) or during inflaton-quanta domination with k = 4, this yields a flat energy-density spectrum where the potential is evaluated at the field value ϕ * reached at the time the CMB pivot scale crosses the horizon during inflation which can be inferred from Eq. (2.5) for a given k. 8 For scales p > p frag entering during an inflaton-quanta dominated universe, the spectrum becomes blue tilted from the momentum dependence of the transfer function up to the largest scales that experienced a super-horizon phase during inflation k end which are represented in Fig. 4. Switching to present-day frequency via Eq.(3.5), the total spectrum can therefore be expressed as with f frag and f end given by Eq. (3.8) for each value of k.
Long wavelengths.On much larger physical scales, such as the CMB fiducial scale k * = 0.05 Mpc −1 , primordial tensor perturbations are customarily parametrized in terms of the tensor-to-scalar ratio which can be expressed as [20] Such tensor perturbations could leave traces in B−polarization modes of the CMB and could potentially be detected by a future experiment.The predicted values of r from Eq. (3.25) appears to be close to the sensitivity reach for the upcoming Simons Observatory (SO) r < 6 × 10 −3 [85].Moreover, future missions such as LiteBIRD and CMB-S4 should allow to disprove or confirm this model by reaching sensitivities as large as r < 2 × 10 −3 [86] and r < 10 −3 [87] respectively.

Additional contributions to gravitational wave production
There have been substantial developments in the computation of gravitational wave production from the reheating dynamics where various contributions have recently been estimated: • GW from graviton bremsstrahlung during reheating, from inflaton decay or inflaton annihilation [46][47][48][49][50].These contributions can be important for large reheating temperatures but suppressed at low reheating temperatures.We expect them to be subdominant in our case.
These contributions are thus expected to be subdominant compared to the signal estimated in this work.

Results and discussion
The gravitational wave spectra resulting from tensor modes excited during inflation, i.e.Eq. (3.24), are represented for each value of k = 4, 6, 8, 10 in Fig. 6 as dotted lines.The gravitational wave signals sourced by non-linearities from the fragmentation process, i.e. the red lines in Fig. 5 right panels, are in turn shown as solid lines.We also represented sensitivity prospects for future gravitational wave observatories, recently reviewed in [88]: the Big Bang Observer (BBO) [89,90], the (ultimate) DECihertz Interferometer Gravitational wave Observatory (uDECIGO) [91][92][93], the Laser Interferometer Space Antenna (LISA) [94], the Einstein Telescope (ET) [95][96][97] in addition to estimate from resonant cavities (Res.cav.) of Ref. [98].This figure illustrates the complementary between O(Hz) detectors such as BBO or DECIGO and O(GHz) cavities in probing gravitational wave production in our setup.
Hz detectors.BBO should almost reach the sensitivity required to detect the scale-independent part of the spectrum at low frequency for k = 4, 6, 8 but should be able to detect the Ω GW,0 h 2 ≃ 10 −15 (f /Hz) 2/3 signal for k = 10.The sensitivity of uDECIGO should be large enough to probe the scale invariant part of the spectrum for all values of k considered in addition to the Ω GW,0 h 2 ≃ 8 Given our normalization choice for λ, the value V (ϕ * ) is the same for all k.
GHz detectors.Electromagnetic cavities have recently garnered interest as a possibility to detect high-frequency gravitational waves.Future experimental setups could pave the way for a new frontier in GW detection [99].Potential high frequency GW detectors with a sensitivity as high as proposed in Ref. [98] should be able to resolve the peaked structure of the spectrum generated for k = 4, pointing undoubtedly towards an inflaton potential quartic about the minimum.Moreover, such experiment could detect tensor modes excited during inflation and re-entering the Hubble radius during the coherently-oscillating inflaton era for both k = 10 at Ω GW,0 h 2 ∼ 10 −10 and k = 8 at Ω GW,0 h 2 ∼ 10 −13 − 10 −12 .Since both signals are also in the reach of uDECIGO at much lower frequencies, the combination of results from both experimental facilities could be used to identify or disprove the underlying inflation scenario.Nevertheless, fragmentation signals for k = 6 and k = 8 appears not to be in the reach of future experimental facilities.

Summary and conclusion
In this work we have investigated a consistent picture of the universe describing its early states, starting from inflation and arriving to a phase dominated by a bath of ultra-relativistic particles.We considered for definiteness the T-model of inflation, characterized by an exponent k that we assigned to k = 4, 6, 8, 10 as specific realizations.While cosmic inflation, with an equation of state w ≃ −1, is sensitive to the asymptotically flat part of the inflaton potential at large field values, the exponent k controls the expansion history close to the minimum, where V (ϕ) ∼ ϕ k .This post-inflation history is additionally determined by the inflaton coupling to states that eventually make up the primordial radiation bath.In this work we explored the possibility that this coupling is sufficiently small that reheating would only occur after a significant amount of time posterior to inflation, i.e. late reheating.Immediately after inflation, coherent oscillations of the inflaton field about the minimum of its potential generate a phase of expansion with an oscillation-averaged equation-of-state w ≃ (k − 2)/(k + 2).Subsequently, excited by the oscillations, parametric resonant effects induce the growth of inhomogeneities in the form of copiously produced inflaton quanta eventually reaching the non-linear regime.Significant mode-mode couplings trigger the fragmentation of the inflaton condensate, marking the transition between a coherently-oscillation phase to an era dominated by a collection of relativistic inflaton quanta, corresponding to an equation of state w ≃ 1/3.Under the assumption of weak couplings, the transition to a universe dominated by a Standard Model thermal bath in equilibrium (i.e.reheating) is therefore achieved only once fragmentation is completed.This transition leaves the equation of state unaffected ≃ 1/3, therefore the precise value of the reheating temperature does not affect the cosmological history, and the subsequent phenomenological consequences explored in this work.
In this paper, we kept track of the expansion history by solving the equation of motion for the inflaton field during inflation and by simulating the post-inflation history with CosmoLattice.We estimated the time dependence of the equation-of-state during the fragmentation process of the inflaton.While fragmentation is relatively fast for k = 4, occurring at a frag /a end ≃ 180, it takes longer to be achieved as the value of k increases, yielding a frag /a end ≃ 4.5 × 10 4 , 6 × 10 6 , 7 × 10 8 for k = 6, 8, 10 and a longer period of coherent-oscillations for large values of k.
We estimated the gravitational wave signal generated during the cosmological history in our setup and identified two main contributions.The first contribution consists of tensor modes crossing the horizon and being excited during inflation.Such modes could re-enter during a phase of the universe with an equation of state w = 1/3 such as inflaton-quanta or radiation dominated phase yielding a scale invariant signal Ω GW,0 h 2 ≃ 5.6×10 −18 for k = 4, 6, 8, 10.Modes re-entering during the coherent inflaton-oscillation phase would benefit from a stiff equation of state w = (k − 2)/(k + 2) inducing a blue-tilted signal Ω GW,0 h 2 ∼ f (k−4)/(k−1) .While the transition between the scale-invariant and the blue-tilted spectra is determined by the size of the Hubble radius at fragmentation, the ultraviolet cutoff corresponds to the Hubble radius at the end of inflation.The second contribution to the GW signal comes from the large inhomogeneities generated during the fragmentation process sourcing tensor perturbations at second order in cosmological perturbation theory.This results in a high frequency signal for k = 4 with a peaked structure at around f ∼ 10 8 Hz as remnants of the parametric resonance structure.For larger values of k = 6, 8, 10, non-linearities induce broad spectra at frequencies respectively of order f ∼ 10 7 , 10 5 , 10 2 Hz.
We compared our results to prospects for experimental facilities sensitive to gravitational waves.We identified three possibilities of probing the gravitational wave signals in the future.First, the scale-invariant part of the spectrum on very large physical scales could be detected in the form of primordial tensor modes by LiteBIRD and CMB-S4 experiments for all k = 4, 6, 8, 10 considered.Second, uDECIGO should reach a sensitivity sufficient to detect the scale-invariant spectrum for frequencies f ∼ Hz for k = 4, 6 and the blue-tilded component of the spectrum for k = 8, 10 as well as the fragmentation-induced signal for k = 10.At last, resonant electromagnetic cavities sensitive to frequencies f ∼ 10 8 Hz should be able to detect the blue-tilted spectra for k = 8, 10 and the fragmentation signal for k = 4.
Our work shows the ability for future experimental facilities to detect gravitational wave signals emitted in the early universe as well as the complementarity between detectors sensitive to different scales corresponding to different epoch of the universe.Gravitational wave astronomy will be an essential probe of inflation scenario and determinant to discriminate various interpretations in the decades to come.While strong efforts are currently employed to develop high-frequency GW detectors [99], we advocate the necessity of building such ultra-sensitive GW observatories to shed light on the dynamics of inflation and reheating.

Figure 1 :
Figure 1: Time evolution of the normalized periodic component of the inflaton field over an oscillation period.

Figure 2 :
Figure 2: Energy density of the inflaton in the form of condensate (ρ φ, blue), in free quanta (ρ δϕ , yellow) and the sum of both contributions (ρ ϕ , green) in the top panels.Oscillation-averaged equation of state in the bottom panels as a functions of the scale factor for k = 4, 6, 8, 10.

Fig. 3 .
In this same Figure we show for comparison the analytical estimates of ε(a) based on Eq. (2.35

Figure 3 :
Figure 3: Analytical (black dashed lines) and numerical (red lines) evaluation of the redshift parameter ε defined in Eq. (2.31) as a function of the scale factor for k = 4, 6, 8, 10.

Figure 4 :
Figure 4: Comoving Hubble radius and relevant scales normalized to the value at the end of inflation H −1 end at a = a end for k = 4.The corresponding gravitational wave frequency at the present day is represented on the vertical axis on the right.Colored bands represent the scales prone to parametric resonances that are typically excited during the fragmentation process (see (3.6)).

Figure 5 :
Figure 5: Normalized gravitational wave energy density generated solely during the fragmentation process: as a function of the scale factor (left panel) and as a function of present-day frequency (right panel).The color on the right plots corresponds to an evaluation time that matches the color on the left plots.

Figure 6 :
Figure 6: Total stochastic gravitational wave background as a function of the present-day frequency.Purple, green, blue and red curves correspond to k = 10, 8, 6, 4 respectively.The solid lines represent the GW signal from non-linearities during the fragmentation process and dotted lines GW resulting from tensor mode super-horizon excitation during inflation.Orange-dashed line are sensitivity prospects for future experimental facilities, see Sec. 3.4 for details.