Reheating after inflaton fragmentation

In the presence of self-interactions, the post-inflationary evolution of the inflaton field is driven into the non-linear regime by the resonant growth of its fluctuations. The once spatially homogeneous coherent inflaton is converted into a collection of inflaton particles with non-vanishing momentum. Fragmentation significantly alters the energy transfer rate to the inflaton's offspring during the reheating epoch. In this work we introduce a formalism to quantify the effect of fragmentation on particle production rates, and determine the evolution of the inflaton and radiation energy densities, including the corresponding reheating temperatures. For an inflaton potential with a quartic minimum, we find that the efficiency of reheating is drastically diminished after backreaction, yet it can lead to temperatures above the big bang nucleosynthesis limit for sufficiently large couplings. In addition, we use a lattice simulation to estimate the spectrum of induced gravitational waves, sourced by the scalar inhomogeneities, and discuss detectability prospects. We find that a Boltzmann approach allows to accurately predict some of the main features of this spectrum.


Introduction
Ever since the inception of cosmic inflation as a potential solution to the initial condition problems of standard Big Bang cosmology, the question of how the universe transitions from a cold, empty, quasi-de Sitter state to a radiation dominated stage in thermal equilibrium, has been an active research topic.The original inflationary proposal (old inflation) was in particular abandoned due to its impossibility to lead to a phenomenologically successful reheating of the universe [1][2][3][4].The overwhelming majority of the subsequent proposals, based on the slow-roll of an elementary scalar field called the inflaton, incorporate reheating mechanisms tied to the coherent oscillation of this inflaton field about its minimum. 1t the perturbative (Boltzmann) level, reheating is modeled as the dissipation of the oscillating inflaton into elementary particles, with a decay rate determined by the quantum mechanical transition amplitude from the time-dependent vacuum to the corresponding particle states.This amplitude is evaluated by averaging the transition rate over the fast oscillation of the inflaton about its minimum, assuming therefore that particle production evolves adiabatically following the expansion of the universe.For two-body decays, the dissipation rate is mostly sensitive to the time-dependence of the effective mass of the oscillating inflaton, which can lead to time-dependent kinematic blocking effects, and in general to a decreasing (or increasing) decay efficiency depending on the shape of the minimum of the inflaton potential [7][8][9][10][11][12].
It is known, however, that the perturbative picture is insufficient to describe the dynamics of reheating in the presence of strong couplings, or when the short time-scale oscillation of the inflaton field leads to the resonant enhancement (bosons) or suppression (fermions) of the growth of the quantum fields associated with the decay products.When active, these collective effects, known as preheating, generically lead to a significant modification of the reheating dynamics [7,[13][14][15][16][17][18].In scenarios with a quadratic inflaton potential, the parametric resonance of the momentum modes of the produced fields can lead to the exponential growth of perturbations in the early stages of reheating.This growth is transient and disordered, due to a mismatch between the particle production rates and the expansion rate of the universe.Nevertheless, if the resonance is sufficiently strong and maintained for a sufficient amount of time, it can bring the growth of fluctuations into the non-linear regime.The inflaton condensate is fragmented in favor of a quasi-thermal bath of inflaton particles and its decay products [19][20][21][22][23].In some scenarios, the result is that the Boltzmann approximation correctly describes the production of daughter fields only in a very narrow range of the available parameter space [24].
Preheating in non-quadratic minima has also been studied extensively [21,[25][26][27][28][29][30][31][32][33][34][35][36][37].Notably, the presence of inflaton self-interactions, even if small, can accumulate over several oscillations and source the resonant growth of non-zero momenta inflaton modes, eventually backreacting with the homogeneous component.This self-fragmentation of the inflaton not only re-distributes the initial energy density of the field, but it can also lead to the formation of localized, soliton-like objects such as oscillons [28,31,32,[38][39][40].The large inhomogeneities triggered by the self-fragmentation of the inflaton can also efficiently source gravitational wave (GW) production.The typical frequency for the resulting stochastic gravitational wave background is larger than its counterpart sourced during inflation f ≳ MHz [41][42][43].Such high frequency is beyond the reach of current future earth or space based interferometers such as LISA or the DECIGO but might be in the reach of future dedicated experiments [44].The resulting gravitational wave spectrum would carry precious information about the dynamics of the inflaton in the first instants following inflation and shed light on the reheating epoch. 2he depletion of the coherent condensate eventually shuts down the parametric resonance.Following fragmentation, one must therefore follow the population of the relativistic bath of decay products by perturbative means, accounting for both the inflaton zero and non-zero mode dissipation.The decay rates will differ from those prior to the backreaction epoch, significantly modifying any estimates for the reheating temperature of the universe.A formalism to estimate the evolution of the radiation energy density during reheating in the pre-and post-fragmentation regimes, including the reheating temperature is the main goal of this work.To showcase the differences with respect to the purely perturbative evolution of the inflaton-radiation sector we specialize to the case of a quartic inflaton potential near its minimum.
This paper is organized as follows.In Section 2 we discuss in detail the evolution of an inflaton field in a quartic potential during reheating, in the absence of interactions with other fields.Sec.2.1 is devoted to the dynamics of the background evolution of the homogeneous inflaton condensate.In Sec.2.2 we study the parametric growth of inflaton fluctuations at linear order, and in Sec.2.3 we study the non-linear regime of their growth, including backreaction and fragmentation effects.For Sec. 2.4 we include the perturbative Boltzmann determination of the spectrum of inflaton fluctuations.The inflaton-matter/radiation couplings are introduced in Section 3. The analysis is divided into the study of the dissipation of the coherent inflaton condensate (Sec.3.1.1)and the decay of the inflaton particles after fragmentation (Sec.3.1.2).Our main results, the corresponding reheating temperatures, are presented in Sec.3.2.Finally, in Sec. 4 we determine the induced gravitational waves from the non-linear scalar dynamics, and discuss their potential for detectability.Our conclusions are presented in Sec. 5.

Post-inflationary dynamics in a quartic potential
The current measurements of the primordial curvature power spectrum, and bounds on the tensor power spectrum, are compatible with the presence of a single, slowly rolling neutral scalar field which drives inflation.We denote this inflaton by ϕ, and for it we will assume the following form for the action, Here g ≡ det(g µν ) is the metric determinant of a flat Friedmann-Robertson-Walker metric with scale factor a, V (ϕ) denotes the inflaton potential, and L int denotes its interactions with the rest of the (extended) Standard Model.The potential V (ϕ) must be chosen in order to match the measured amplitude and tilt of the scalar power spectrum, A S * ≃ 2.1 × 10 −9 , n s ≃ 0.966, and to avoid the upper bound on the tensor-to-scalar ratio, r < 0.036 [46][47][48].Moreover, in this work we impose the condition that this potential is quartic during the post-inflationary reheating.Among several suitable candidates, we choose for definiteness the quartic T-model attractor [12,49], where M P = 1/ √ 8πG N ≃ 2.435 × 10 18 GeV denotes the reduced Planck mass.The constant λ is determined by the scalar spectrum amplitude, evaluated at the horizon exit time of the Planck pivot scale, k * = 0.05 Mpc −1 .For a generic potential, the number of e-folds before the end of inflation at horizon exit time are determined from the expression [50,51] Here H 0 = 67.36km s −1 Mpc −1 [46], T 0 = 2.7255 K [52] and a 0 = 1 denote the present Hubble parameter, photon temperature and scale factor, respectively.The energy density at the end of inflation is denoted by ρ end , and the energy density at the beginning of the radiation dominated era by ρ rad .The effective number of degrees of freedom during reheating is denoted by g reh .The e-fold averaged equation of state parameter during reheating corresponds to As is well know, and we show explicitly in the following section, for quartic reheating w ≃ 1/3.The last term in (2.4) vanishes, and therefore the number of e-folds is uniquely determined, N * ≃ 56 for T-model inflation, and so is the coupling constant λ ≃ 3.3 × 10 −12 .We now study the postinflationary evolution of the inflaton sector, from the coherent oscillation stage, to the backreaction and post-fragmentation phase.

Coherent oscillations
We first consider the dynamics of the inflaton field at the background level after the end of inflation.Inflation ends when ä = 0, or equivalently when φ2 = V (ϕ), which for the T-model (2.2) corresponds to ϕ end ≃ 1.52 M P , ρ end ≃ (4.5 × 10 15 GeV) 4 . (2.6) Our focus here will be the early stage of reheating, with a/a end ≲ O( 102 ).Importantly, we will assume that the main decay channel of the inflaton is to fermionic fields, which do not manifest the exponential growth from the parametric resonance of their mode equations.This allows us to assume that the coherence of the inflaton oscillation is maintained, and that the third term in (2.1) can be disregarded at early times [11,12].If this is the case, variation of the action with respect to the homogeneous inflaton field and metric yields the equations of motion where H = ȧ/a is the Hubble parameter, and a dot represents differentiation with respect to cosmic time.These equations describe an underdamped anharmonic oscillator, which can be parametrized in terms of an envelope function ϕ 0 (t), encoding the redshift due to expansion, and a quasi-periodic function P(t), which encodes the short time-scale oscillation, ϕ(t) ≃ ϕ 0 (t) P(t) . (2.9) To determine the time-dependence of the envelope function, we multiply Eq. (2.7) by ϕ and average over one oscillation.This yields, implying that the oscillation-averaged energy and pressure densities can be written as [12] where we have used ⟨P 4 ⟩ = 1 3 , as it can be checked from the solution for P, see (2.15) below.Hence, the equation of state parameter is w ϕ ≃ 1/3, and the equation of motion (2.7) can be rewritten as ρϕ + 4Hρ ϕ ≃ 0 .
(2.13) Thus, as is well known, the energy density of the inflaton redshifts as radiation, ρ ϕ ∝ a −4 .As a result, we determine that the decaying envelope of the oscillation redshifts as ϕ 0 (t) ≃ ϕ end a end a(t) . (2.14) The time-dependence of the quasi-periodic function P(t) can be obtained by exploiting the underdamped nature of the oscillations.Approximating the envelope as constant during one oscillation, the equation of motion (2.7) can be written as where the effective mass of the inflaton is defined as where the sum is taken over the harmonic modes of the ϕ oscillation with the fundamental angular frequency given by In order to account for expansion, we can indirectly exploit the conformal invariance of the model.Since a ∝ t 1/2 , conformal time is τ ≡ a −1 dt ∝ a, and m ϕ t = m end (τ − τ end ).We can therefore characterize the stage of coherent oscillations of the inflaton by means of the approximate solution

Parametric resonance
In the Introduction we have made the case for the study of the growth of the inflaton fluctuations during reheating, driven by the self-interaction.We now proceed to study the resonant growth of the inflaton momentum modes in the linear approximation.For these purposes, we denote the nonhomogeneous inflaton perturbation as δϕ(t, x).Variation of the action (2.1), disregarding couplings of ϕ to other degrees of freedom, leads to the following dynamical equation, Introducing the canonically normalized fluctuation X ≡ a δϕ, and switching to conformal time, we can write the canonically quantized field as where k denotes the comoving momentum, and âk and â † k are the annihilation and creation operators, respectively, satisfying the canonical commutation relations We ensure that the corresponding canonical commutation relations between the field, X k , and its momentum conjugate, X ′ k , are fulfilled by imposing the Wronskian constraint In this decomposition, the equation of motion (2.20) reduces to where ′ denotes the derivative with respect to conformal time.Examination of the factor in parenthesis shows that the second term decreases over time; for it we have where g 1 (z) and g 2 (z) are periodic functions, and µ k is a complex number, called the Floquet exponent.Exponentially growing solutions are therefore found if Re µ k > 0. In order to determine the Floquet exponents, we follow the eigenvalue method described in detail in [22].The result of this numerical computation is shown in Fig. 1.The simplicity of the Floquet chart is noteworthy, as only two main resonance bands appear.One at small momentum values, k/m end ≲ 1.8 × 10 −4 , and a second one, with a narrow width and Floquet exponent more than two orders of magnitude larger, at 0.71 ≲ k/m end ≲ 0.76.Due to the conformal nature of the quartic potential, the Floquet exponents are insensitive to the expansion of the universe, and a mode within a resonance band will stay in the band while the oscillation of the inflaton lasts.Therefore, modes with k/m end ≃ 0.7 will rapidly grow at the onset of reheating, and become non-linear after a few e-folds, even for a small coupling constant λ.

Backreaction and fragmentation
At linear order, the growth of the momentum modes of X can be tracked by solving Eq. (2.22) with the appropriate initial conditions, which we take as the positive-frequency Bunch-Davies vacuum where (2.27) Our main quantity of interest will be the phase space distribution (PSD) of the fluctuations, which coincides with the comoving occupation number of ϕ.Its UV-finite form is given by [18] Two other important quantities can be obtained from this PSD, the number density and the energy density of the fluctuations.Their UV-regular forms are computed as [18,23] ) (2.30) The resulting form of the PSD in the linear approximation is shown as the dashed lines in the left panel of Fig. 2. As expected, the PSD is peaked around k/m end = 0.73, in full agreement with the results of the Floquet analysis of the previous section.The height of this peak increases with time, as the resonant growth accumulates during the oscillation of ϕ.We only show three snapshots of the PSD for a ≤ 130, since the system rapidly evolves toward the non-linear regime, making the analysis based on the solution of (2.22) unsuitable.
The impact of the rapid growth of the resonant modes on the energy density of the inflaton fluctuations can also be appreciated in Fig. 3.We follow the scale factor dependence of ρ δϕ (the orange curve) by means of (2.30)only up to a/a end = 10 2 .Within this range, the growth of the energy density in fluctuations is clear, ramping up at a/a end ∼ 30, and growing exponentially fast afterwards.The transition from the linear to the non-linear regime is clear.
From this discussion it follows that, despite the smallness of λ, the growth of the momentum modes in the main resonance band is strong enough to eventually be able to draw an O(1) fraction of the energy density of the zero-mode of ϕ.The inflaton "fluctuation" is no longer such, and the system enters the non-linear regime.Mode-mode couplings become important, redistributing the energy of the resonant mode into other modes (rescattering), leading to large configuration space gradients, and effectively fragmenting the once homogeneous condensate.This phenomenon is known as backreaction, and it does not allow for a straightforward spectral analysis.Instead, the full non-linear configuration-space equation of motion for the full operator ϕ(t, x) must be solved, which is also a non-trivial task.To get around some of the complications, it is argued that, since in the non-linear regime occupation numbers are large, n k ≫ 1, the dynamics of the system can be adequately tracked by approximating ϕ as a classical field.The full non-linear PDE is then solved over a configuration-space lattice.The energy density of the inflaton is computed from the spatial average of the energy-momentum tensor of ϕ, which we denote with an over-bar, Spectral data in turn is obtained upon Fourier transformation of configuration-space quantities.The quartic inflaton system has been extensively studied in the past by means of lattice methods, see e.g.[21,[30][31][32][33][34][35]37].For our analysis we use the publicly available code CosmoLattice [60,61].
Our lattice results are summarized in Figs. 2 and 3.The left panel of Fig. 2 shows, as continuous curves, the form of the PSD of the inflaton as computed by the lattice code.At early times it shares the features of the linear approximation, showing a main growing peak, located at the momentum band predicted by the Floquet analysis.For the earliest times, the lattice PSD coincides with the linear prediction.Nevertheless, for a/a end ≳ 10 2 , the effect of rescattering becomes evident.The growth of the resonant mode is traded for the population of the PSD for lower momentum modes, and larger momentum adjacent modes.As time passes, the redistribution of energy becomes more efficient, and a UV tail in the distribution appears.At first this tail is populated more efficiently at some discrete values of k.Notably, these wavenumbers can be estimated by means of the Boltzmann approximation, as we show in Section 2.4.In any case, these features are quickly erased, and a smooth limiting distribution emerges, shown in red.The right panel of Fig. 2 conveys this asymptotic behavior of the PSD, in terms of the comoving number density of ϕ.At first n δϕ grows rapidly, driven by the parametric resonance.However, for a/a end ≳ 300, this growth stops, the comoving number density freezes, and the effect of rescatterings is merely the kinetic redistribution of energy among modes.
Fig. 3 contains the lattice results for the total energy density of ϕ, the density of the condensate, and that of free particles, for a/a end > 10 2 .For definiteness, and in order to connect with the results of the spectral analysis, we define the condensate component of the energy density as follows [23], that is, the energy density of the spatially averaged inflaton field.In turn, we take ρ δϕ = ρ ϕ − ρ ϕ .This definition of the fluctuation energy density matches excellently with the spectral result, and is shown as the orange curve in Fig. 3. Its value rapidly grows until a/a end ≃ 180, point at which the backreaction in the oscillating condensate is noticeable.The total energy density of ϕ redshifts as radiation even after rescattering becomes important, but the energy density of the zero mode is now efficiently transferred to fluctuations.Interestingly, despite the initial rapid growth in ρ δϕ , the dissipation of ρ ϕ after entering backreaction is gradual, ρ ϕ ∝ a −5.3 , unlike the exponential decrease of preheating scenarios in quadratic potentials [23].The survival of this coherent component will be an important ingredient in our exploration of the decay of ϕ of Section 3.

The Boltzmann approximation
In the previous sections we have shown that a combination of spectral and lattice methods is necessary to correctly determine the distribution of the inflaton fluctuations.We now take the opportunity to show that, in the linear regime, it is possible to extract non-trivial information for the PSD by means of the integration of the Boltzmann equation for the inflaton fluctuations.Inflaton quanta δϕ sourced from the inflaton background ϕ are induced by the interaction Lagrangian (2.34) Following [9,12,24], the Boltzmann equation in the presence of anharmonic oscillations of the inflaton background takes the form Here we denote the physical four-momenta of the fluctuations by P and P ′ .K n = (E n , 0) is the physical four-momentum of the inflaton condensate, where E n = n ω ϕ denotes the energy of the n th oscillating mode.M n represents the transition amplitude corresponding to the production of a pair of inflaton quanta |f ⟩ = |δϕδϕ⟩ from the n th Fourier mode of the coherently oscillating inflaton background field.
where Vol 4 denote the space-time volume.From Eq. (2.34), the matrix elements M n are found to be (2.37) The mean squared amplitude |M n | 2 accounts for a factor of 2 for identical particles in the final state.The Pn are the Fourier coefficients of the harmonic decomposition over one oscillation of the square of the quasi-periodic function P(t) defined in (2.9), The spatially homogeneous background inflaton condensate distribution can be expressed as the zero-mode Neglecting the backreaction of the inflaton quanta onto the condensate, the Boltzmann equation takes the following form where denotes the (time-independent) kinematic suppression factor.The quantum mechanical Bose enhancement factor, which appears in the right side of Eq. (2.40), can be scaled out by defining the "classical" distribution f c χ [62][63][64]: The Boltzmann equation reduces then to As the right-hand side of Eq. (2.43) is f c δϕ -independent, we can immediately integrate this equation to get [65,66] with q ≡ a(t)|P |/(a end m end ) the re-scaled comoving momentum referred to the end of inflation, denoted as k/m end in Fig. 2.
Peak structure.The distribution function corresponds to a collection of Dirac delta functions, located at different momenta (1/2)ncβ n .The β n coefficient is a decreasing function of n whose first (real) non-vanishing value is achieved for n = 5.By denoting the location of the i th peak as qi , the distribution can be equivalently expressed as as sum of contributing modes via with time-dependent coefficients F c i (t) and peak locations determined by where i = n − 4.This ratio is remarkably independent of any model parameter.It is straightforward to check that the first peak q1 is predicted to be at q1 ≃ 0.7, in excellent agreement with the Floquet analysis and the linear and lattice results.The relative location of peaks at larger momenta can be straightforwardly inferred from Eq. (2.46) and corresponds to the following ratio qi qj = (i + 4)β i+4 (j + 4)β j+4 . (2.47) Limitation of the Boltzmann approach.The previous computation resulted in a distribution peaked at discrete values of momenta, with an infinitesimally narrow width and a large amplitude peak characteristic of a Dirac delta function.In reality however, as we have seen in the previous sections, the distribution is spread around the peaks.Immediately after the end of inflation, from our linear approach we found a distribution function at the first peak f δϕ (q 1 , t end ) ∼ 0.2.Moreover, the mode amplitudes are expected to grow continuously F c i (t) ∝ a(t)/a end during reheating, becoming rapidly even larger.For such large amplitudes the Boltzmann approach breaks down almost immediately [24,63].
The amplitude of the peaks, therefore, cannot be accurately estimated from the Boltzmann approach starting from the first instants after the end of inflation.However, we find that this approach gives a qualitatively good description of the location of the PSD peaks appearing at the onset of backreaction, and notably, for the peaks in the induced gravitational wave spectrum.We postpone the detailed comparison to Section 4.

Reheating in a quartic potential
We now turn to the study of the decay of the inflaton into light degrees of freedom, necessary to complete reheating and populate the universe with a relativistic plasma in thermal equilibrium.As mentioned in the Introduction, our analysis will be based on the assumption that fields coupled to the inflaton are not strongly sourced via parametric resonance.To satisfy this assumption, we will implicitly assume that the main decay products of ϕ are spin 1/2 fermions, denoted as ψ.Moreover, we will assume for simplicity that the decay of the inflaton occurs at tree-level, and produces two particles in the final state.Thus, generically, we consider It is worth noting that the Yukawa term in the previous expression can lead to fermion preheating [67][68][69][70][71], for which the production of particles is resonantly suppressed, a consequence of Fermi-Dirac statistics.Perturbatively, the oscillating ϕ induces an effective mass for ψ than can kinematically block the decay into ψ very efficiently [12].The exploration of these interesting effects is not the main purpose of the present work, and will be postponed for a future study.We will therefore disregard these mechanisms, a possibility that arises for g ′ ≫ g.

Production rates
Schematically, the Boltzmann equation that must be integrated to track the phase space distribution of the decay products can be written as with C denoting the collision term.For generality (and future applications) let us for now assume a non-vanishing, and potentially time-dependent mass m ψ (t) for the decay products.Applying the operator d 3 P P 0 /(2π) 3 on both sides of (3.2), we obtain the following general form for the continuity equation for the energy density ρ ψ , ρψ + 4Hρ ψ − Here we introduce R(t), the radiation-energy production rate per unit of volume.
In the following we assume that the states ψ eventually thermalize among themselves and the rest of the Standard Model states.We therefore identify ρ ψ = ρ R , where the later denotes the energy density of the relativistic radiation plasma.The evaluation of the production rate in Eq. (3.3) will be split in two.We estimate the contributions to the production rate from the oscillating inflaton condensate ϕ and from the fragmentated inflaton quanta δϕ R(t) ≡ R ϕ (t) + R δϕ (t) . (3.4)

Decay of the coherent oscillations
Let us first briefly discuss the decay of the inflaton at early times, when it is composed mainly of the coherent, oscillating zero-mode.Under the assumption that quantum statistics do not play a role in the dissipation process, the Boltzmann equation which describes the decay of the inflaton can be written as [8,10,12] ρϕ + 3H(1 + w ϕ )ρ ϕ = −R ϕ (t) . (3.5) The rate in this equation is evaluated as where the decay rate Γ ϕ for a generic process ϕ → A + B is Here the sum is taken over the harmonic modes of the function P(t) over one oscillation, with associated energy E n = nω ϕ (see Eq. (2.17)).M n denotes the transition amplitude in one oscillation for each mode from the coherent state to the two-particle state that can be defined analogously to Eq. (2.36), with |f ⟩ = |A, B⟩ and L I = L int , cf. (2.1).In the case of an inflaton oscillating about a quartic potential, assuming a fermionic decay for which the masses can be disregarded, the decay rate Γ ϕ can be evaluated explicitly, taking the form [12] Γ ϕ→ ψψ (t) = α 2 y 2 8π m ϕ (t) , where y = {g, g ′ }, m ϕ is the envelope inflaton mass (2.16), and α ≃ 0.71 is an efficiency factor that codifies the anharmonicity of the ϕ oscillation [12].The continuity-Friedmann system that completes Eq. (3.5) corresponds to ρR + 4Hρ R = (1 + w ϕ ) Γ ϕ ρ ϕ , (3.9) with w ϕ = 1/3.

Decay of the fragmented inflaton
In Section 2.3 we discuss how the inflaton is fragmented after it enters the backreaction regime.When this occurs, the production of particles can no longer be characterized by (3.9) with decay rate (3.8).Instead, one must now determine the form of the continuity equation for ρ R from the microscopic Boltzmann equation, assuming a population of inflaton particles with PSD f δϕ (K).Disregarding inverse decays and Pauli blocking / Bose enhancement, the collision term in the Boltzmann equation (3.2) has the form (3.11) The particle production rate can then be evaluated in a straightforward way, and yields where denotes the decay rate of free inflaton particles δϕ.In the case corresponding to (3.1), Γ δϕ→ ψψ = y 2 m ϕ /8π for negligible masses for the outgoing states.

Total rate
Summarizing our previous findings, the total rate (3.4) accounting for production from the inflaton condensate and particle rates reads The second line corresponds to the two-body fermionic decay channel, in the limit where m ψ ≪ m ϕ .Fig. 8 shows the scale factor dependence of the effective mass of the inflaton, inherited from its coherent oscillation (left), and that of the rate R(t) decomposed into the condensate (R ϕ (t)) and particle (R δϕ (t)) contributions (right).In both panels we have multiplied the corresponding variable by its redshift factor prior to ϕ fragmentation.For the effective mass, we have m ϕ (a/a end ) ≃ const.quickly after the beginning of reheating, until a/a end ≃ 160, at the onset of strong rescattering.Due   to the efficient conversion of the inflaton zero-mode into finite momentum particles, m ϕ is further reduced, first in a non-monotonic manner while strong backreaction takes place, and monotonically for a ≳ 300 a end , with m ϕ ∝ (a/a end ) 1.32 , as determined from the lattice data.
In the right panel of Fig. 8 we show the condensate (black) and particle (orange) contributions to the total rate R(t), rescaled by its scale factor dependence and by the amplitude parameter y 2 /8π.Similarly to the effective mass, for a/a end ≲ 160, the rate is solely determined by the anharmonic oscillation of the inflaton, R(t) ∝ m ϕ ρ ϕ ∝ (a/a end ) −5 .However, at later times, the parametric growth of non-vanishing k modes seeps energy from the classical ϕ, reducing accordingly the corresponding rate.For the condensate component, at late times, we find R ϕ ∝ m ϕ ρ ϕ ∝ (a/a end ) −6.6 .On the other hand, following the orange curve, we see that the contribution to R from the fragmented inflaton becomes quickly important, and in fact dominates the dissipation process for a/a end ≳ 200.Since at this stage the comoving number density of δϕ is conserved, we have R δϕ ∝ m 2 ϕ n δϕ ∝ a −5.65 .

Reheating temperatures
With the Boltzmann rate R(t) at hand, we can now estimate the effect of fragmentation on the instantaneous temperature of the primordial plasma during reheating, and its value at the beginning of the domination by the thermal bath, which we denote by T reh .For simplicity we assume the instantaneous thermalization of the inflaton decay products, and for definiteness we assume an effective number of degrees of freedom g reh = 427/4 above the electroweak scale, which is the Standard Model value.

Condensate decay
If we neglect the effect of backreaction, or assume a large coupling constant so that the decay occurs prior to fragmentation, we can use the condensate form for R(t).Fixing for simplicity a end = 1 for   3.3) assuming no fragmentation of the inflaton condensate.For those curves, the end of reheating is marked with a triangle.Continuous curves correspond to the solution of (3.3) accounting for the fragmentation of ϕ.The end of reheating in this case is marked by a star.now, Eq. ( 3.3) can be written as in general, and as for R = R ϕ .In the second equality we have used the fact that, to a good approximation, ρ ϕ a 4 ≃ ρ end ≫ ρ R a 4 up until the end of reheating.Straightforward integration then yields that is, T ∝ a −3/4 during reheating.The end of reheating is determined by the inflaton-radiation equality condition, ρ R = ρ ϕ .Solving with the above expression for ρ R , one obtains [11,12] a reh a end ≃ πρ Fig. 5 shows the evolution of the instantaneous temperature during and after reheating in the condensate approximation, shown as the dashed lines, for y = {10 −1 , 10 −2 , 10 −3 }.This temperature rapidly rises from zero at the end of inflation to reach a maximum T max , which may be approximated as [12] T max ≃ 6 × 10 14 y 1/2 GeV .
Soon after reaching this maximum temperature, the relativistic bath of inflaton decay products continues being populated by these decays, while being redshifted by expansion, following the relation (3.19).Finally, upon reaching the end of reheating, given by (3.20), the thermal plasma simply redshifts, T ∝ a −1 during the radiation dominated epoch.Fig. 6 shows the reheating temperature T reh , as a function of the effective coupling y, corresponding to the orange dashed curve in the pure condensate approximation.In the range of couplings selected the reheating temperature is always well above the lower bound imposed by successful big bang nucleosynthesis (BBN), T BBN ∼ 1 MeV [72].

Fragmented inflaton decay
Taking now into account the preheating phase, we can approximate the total rate R ≃ R δϕ at late times, a/a end ≫ O( 102 ).Parametrizing with γ ϕ ≃ 2.49 × 10 −15 and x ≃ 0.65 (as determined from the lattice, cf.Fig. 8), we can immediately integrate (3.16) to obtain during reheating.This implies that the instantaneous temperature redshifts as T ∝ a −3/4−x/4 ≃ a −0.91 , that is, more rapidly than in the pure-condensate scenario.In this case, reheating ends when and which is the main result of this work.In particular, we note that due to fragmentation, T reh ∝ y 5.68 .Therefore, low reheating temperatures are expected even for comparatively large values of y with respect to the pure condensate scenario.
The effect of the self-resonance of the inflaton on the evolution of the instantaneous temperature of the relativistic plasma can be clearly seen in Fig. 5.Here the continuous curves correspond to the fragmentation scenario.Immediately after the beginning of reheating the temperature follows the condensate result, as the bulk of ρ ϕ is still contained in the coherent component of ϕ.However, after the onset of backreaction, the condensate is depleted and the decay process is dominated by the decay of k ̸ = 0 modes, with a rate that redshifts faster than that of the pure condensate case (cf.Fig. 8).The result is a faster redshift of T , and consequently a reduced reheating temperature, shown in the figure as stars.The difference in reheating temperatures can be better appreciated in Fig. 6.For y ≳ 0.17 the decay of ϕ occurs earlier than the onset of backreaction, and the full result is indistinguishable from the pure condensate scenario.On the other hand, for smaller couplings, T reh is smaller relative to the condensate result, with the difference increasing for decreasing y.The small kinks that may be noted in the T reh vs y curve correspond to the effect of the change in the effective number of relativistic degrees of freedom.We finally note that only for y ≳ 2.8 × 10 −4 can the BBN constrain be averted.

Gravitational waves
In cosmological perturbation theory, tensor metric perturbations are sourced by scalar metric perturbations, starting at second order.This contribution, sensibly negligible prior to fragmentation, would become significant once non-linearity induces sizable mode-mode couplings, efficiently sourcing gravitational waves.In this section we compute the gravitational wave production associated to the large inhomogeneities generated posterior to the end of inflation by the self-fragmentation of the inflaton condensate.Accounting for tensor metric perturbations, the isotropic and homogeneous Friedmann-Robertson-Walker metric is modified to where the transverseness and tracelessness (TT) conditions ∂ i h ij = 0 and h i i = 0 are satisfied.The transverse-traceless component of the anisotropic stress Π TT ij = ∂ i ϕ∂ j ϕ TT appears as a source term in the equation of motion for the tensor metric perturbations 4 where the anisotropic stress projected onto the transverse-traceless component in Fourier space is [41] Π TT ij (k, τ ) = P iℓ ( k)P jm ( k) − with P ij = δ ij − ki kj and ki = k i /k and Π ℓm (k, τ ) = d 3 p (2π) 3/2 p ℓ p m ϕ(p, τ ) ϕ(k − p, τ ) . (4.4)

Simulated spectrum
We use CosmoLattice to simulate the production of gravitational waves for our model.In practice, a TT projector is implemented in a discretized version on the lattice.The non-uniqueness of the definition of such a projector can induce small differences in the GW spectrum.However, it has been argued it only marginally affects the UV part of the GW spectrum [60,73].On the lattice, the GW energy density is computed as a volume average with hij = ah ij .This sum is performed by CosmoLattice, after discretization, over a finite lattice volume V (details can be found in Ref. [61]).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 where we normalized the scale k to the resonant value from the Floquet analysis k ≃ 0.7 m end .ρ rad,0 is the radiation energy density at the present time.The GW energy density estimated from our lattice simulation is shown in Fig. 7, evaluated a different times on the left panel, as a function of the frequency, and integrated over frequency as a function of time on the right panel.From this figure one can see that the GW energy density significantly increases at around τ /τ end ≃ 60 − 80 when scalar inhomogeneities, triggered by the parametric resonance discussed in Sec.2.2, start to increase, as seen in Fig. 3.The GW energy density relative to the critical energy density stabilises at around τ /τ end ≃ 175 when fragmentation is achieved, asymptoting smoothly towards a constant value ρ GW /ρ c ≃ 7.3 × 10 −6 at larger τ /τ end .One of the most striking features of the final GW spectrum, i.e. the red curve in Fig. 7 on the left panel, is the peak structure. 5In the red curve, the dominant peak is located at a frequency f ≃ 1.9 × 10 8 Hz slightly larger than the predicted resonant comoving scale k/m end ≃ 0.7 from the Floquet analysis from Eq. (4.7).From the left panel of Fig. 7, one can see that the initial GW spectra dominant peak was located at a slightly lower frequency and progressively displaced towards higher frequencies, explaining such difference.Disregarding the smaller peak around ∼ 10 8 Hz that seems to appear at the onset of fragmentation (originated by non-linear effects), we can compare the peak structure of the spectrum to those of the scalar spectra and predictions from the Boltzmann approach, cf.Eq. (2.47).A selection of numerical results are presented in Table 1, and are shown graphically in Fig. 8 for the inflaton fluctuation PSD and the gravitational wave spectrum.Even if the peak structure of the scalar spectra are smoothed out by the end of reheating, the structure of the peaks remains imprinted in the tensor spectra.The Boltzmann approach allows for a very decent prediction for the location of the peaks in the scalar spectra.More remarkably, the Boltzmann approach allows to explain the spacing between the first peaks with a precision at the < 10% level.q1 q1 /q 2 q1 /q 3 q1 /q 4 q2 /q 4 q2 /q 5 Lattice (scalar) 0.

Signal and detection prospects
There are essentially two distinct GW signals corresponding to two different physical processes predicted by this model.First, during inflation quantum fluctuations of tensor modes are stretched on macroscopic super-horizon scales.Once such scales re-enter the horizon in a subsequent phase of the universe, tensor perturbations manifest as a stochastic background of gravitational waves.Primordial tensor perturbations, customarily parametrized in terms of the tensor-to-scalar ratio can be expressed for the potential of Eq. (2.2) as [12] r ≃ 12 evaluated at the CMB fiducial scale k * = 0.05 Mpc −1 .Such value for the tensor-to-scalar ratio appears to be close to the sensitivity reach for the upcoming Simons Observatory (SO) r < 6 × 10 −3 [74].Moreover, future missions such as LiteBIRD and CMB Stage-4 (CMB-S4) should reach sensitivities as small as r < 2 × 10 −3 [75] and r < 10 −3 [76] respectively, allowing to disprove or confirm this model.The second GW signal is generated during fragmentation from large inflaton inhomogeneities which we simulated using CosmoLattice.We can estimate the current GW energy density from where we took g reh = g s,reh = 106.75,g s,0 = 3.909 and g 0 = 3.363.The present day normalized GW energy density can thus be expressed as (4.11)We find that for τ /τ end > 300, as illustrated on the right panel of Fig. 7, the produced GW energy density asymptotes during reheating to the value ρ GW /ρ c ≃ 7.3 × 10 −6 .This allows to estimate the contribution to the effective number of non-relativistic degrees of freedom which is far below the current sensitivity from Planck N eff = 2.98 +0.39 −0.38 at 95% confidence level (TT, TE, EE+lowE+lensing+BAO) [46] as well as future CMB missions CMB-S4 [77] and CMB-HD [78] with an expected precision for ∆N eff of 0.06 and 0.027 respectively at 95% confidence level.In order to facilitate comparison with sensitivity prospects for future experiments, we can define the dimensionless characteristic strain as Fig. we compare our estimate of the present day strain as a function of the frequency, in comparison to the sensitivity estimate for high-frequency resonant electromagnetic cavities of Ref. [44].The frequency range of the GW signal emitted from the fragmentation process falls in the appropriate sensitivity range of the proposal of Ref. [44].Remarkably, the reconstitution of the peak structure of such high-frequency GW signal, in addition to measurements of the tensor-to-scalar ratio would help break the ambiguities in determining the realization of inflation.

Summary and conclusions
In this work we have investigated the post-inflationary dynamics of the inflaton field and its decay products, under the assumption of a potential with a quartic minimum.For definiteness, we considered the T-model of inflation, compatible with current constraints on inflation.After inflation, the inflaton condensate coherently oscillates about its minimum.Triggered by parametric resonance effects sourced by the oscillating background, exponential growth of inflaton inhomogeneities result in the fragmentation of this condensate.In this work we explored the role and consequences of such non-linear effects on the post-inflationary history and the successful reheating of the universe.We summarize in the following the key aspects and results of our analysis.
Parametric resonances and fragmentation.First, we characterized the parametric resonant effects triggering fragmentation using a standard Floquet approach.We found that due to the conformal symmetry of the system, only Fourier scales located at k/m end ≃ 0.7 would experience a significant exponential growth.By solving numerically the equation of motion for the inflaton fluctuations at linear order in perturbation theory, we were able to observe a growth of modes with k/m end ≃ 0.7, therefore confirming predictions from the Floquet analysis.To go beyond the linear analysis, we simulated the evolution of the space-time dependent inflaton field configuration with the nonperturbative code CosmoLattice.By computing the occupation number for the inflaton perturbations, we recovered results from the linear approach in the first instants after the end of inflation.We find that at around 5 e-folds after the end of inflation, perturbations enter the nonlinear regime, efficiently backreacting on the inflation condensate, leading to fragmentation.
Before rescatterings fully redistribute the energy density into other modes, the occupation number of inflaton quanta features several peaks located at higher momentum than the resonant k/m end ≃ 0.7 Floquet mode.We estimated analytically the position of these peaks by means of a Boltzmann approach.Notably, despite the validity of this formalism only in the linear regime, we find that it not only accurately predicts the dominant Floquet peak, but it also matches the lattice peak location within a few percent, at least before they are fully washed-out by rescatterings.
Importantly, we found that the redistribution of the inflaton energy density does not fully erase the coherent zero mode.Upon the onset of fragmentation, the energy budget of the universe becomes dominated by a collection of inflaton quanta redshifting as radiation ρ δϕ ∝ a −4 , but a leftover inflaton condensate persists with ρ ϕ ∝ a −5.3 .
Reheating.The main result of this work consists in the exploration of the consequences of fragmentation to successfully achieve reheating; that is, the transition to a universe dominated by radiation in thermal equilibrium.By using a Boltzmann approach, we computed the contributions to the production rates of relativistic fermionic states coming from the oscillating inflaton condensate R ϕ , and also from the fragmentated population of inflaton quanta R δϕ .Prior to backreaction, the rate at which energy is injected into the thermal bath per unit volume per unit time, scales as R ϕ ∝ a −5 .After fragmentation, we find that the condensate rate R ϕ ∝ a −6.6 , and the inflaton particle rate R δϕ ∝ a −5.65 , both redshift faster, due to the rapid decrease in the time-dependent induced inflaton effective mass.As a consequence, reheating is less efficient after fragmentation.The instantaneous temperature falls as T ∝ a −0.91 (compared to T ∝ a −3/4 without backreaction).For fermionic decays, the reheating temperature scales as T reh ∝ y 2 for y ≳ 10 −1 , and as T reh ∝ y 5.68 for y ≲ 10 −1 , where y is the inflaton-matter Yukawa coupling.Requiring the reheating temperature to be larger than the BBN temperature T BBN ∼ MeV necessitates couplings larger than y > 2.8 × 10 −4 .
We must emphasize that these results are valid only for a quartic inflaton minimum.Nevertheless, the lattice+Boltzmann formalism developed in this work is general, and can be applied to a more diverse collection of inflation models, with non-quartic potentials.The precise details of the evolution of energy densities and temperatures will vary on a model-by-model basis, and will be the subject of future follow-up work.
Tensor perturbations.The large mode-mode couplings for the inflaton inhomogeneities inducing fragmentation can also source sizable tensor perturbations.We estimated the production of gravitational waves from a simulation with CosmoLattice.The gravitational wave spectrum features several peaks located at high frequencies f ∼ 10 8 − 10 9 Hz whose precise locations are inherited from the scalar spectrum.The ratios of the location for the various peaks matches relatively well predictions from the Boltzmann approach.We find that the GW energy-density contribution to the effective number of relativistic species ∆N eff ∼ 10 −5 is beyond the reach of any upcoming experiments.However, the frequency range for the GW spectrum appears in sensitivity-reach estimate for resonant electromagnetic cavities [44].

Figure 1 :
Figure 1: Floquet chart for Eq.(2.24).Due to the conformal nature of the model (2.7), a mode that begins inside a resonance band (shown in gray) stays indefinitely in the band, regardless of the expansion of the universe.

Figure 2 :
Figure 2: The particle phase space distribution (PSD) for the inflaton fluctuations for a selection of scale factors, coded by color (left), and the corresponding scale factor dependence of the comoving number density (right).The lattice results are shown in both panels as solid curves.In the left panel, only for the first three values of a we also show the PSD as computed from the linear spectral analysis, as dashed lines.The main Floquet band is shown in light gray.

Figure 3 :
Figure3: Top: Pre-and post-fragmentation evolution of the total inflaton energy density (blue), the energy density of the inhomogeneous component of ϕ (orange) and the energy density of the homogeneous component (green).For ρ δϕ , the linear approximation (2.30) is depicted up to a/a end ≤ 10 2 .At later times the lattice result ρ δϕ = ρ ϕ − ρ ϕ is shown.Bottom: Oscillation-averaged inflaton equation of state parameter.

Figure 4 :
Figure4: Scale factor dependence of the effective inflaton mass (left), and of the effective Boltzmann rate for the relativistic decay products R(t) (right).The total rate is the sum of the condensate contribution (black) and the particle contribution (orange).

Figure 5 :
Figure5: Scale dependence of the instantaneous temperature of the relativistic plasma during and after reheating, for a selection of effective inflaton-fermion couplings.Shown as dashed curves is the numerical integration of Eq. (3.3) assuming no fragmentation of the inflaton condensate.For those curves, the end of reheating is marked with a triangle.Continuous curves correspond to the solution of (3.3) accounting for the fragmentation of ϕ.The end of reheating in this case is marked by a star.

Figure 6 :
Figure6: Reheating temperature as a function of the inflaton-fermion coupling, accounting for ϕ fragmentation (continuous), and without accounting for fragmentation (dashed).Shown in gray is the region forbidden from big bang nucleosynthesis considerations.

Figure 7 :
Figure 7: Left: Differential gravitational wave energy density per logarithmic (present day) frequency interval, normalized to the critical density, evaluated at selected conformal time.The color code corresponds to the evaluation time of the right panel.Right: Integrated differential GW energy density evolution with respect to conformal time, normalized to the critical energy density.

Figure 8 :
Figure 8: Comparison between the peak structure of the scalar phase space distribution at the onset of fragmentation (left), and of the final gravitational wave spectrum (right), with the predictions of the Boltzmann approach (vertical gray lines).The Boltzmann lines are only indicative of the value of the comoving momenta, the ordinate values are matched with the corresponding black curve.
[41] P .The GW frequency f at the present epoch can be related to the comoving Fourier scale k via[41]

Table 1 :
Comparison of the peak structure of the scalar and tensor spectra where qi corresponds to the location of the i th peak in terms of the rescaled momentum q ≡ k/m end .The Boltzmann predictions are based on Eq. (2.46) and Eq.(2.47).
[44]lated amplitude of the present-day GW signal amplitude as a function of the frequency (solid red line) and sensitivity estimate for resonant electromagnetic cavities of Ref.[44](orange dashed line).Such gravitational wave signal would induce a deviation to the effective number of relativistic species ∆N eff ≡ N eff − N SM eff with N SM eff = 3.046.This contribution can be expressed as