Impact of Dark Sector Preheating on CMB Observables

The prediction of a nearly scale-invariant spectrum of curvature and tensor fluctuations is among the main features of cosmic inflation. The current measurements of the primordial fluctuations in the cosmic microwave background (CMB) provide tight constraints on the amplitude of the scalar and tensor spectra, and the scalar tilt. However, the precise connection between these observables and a given inflationary model, depends on the expansion history between the end of inflation and the beginning of the radiation dominated era, which corresponds to the reheating epoch. This mapping between horizon exit and reentry of fluctuations, parametrized by the number of $e$-folds $N_*$, can therefore be affected by the presence of a transient epoch of non-perturbative particle production during reheating ({\em preheating}). Using a combination of perturbative and lattice computations, we quantify the impact of preheating in a non-equilibrated dark matter sector on the CMB observables, under the assumption of a simultaneous perturbative decay of the inflaton into Standard Model particles. Combined with structure formation constraints, this allows us to impose stringent bounds on the post-inflationary reheating temperature.


Introduction
The successful solution of the classical problems of the Standard Big Bang cosmological scenario (the flatness, horizon, and monopole problems), and the prediction of a near scale-invariant spectrum of primordial fluctuations, are among the main features of the inflationary paradigm [1][2][3][4].In its arguably simplest realization, inflation is realized by the slow-roll of a real, scalar field known as the inflaton, down a flat potential.Quantum fluctuations in the inflaton field and the spacetime metric are also stretched by the expansion beyond the horizon scale, turned into curvature and tensor fluctuations.Their primordial spectra leave an imprint in the form of temperature and polarization anisotropies of the cosmic microwave background (CMB).Their measurement and characterization by WMAP, Planck, and BICEP/Keck provide now a stringent bound at the pivot scale k * = 0.05 Mpc −1 , on the ratio of the tensor power spectrum relative to the scalar one, r < 0.036, at the 95% CL.Moreover, the tilt of the scalar spectrum lies in the range 0.958 < n s < 0.975 for r = 0.004 at the 95% CL [5].These values for the CMB observables can be comfortably accommodated in scenarios with asymptotically flat inflaton potentials, among which the R + R 2 Starobinsky model can be singled out in particular [6].
An additional feature of slow-roll inflation is the automatic graceful exit from accelerated expansion, as the field reaches the minimum of its potential, and commences to coherently oscillate about it.The dissipation of the energy stored in these oscillations into a hot, thermal plasma of relativistic particles reheats the universe.Generically, this reheating not only populates the Standard Model (SM) degrees of freedom, but also the degrees of freedom that make up the dark sector of the universe.At the perturbative level, the dissipation process is modeled as the quantum mechanical decay of inflaton quanta into elementary particles, with a decay rate determined by the adiabatic transition amplitude, ignoring the short time-scale oscillation of the inflaton field [7][8][9][10][11][12][13][14].
There exist, however, scenarios in which the high frequency oscillation of the inflaton is the primary driver for particle production.In particular, integer-spin light fields coupled to the inflaton can undergo an exponential growth of their fluctuations, due to the resonant amplification of their field modes originated from a nonadiabatic change in their effective masses.This nonperturbative particle production mechanism is known as preheating, and can significantly alter the post-inflationary expansion history of the universe [8,[15][16][17][18][19].When the amplification of fluctuations is strong, it can backreact into the classical, coherent inflaton condensate via mode-mode couplings, sourcing large gradients and eventually leading to its fragmentation [11,[20][21][22][23][24].Although this non-perturbative production of particles does not lead in general to a complete reheating and a radiation dominated era, it can play a role in the efficiency of the decay of the inflaton into light degrees of freedom in the early and late stages of reheating, modifying the duration of the reheating process [25].Importantly, despite the fact that super-horizon fluctuations are generically impervious to these complex dynamics [26][27][28], the precise connection between a given observable scale and its horizon-exit during inflation depends on the duration of reheating and the corresponding change in the total energy density of the universe [29][30][31][32][33][34][35][36][37][38][39][40].Therefore, the presence of preheating, on top of a dissipative, perturbative decay of the inflaton, can add to the uncertainty that reheating embodies on the connection between a given inflationary model and its corresponding CMB predictions [41][42][43][44][45].
As studied by one of the authors in previous works, preheating provides an efficient out-ofequilibrium mechanism for generating the relic dark matter abundance of the universe [23,46,47].In this paper we explore the impact of dark sector preheating on the CMB observables n s and r, in the case when the inflaton oscillates about a quadratic minimum of its potential.For definiteness we consider a scalar dark matter (DM) candidate, minimally coupled to gravity, which interacts directly only with a Starobinsky-like inflaton sector.We numerically track the evolution during reheating of the inflaton field, the dark matter fluctuations, and the relativistic thermal plasma which contains the Standard Model particles, which we assume to be produced perturbatively.For large inflaton-DM couplings, we track the instantaneous energy density of the plasma before, during, and after the backreaction regime, during which the oscillating inflaton condensate is fragmented in favor of an inhomogeneous collection of free particles.As a result, we obtain semi-analytical estimates of the change in the number of e-folds, the scalar tilt, the tensor-to-scalar ratio, and the DM relic abundance, in terms of the depletion of the comoving inflaton energy density due to backreaction.From these parameters, we derive strong constraints on the allowed couplings between the inflaton and the visible sector fields.This paper is organized as follows.In Section 2 we revisit the analytical and numerical approximations necessary to estimate to high accuracy the CMB observables for the Starobinsky model of inflation, with emphasis on their dependence on the reheating epoch.Section 3 is devoted to the discussion of the parametric excitation of dark degrees of freedom due to their coupling to the coherently oscillating inflaton during reheating.We present a numerical exploration of the nonlinear backreaction regime, and the resulting fragmentation of the inflaton condensate.Our main results are discussed in Section 4. In Section 4.1 we find a parametrization of the effect of strong DM preheating on the evolution of the inflaton and radiation energy densities during reheating.In Section 4.2 we compute the DM relic abundance for strong preheating, and determine the effect of the backreaction onto the inflaton field.Section 4.3 is devoted to the precise computation of the CMB observables, and to mapping the current observational constraints onto the effective strength of the inflaton-matter couplings.We present a summary and our conclusions in Section 5.

Inflationary dynamics
The field content of our study consists of an inflaton field ϕ, a sterile scalar dark matter field χ, and the Standard Model degrees of freedom, relativistic at early times, to which we will collectively refer in terms of their total energy density ρ R .We write the corresponding action as Here g ≡ det(g µν ) denotes the metric determinant of a flat Friedmann-Robertson-Walker metric with scale factor a, m χ is the DM mass, σ is the dimensionless inflaton-DM coupling, L SM contains all the SM terms, including their coupling to the inflaton, and V (ϕ) denotes the inflaton potential.For definiteness we choose this potential as the Starobinsky potential, where GeV denotes the reduced Planck mass.In its original inception, this model arises as the effective scalar-tensor theory in the Einstein frame of the R 2 generalization of the Einstein-Hilbert action [6].The model also naturally appears in no-scale supergravity inflation models, in which a singularity in the kinetic term of the inflaton is manifested as the asymptotically flat potential (2.2) once the field is translated into its canonically normalized counterpart [48][49][50][51][52].
Disregarding the coupling of the inflaton to other fields during inflation, variation of the action (2.1) with respect to the homogeneous inflaton field and metric yields the equations of motion where dots denote derivatives with respect to cosmic time, and H = ȧ/a is the Hubble parameter.The solution of this system describes the slow-roll of the inflaton field as it rolls down the potential (2.2), and the quasi-exponential growth of the scale factor a.Although the numerical solution of Eqs.(2.3)-(2.4) is a straightforward exercise, analytical approximations for the amount of expansion and the inflationary observables can be obtained via the slow-roll approximation, which we review for the Starobinsky potential below.

Slow-roll approximation
Near the minimum, for ϕ ≪ M P , the potential is quadratic, V (ϕ) ≃ 1 2 m 2 ϕ ϕ 2 , and thus the parameter m ϕ corresponds to the inflaton mass.Its value is determined by the amplitude of the curvature power spectrum, A S ≃ 2.1 × 10 −9 , evaluated at the horizon exit time of the Planck pivot scale k P = 0.05 Mpc −1 [53,54].The connection can be determined by means of the slow-roll approximation.This approximation is valid when the slow-roll parameters are both ≪ 1.In terms of these parameters, the amplitude of the scalar spectrum can be approximated as where the star denotes evaluation at the horizon crossing time of the scale k * .The field value at this time can be computed in terms of the number of e-folds to the end of inflation yielding ϕ end where ϕ end ≃ 0.615 M P , (2.9) is the inflaton field value at the end of inflation, when ä = 0, or equivalently φ2 end = V (ϕ end ) [36,38].Substitution into (2.6)gives the following mass normalization, (2.10) The scalar tilt n s , and the tensor-to-scalar ratio r, can be determined to leading order from the slow-roll parameters as follows [55], (2.12) The number of e-folds to the end of inflation are in turn related to the post-inflationary expansion history via the expression [29,31] 1/4 43 11 where H 0 = 67.36km s −1 Mpc −1 [53], T 0 = 2.7255 K [56] and a 0 = 1 are the present Hubble parameter, CMB temperature and scale factor, respectively.The energy density of the universe at the end of inflation is denoted by ρ end , while g reh are the number of degrees of freedom at the end of reheating.At temperatures above the electroweak phase transition we will consider for definiteness the SM value g reh = 427/4.R rad denotes the reheating parameter [31] R rad ≡ a end a rad where ρ rad and a rad are the energy density and scale factor at the onset of radiation domination after the end of reheating, 1 and a end is the scale factor at the end of inflation.For the Starobinsky potential, substitution of (2.8), (2.9) and (2.10) results in the following non-algebraic equation, We note that for the Planck pivot scale, and SM degrees of freedom after inflation, the previous expression results in the upper bound N * ≲ 55.74, corresponding to instantaneous reheating, R rad = 1.If reheating is not instantaneous, and can be modeled as the perturbative decay of the inflaton, the reheating parameter can be approximated as [38] ln where denotes the two-body inflaton decay rate, with y an effective Yukawa-like coupling.

Exact CMB observables
An accurate determination of the inflationary observables requires the numerical computation of the primordial scalar and tensor power spectra.The necessity for this is two-fold.Firstly, Eqs.(2.11) and (2.12) are accurate only to leading order in the slow roll approximation.Secondly, and more importantly, curvature and tensor spectra do not freeze-out in value immediately upon horizon crossing, and this delay leads to a shift between the analytical and numerical relation between (n s , r) and N * [38].Scalar fluctuations can be characterized in terms of the gauge-invariant Mukhanov-Sasaki variable, which in the Newtonian gauge is defined in terms of the inflaton fluctuation δϕ and the scalar metric perturbation Ψ as [57,58] where the dot denotes a derivative with respect to cosmic time, and â † k and âk denote creation and annihilation operators, respectively, which obey the canonical commutation relations The mode functions Q k satisfy the equation of motion with Bunch-Davies initial condition Q k≫aH = e −ikτ /a √ 2k, being dτ = dt/a the conformal time.Upon solution, the curvature perturbation is determined as with power spectrum that evaluates to in terms of the Q mode functions.The amplitude and tilt are evaluated from their definitions, not at horizon crossing, but at the end of inflation, ensuring that their freeze-out values have been reached,  In the case of tensor perturbations, which are automatically gauge-invariant, we use the standard transverse, traceless perturbation where γ = +, × denotes the two polarization modes, ϵ γ ii = k i ϵ γ ij = 0 and ϵ γ ij ϵ γ ′ ij = 2δ γγ ′ , and bk,γ , b † k,γ denote the creation and annihilation operators, for which Without anisotropic stress sources at linear order, the equation of motion satisfied by the tensor mode functions is given by The corresponding power spectrum is defined as from which the tensor-to-scalar ratio is computed as follows, Fig. 1 shows the form of the scalar and tensor power spectra for the Starobinsky model, from the numerical integration of (2.19) and (2.25).There we assume that the Planck pivot scale left the horizon 55 e-folds before the end of inflation.As is shown, choosing the inflaton mass as m ϕ ≃ 3 × 10 13 GeV we recover the measured amplitude of the curvature power spectrum at k P .For comparison purposes we emphasize the difference in the spectra amplitudes at the WMAP pivot scale, k W = 0.002 Mpc −1 , at which the constraints on the tensor-to-scalar ratio are commonly presented [5,53,54].Fig. 2 shows in turn a comparison between the numerical and analytical approximations to the CMB observables n s , r, here as functions of the number of e-folds.We note a relatively large deviation between the numerical result (black, continuous) and the slow-roll result (2.11) and (2.12) (gray, dotted).The main culprit behind this shift is the non-instantaneous freeze-out of the power spectra upon horizon crossing; the spectra continue evolving until k/aH ≪ 1.In order to make this manifest, we also show in Fig. 2 the red, dashed curves, which correspond to the evaluation of the slow-roll approximation with a time-delay of 2.5 e-folds.In this case, the agreement between the exact results, evaluated at the end of inflation, and the analytical estimates is significantly improved.

Dark sector preheating
We now turn to the description of the dynamics of reheating.Our discussion is based on the assumption of a perturbative decay of the inflaton into SM degrees of freedom, e.g. through the decay of the inflaton into fermions, and the non-perturbative production of scalar dark matter.In this Section we describe the resonant production of dark matter particles, in the linear and non-linear regimes, using a combination of spectral and lattice methods.

Resonant production of dark matter
Starting from the field action (2.1), and upon variation with respect to the DM field χ, we obtain the following equation of motion, This equation can be reformulated in a quantization-friendly form by switching to conformal time and introducing the re-scaled field Here Âp and Â † p denote the annihilation and creation operators, respectively, which satisfy the standard commutation relations [ Âp , Â † The canonical commutation relations between the field operator X and its conjugate momentum are in turn fulfilled if the Wronskian condition X p X * ′ p −X * p X ′ p = i is satisfied.In terms of the X mode functions, Eq. (3.1) takes the form where primes denote differentiation with respect to conformal time, and where the effective angular frequency is defined as The initial conditions for the X p are taken to be the Bunch-Davies vacuum mode functions, At any later time, the UV-regular number and energy densities of the field χ can be computed as follows [19,23] ) where denotes the phase space distribution (PSD) of the field χ, identical to the particle occupation number of the DM.The effective frequency (3.5) can be equivalently rewritten in terms of the instantaneous Hubble parameter as ω 2 p = p 2 + a 2 m 2 eff , where and w = P/ρ denotes the total equation-of-state parameter [47].For sufficiently small bare mass m χ , coupling σ, and momentum p, the last term in (3.10) might dominate during inflation (w ≃ −1), leading to the tachyonic excitation of super-horizon modes.Assuming a DM mass smaller than the Hubble parameter during inflation, m χ ≪ H, the tachyonic instability would occur if In other words, for where the non-perturbative growth of low-momentum modes during inflation can be safely ignored.For the Starobinsky model, evaluation of (3.12) at the end of inflation leads to σ/λ ≳ 0.5.In what follows, we will assume that the effective coupling σ/λ ≫ 1, and therefore we will disregard the evolution of the χ mode functions prior to the beginning of inflation.
In the large coupling regime, the dominant term in the effective mass corresponds to σϕ 2 , which is quasi-periodic in time, leading to the parametric excitation of momentum modes.The solution to the equation of motion (2.3) near the minimum of the potential is approximately given by with ϕ 0 ≃ ϕ end m ϕ t denoting a decreasing envelope.It is then convenient to introduce the following variables, where denotes the dimensionless comoving momentum relative to the end of inflation, rescaled by with respect to the inflaton mass.Eq. (3.4) takes then the Mathieu form which presents parametric resonance.More precisely, Floquet's theorem states that the solutions of this equation have the form [59] x q (z) = e µqz g 1 (z) + e −µqz g 2 (z) .
Here g 1 (z) and g 2 (z) are periodic functions, and µ q is called the Floquet exponent.If Re µ q ̸ = 0 exponentially growing solutions exist.The Floquet exponents, and the corresponding resonance bands as a function of A q and κ can be found following the eigenvalue method described in detail in [11].The left panels of Fig. 3 show the structure of the resonance bands of the Mathieu equation (3.18), with a color coded Floquet exponent.Since the parameters A q and κ depend on redshifting quantities, each momentum mode flows on this chart, from an initial point determined by ϕ end and the effective coupling σ/λ, to the origin.The figure shows the corresponding flows for a few selected (re-scaled) momenta q for σ/λ = 10 2 , 10 4 .For the first case, low-momentum modes only cross a couple of resonance bands, which is manifest in the behavior of the energy density of the DM field, shown in the right panel.The resonant growth of ρ χ is not accumulated over several oscillations, and the χ field enters the pure redshift regime, with ρ χ ∝ a −4 , for a/a end ≳ 4. On the other hand, for σ/λ = 10 4 , each low momentum mode can cross a large number of resonance bands, and the resonance can efficiently accumulate.As shown in the right panel, this results in the quasi-exponential growth of the energy density of the DM, which becomes quickly comparable with ρ ϕ .The linear description in terms of Eq. (3.4) breaks down and the system then enters the backreaction regime, which we describe below.

Backreaction and fragmentation
As demonstrated above in Fig. 3, for large effective couplings σ/λ ≫ 1, the resonant growth of the χ mode functions can take them into the backreaction regime.The energy density of χ becomes comparable to the energy density of the inflaton, and the DM can no longer be treated as an spectator field.More importantly, the explosive production of DM can disrupt the homogeneous inflaton condensate by means of the rescattering of χ particles into ϕ, i.e. the production of inflaton particles.Moreover, the free DM particles can also scatter inflaton fluctuations away from the condensate.At the level of fields, these dynamics are manifested as the growth of spatial gradients of the inflaton field, eventually leading to the "fragmentation" of the homogeneous condensate [11,[20][21][22][23].
In terms of the field mode functions, these effects manifest as mode-mode couplings, since now the inflaton field in (3.1) is to be replaced by ϕ(t) → ϕ(t) + δϕ(t, x), where the fluctuation δϕ is also quantized.This of course leads to a set of non-linear Heisenberg equations of motion for the ϕ, χ quantum fields, which is challenging to solve.For this reason, spectral methods are not the tool of choice in the backreaction regime.Instead, we resort to the solution of the classical equations of motion for the coupled fields.At large couplings the occupation numbers for the resonantly excited momentum modes are f χ,ϕ ≫ 1, justifying this classical approach.The system of non-linear partial differential equations with is solved by means of finite-difference techniques on a spatial lattice.Among the variety of publicly available codes [60][61][62][63], we favor CosmoLattice [64,65] due to its ease of use.The lower right panel of Fig. 3 shows precisely the lattice solution of the system (3.20)-(3.23)for σ/λ = 10 4 .For definiteness we have associated the presence of backreaction with the condition ρ χ ≳ 0.1 ρ ϕ [23,46].Since this system is not dissipative, the inflaton and the DM exchange energy when the backreaction regime is reached, entering a state of quasi-equilibrium until the expansion of the universe overcomes the rescattering rates, and the system enters the pure redshift regime, with ρ ϕ ∝ a −3 and ρ χ ∝ a −4 .

Distributions and densities
By means of a combination of spectral and lattice methods we can track the time evolution of energy and number densities of fields, and of their corresponding PSDs.The top panels of Fig. 4 show the PSDs of the DM and the inflaton as functions of the scale factor, for the coupling σ/λ = 10 4 (and σ/λ = 10 5 for the bottom panels). 2 Focusing on χ, at early times, a/a end ≲ 9 (a/a end ≲ 5), the growth of the distribution follows the spectral result as given by (3.9).The distribution clearly shows a peak pattern produced by the resonant growth owing to the coupling of χ to the oscillating inflaton.This growth "in bursts" can also be appreciated in the DM energy and number densities, which are displayed in Figs. 5 and 6.At later times, the backreaction regime is reached, and the efficient energy exchange between the dark and inflaton sectors rapidly erases the resonance peaks of the distribution.The distribution quickly evolves towards the UV, with two plateau-like regions owing to Kolmogorov-like turbulence [66,67], and a quasi-thermal exponential tail.Relativistic modes come to dominate the energy density of χ, and eventually the pure redshift regime ρ χ ∝ a −4 is reached.This signals the end of the backreaction regime, and the return to dominance of ϕ.The PSD freezes, and so does the comoving number density, as shown in the right panels of Figs. 5 and 6.
In the absence of DM interactions other than its coupling to ϕ, the number density produced during reheating is only redshifted to become the present DM relic abundance.We postpone this discussion to Section 4.2.Focusing now on the inflaton ϕ, the right panels of Fig. 4 show how the production of nonzero modes of ϕ (i.e.particles) is induced by rescattering.Even before backreaction is reached we have f ϕ ≫ 1 for low momentum modes.Initially, the production of non-relativistic modes is most efficient, until a maximum in the PSD is reached when backreaction is reached.Posterior to this, the cascading of energy toward the UV is efficiently mediated by rescattering, until the PSD freezes when ρ χ ≪ ρ ϕ .The comoving number density of inflaton quanta remains also frozen, as it can be observed in the right panels of Figs. 5 and 6.The depletion of the inflaton zero-mode in favor of q ̸ = 0 modes can be more explicitly visualized in the upper left panels of Figs. 5 and 6.In them, three energy densities for the inflaton have been identified.The total energy density ρ ϕ is computed from the spatial average of the full energy momentum tensor of ϕ, On the other hand, we identify the condensate component of the energy density of ϕ as the energy density of the spatially averaged field [23,25,68], The energy density of inflaton fluctuations can then simply be defined as whose value is consistent at early times with the spectral definition of Eq. (3.8), replacing χ by ϕ [25].In terms of these "components" of ρ ϕ it is easy to see the onset of backreaction in the oscillating condensate.Non-zero modes are efficiently populated by rescatterings, draining energy density from φ until ρ δϕ > ρ φ.At this point, the inflaton can be considered as fragmented in favor of free quanta.As the bottom panels of Figs. 5 and 6 also show, the rapid population of relativistic ϕ and χ modes drives the oscillation-averaged equation-of-state parameter ⟨w⟩ close to its radiationdomination value.Importantly, the depletion of ρ φ is rapid and complete.After fragmentation the inflaton is fully replaced by free particles, which eventually are redshifted by expansion and, in the absence of other interactions, lead again to a matter dominated era, with ρ ϕ ≃ ρ δϕ ≃ m ϕ n ϕ .
The end of backreaction, and the subsequent dominance of nonrelativistic free inflaton quanta, indicate that reheating cannot be completed by dark sector preheating.The couplings contained in L SM in (2.1) are crucial to fully dissipate ϕ and reheat the universe.In the next Section we study the decay of the fragmented inflaton, and the effect that the nonlinear dynamics during preheating have on the duration of reheating, the dark matter abundance and the inflationary CMB observables.

Signatures of preheating
Having characterized the effect of dark sector preheating on the evolution of the background metric and inflaton energy densities, we now proceed to quantify its effect on the reheating process, understood here as the dissipation of the inflaton energy density into SM relativistic degrees of freedom.We begin by reviewing the Boltzmann formalism for perturbative reheating, to continue with the determination of DM relic abundances, and the effects of preheating in CMB observables.

Reheating after inflaton fragmentation
As we briefly mentioned in the Introduction, we will assume that the decay of the inflaton into the visible sector is driven by the perturbative decay of the inflaton quanta.For simplicity, we will model this process as a two-body decay into a pair of fermions f , with effective Yukawa coupling y, If the mass of the final state particles can be neglected, and assuming that quantum statistics can be disregarded, the Boltzmann equation that determines the population of the primordial radiation plasma with energy density ρ R is given by [25] ρR where Γ ϕ is the vacuum inflaton decay rate (2.17).The first term in the right hand side corresponds to the contribution from the dissipation of the coherent inflaton condensate, with w ϕ = 0.The second term represents the decay of the free particles that compose the fragmented inflaton after backreaction.It is worth recalling here that prior to fragmentation (or in the absence of it), ρ ϕ ≃ ρ φ, while after fragmentation, ρ ϕ ≃ m ϕ n δϕ , as shown in Fig. 5.
Disregarding the coupling of ϕ to χ, which is valid well before or well after the backreaction regime with w ≃ 0, Eq. (4.2) is complemented by the system of equations from which the dynamics of the background quantities can be fully determined.The end of reheating is defined as the crossover time t reh , when ρ ϕ (t reh ) = ρ R (t reh ).In turn, the reheating temperature T reh is defined from the thermodynamic relation evaluated at t = t reh .In the absence of dark sector preheating, Eqs.(4.2)-(4.4)can be solved approximately during reheating.For t ≪ t reh , ρ ϕ ≃ ρ end (a end /a) 3 , and (4.2) can be rewritten as   Inflaton and radiation energy densities after the end of inflation, in the presence of dark sector preheating (continuous) and in its absence (dashed).For the continuous curves, the full inflaton energy density is shown as the blue line, while the contribution of non-relativistic modes is shown in purple.Their difference is highlighted in the inset.The horizontal dotted line corresponds to the energy density at the end of reheating, with solution [14] ρ R (a) ≃ 2 5 3M This expression implies that the maximum temperature of the primordial plasma, assuming instantaneous thermalization, is given by In turn, the reheating temperature can be approximated by , at a reh a end pert ≃ 25 12 As is well know, T reh depends only on the couplings and number of degrees of freedom of the theory if a reh ≫ a end .Therefore, we do not expect this temperature to be dependent on the details of preheating, as long as the complete decay of the inflaton does not occur during fragmentation.
In the presence of dark sector preheating, the predictions for the duration and final temperature of reheating are modified, as is shown in Fig. 7 for σ/λ = 10 4 .The perturbative approximation discussed above is shown as the dashed blue and red lines, labeled as ρ (pert) ϕ and ρ (pert) R , respectively.On the other hand, the energy densities of the inflaton and the radiation in the presence of non-perturbative DM production are shown as the continuous blue and red curves, respectively.We observe that, prior to the onset of backreaction, for a/a end ≲ 8, the exact energy densities follow closely the perturbative approximation, with a small deviation originating from the transition from the quasi-de Sitter epoch to matter domination within the first oscillation of ϕ.Notably, the maximum of ρ R is essentially identical in the two cases, and therefore the approximation (4.8) remains valid in the presence of a strong ϕ-χ coupling.
As it can be appreciated in Fig. 7, and also Fig. 5, χ backreaction is important for 8 ≲ a/a end ≲ 100.During this time, a significant fraction of the energy density of the inflaton condensate is transferred to DM and free inflaton particles.For a/a end ≳ 10, ρ δϕ > ρ φ, and inflaton particles make up almost the entirety of ρ ϕ .Nevertheless, only a fraction of them are non-relativistic, due to the UV-cascading induced by mode-mode couplings (see Fig. 4).The result is the hierarchy ρ (pert) ϕ > ρ ϕ > m ϕ n δϕ ≫ ρ φ, which suppresses the effective particle production rate, as per (4.3).Therefore, a reduction in the radiation energy density is expected, and it can be appreciated in the continuous red curve.In this regime, the time-dependence of all background quantities must be determined numerically, from the lattice results and the continuity equation (4.2).
Backreaction ends when rescattering is strongly suppressed by the expansion of the universe.For σ/λ = 10 4 , this corresponds to a/a end ≳ 100, as shown in Fig. 7.In this regime, ρ ϕ ≃ ρ δϕ ≃ m ϕ n δϕ .The energy density of DM is no longer comparable to ρ ϕ , and the main particle production process is now the perturbative decay into radiation.Therefore, the system (4.2)-(4.4)with ρ φ ≃ 0 accurately describes the evolution of the inflaton-radiation mixture even beyond the end of reheating.Let us denote by the loss in the comoving inflaton energy density between the end of inflation (a end ) and the end of backreaction (a m ), due to χ-preheating.Then, we can modify (4.9) to estimate the change in the reheating temperature and the duration of reheating due to the depletion in ρ ϕ .We note that the reheating temperature, being independent of ρ end is insensitive to preheating, while the duration is modified to a reh a end preh ≃ 25 12 For the parameters in Fig. 7, our analytical estimates give ρ reh ≃ 1.21 × 10 −25 M 4 P , a reh /a end ≃ 6.1 × 10 4 (no preheating) and a reh /a end ≃ 3.7 × 10 4 (χ preheating).Numerically, we find ρ reh ≃ 1.26×10 −25 M 4 P , a reh /a end ≃ 4.3×10 4 (no preheating) and a reh /a end ≃ 2.6×10 4 (χ preheating).The analytical estimates for a reh are shifted by a factor of ∼ 1.4 since we have neglected the exponential decay of ρ ϕ at the end of reheating [36,69].Nevertheless, there is a few percent level agreement for ρ reh and a sub-percent agreement for the shift factor ∆ 1/3 ≃ 0.60 for a reh in the presence of resonant DM production.
With the duration of reheating, and the dilution of the energy density, we can now evaluate the shift in the number of e-folds for a given pivot scale k * .Substitution of (4.11) in (2.13), with a rad ∝ a reh , gives in excellent agreement with numerical results, shown in the first panel of Fig. 8.The second panel of the figure demonstrates the independence of T reh with respect to dark sector preheating, for perturbative Yukawa couplings.It is worth noting that Eq. (4.12) has the same functional dependence for the shift in N * in the case of a late-time injection of entropy by a factor ∆, as is the case for embeddings of Starobinsky inflation in flipped SU(5) grand unified theories [38,70].

Dark matter relic abundance
In the absence of DM interactions other than its coupling to ϕ, the number density produced during reheating is only redshifted to become the present DM relic abundance.Assuming the χ field is nonrelativistic at the present time, the dark matter relic abundance can be computed as where ρ c = 3H 2 0 M 2 P is the present critical density of the universe, and n χ (a m )(a m /a end ) 3 is the comoving DM number density, evaluated after backreaction and therefore constant in time.Following a standard computation, the amount of expansion between the end of inflation and today can be approximated as [29,31,46]   where in arriving to (4.This expression for the primordial DM closure fraction is an improvement over the computation presented previously in [46], as it takes into account the more prompt reheating due to DM backreaction.Fig. 9 shows the reheating temperature necessary to saturate the observed value of Ω χ , as a function of the DM mass.The continuous lines show the (inverse) relation accounting for the shift in a reh , ∆ ̸ = 1, while the dashed curves ignore this effect.As expected, the effect of ∆ is to enhance the primordial abundance given a fixed m χ and T reh .Therefore, lower masses/temperatures are required compared to the naive result with ∆ = 1.This shift corresponds simply to a reduction by a factor of ∆ in T reh for a given m χ , or vice-versa, given that the logarithmic correction to N * is negligible in comparison (see Eq. (4.12)).We note that high reheating temperatures require small (sub-keV) masses.In such cases, the χ field can behave at late times as non-cold dark matter, resulting in a suppression of clustered structure overdensities on galactic scales.This suppression is absent in Lyman-α measurements of the matter power spectrum above scales k ≃ 15h Mpc −1 .For warm dark matter candidates, this constraint can be translated to a lower bound in the mass, m WDM = (1.9− 5.3) keV [71][72][73][74][75][76][77].In the present case of non-thermal, resonant DM production, this lower bound can be translated upon matching the equation-of-state parameter of χ with the equationof-state parameter of a thermal candidate with mass m WDM .The full details of this procedure can be found in [78].In Refs.[46,47] this mapping has been applied to the model at hand (2.1).The result is a quasi-universal lower bound in the backreaction regime, of m χ ≳ 30 eV.The shaded gray region in Fig. 9 corresponds to the parameter space ruled out by structure formation considerations.It is worth noting that the shortened reheating due to backreaction lowers the allowed range of temperatures.The maximum reheating temperatures correspond to For more details on the phenomenology of DM production from preheating, including relic abundance, structure formation, isocurvature and reheating constraints, in the weak and strong coupling regimes, see [46,47].

CMB observables
After having determined the effect of DM preheating on the number of e-folds after horizon crossing, we are now in position to compute the effect in the CMB observables n s and r.A straigthforward substitution of the numerical results in Fig. 8 onto the exact spectral parameters in Fig. 2 produces the curves seen in Fig. 10.The left panel shows the dependence of the spectral tilt on the Yukawa coupling y in the absence of preheating (solid, red), and for the DM-inflaton couplings σ/λ = {10 4 , 10 5 , 10 6 } (dashed, dot-dashed and dotted, respectively).In this same panel, for reference, we show the 1σ and 2σ Planck+BK18+BAO CL contours, marginalized at r 0.002 ≃ 0.004 [5].As it is clear, the result is an increase in the tilt of the scalar power spectrum, larger for stronger ϕ-χ couplings.This shift in n s is present at all y < 1. Combining expressions (2.11) and (4.12), and the correction to the slow-roll prediction shown in Fig. 2, we can parametrize the effect of DM preheating on n s as This result is in excellent agreement with the numerical data, indistinguishable from it at the subpercent level.The effect on n s from non-linear dynamics is therefore always small, δn s ≲ 10 −3 .Such a change is smaller than the forecasted sensitivity of future CMB experiments such as CMB-S4 [79], the Simons Observatory [80] or LiteBIRD [81,82], although it could in principle be probed by EUCLID [83] or SKA [84,85].Nevertheless, the shift in n s is important for constraining the allowed inflaton-SM couplings.Specifically, we find that the minimum coupling allowed by observations can differ by more than an order of magnitude at the 68% and 95% CL, as shown in Table 1.A similar result applies for the minimum reheating temperature, also shown in Table 1.At σ/λ = 0 the bounds presented here are consistent with those in [38], accounting for a difference in the number of degrees of freedom g reh at high temperatures (SM vs MSSM), and a slight modification in the numerical method.
The right panel of Fig. 10 shows the numerically computed values for the tensor-to-scalar ratio r as a function of the Yukawa coupling, in the absence and presence of preheating, analogously to the n s panel.Here we do not depict the Planck+BICEP/Keck constraints, as for all couplings the tensor-to-scalar ratio is below the 1σ upper bound r ≲ 0.035 [5].We note in this case that the effect of preheating is a reduction in the tensor amplitude relative to the scalar spectrum.This effect can  be justified upon combining the slow-roll approximation (2.12), together with (4.12) and the shift N * → N * + 2.5 which corrects the slow-roll result (see Fig. 2), r ≃ r (pert) + 8 ln ∆ (N * + 2.5) 3 . (4.20) Similarly to n s , the effect on r is far below the sensitivity threshold of the next generation of CMB observatories.Moreover, the constraints on the inflaton couplings to other fields are largely independent of r for small values of it.Therefore, the effect of preheating on gravitational waves sourced around CMB frequencies can be neglected for plateau-like potentials.Nevertheless, as is well known, the non-linear growth of scalar inhomogeneities during preheating can efficiently source gravitational waves at high frequencies, f ∼ 10 8 Hz [86][87][88][89][90][91][92].At CMB scales, the main effect would correspond to an enhancement of the effective number of relativistic species, ∆N eff = N eff − N SM eff , with N SM eff = 3.046.However, the resulting enhancement is typically ∆N eff ≲ 10 −4 [92,93], well below the expected precision at the 2σ level of CMB-S4 (∆N eff < 0.06) [79], CMB-HD (∆N eff < 0.027) [94], or COrE and Euclid (∆N eff < 0.013) [83,95].

Summary
One of the largest theoretical uncertainties in mapping the predictions of any inflationary model to measurements of the primordial power spectra is the precise characterization of reheating.Under the assumption of perturbative inflaton decays, the parametrization of the duration of reheating, and the amount of dilution of the energy density, encoded in the reheating parameter R rad , can be simply determined in terms of a single effective inflaton-matter coupling, often parametrized in terms of the reheating temperature T reh [31][32][33][34][35][36][37][38][39][40].The result is a straightforward connection between the main CMB observables n s and r, and the microphysical inflaton couplings to other fields (c.f.Eq. (2.16)).However, when particle production is enhanced (or suppressed) by collective phenomena, the perturbative approximation fails, and a simple analytical description of reheating dynamics is not available.This in turn prevents a simple relation between observables and field couplings.
In this paper we have determined the impact of resonant dark matter production during reheating on the evaluation of CMB observables.Since generality is broken in the presence of nonlinearity, we have limited our numerical study to the Starobinsky model of inflation, a scalar dark sector, and perturbative inflaton decays into Standard Model particles.Specifically, we have considered the four-legged interaction L ⊃ −σϕ 2 χ 2 /2 between the inflaton ϕ and the DM χ, whose strength can be characterized by the dimensionless ratio σ/λ, with λ = m 2 ϕ /2M 2 P .For σ/λ ≳ 5 × 10 3 [46], parametric resonance drives an exponential growth of the DM mode functions, into the backreaction regime.These nonlinear dynamics not only result in a reduction of the comoving energy density of the inflaton by an O(10) factor, but also in the conversion of the remaining coherent, homogenous inflaton field into free inflaton quanta.The efficiency of the decay of the inflaton into SM particles is in consequence reduced, although the transition time from reheating to radiation domination is also reduced.
In order to map the aforementioned shift on the duration of reheating to CMB observables, we have accurately tracked the energy and number densities of the DM field and the inflaton throughout the backreaction regime, with the help of CosmoLattice [64,65].With them, we computed the perturbative SM particle production rate, following the evolution of the energy density of the primordial plasma before, during and after the transient period of preheating.We have found that reheating temperatures are impervious to dark sector preheating, while the duration of reheating, and the number of e-folds N * , can be simply parametrized in terms of the loss in comoving energy density of the inflaton during backreaction, denoted by ∆.This has allowed us to find a functional form for the correction to the CMB observables, n s (∆), r(∆).Although we find that the resulting increase in n s and the decrease in r are too small to be comparable to the sensitivity of next-generation CMB experiments, they modify the favored range of inflaton-SM couplings by up to an order of magnitude.In terms of the reheating temperature, T reh < 203 GeV are excluded at the 95% CL in the absence of preheating, while in the presence of strong preheating T reh can be as low as 12 GeV.
Additionally, we have revisited in this work the calculation of the primordial DM relic abundance produced from the parametric excitation of the dark scalar field, previously explored in [46,47].We have improved the calculation of the DM closure fraction by including the shift in the duration of reheating due to the depletion of the comoving inflaton energy density due to DM preheating.We find that the shortened duration of reheating enhances the relic abundance as Ω χ → Ω χ /∆.The allowed range of reheating temperatures for a given inflaton-DM coupling and DM mass is therefore reduced by the same factor.Table 2 summarizes the constraints on T reh and m χ from structure formation and CMB considerations.Our first conclusion is that only light DM masses are allowed, ≲ 40 GeV at the 1σ level.Our second and main conclusion is that, at the CMB 68% CL, only a narrow range of reheating temperatures is allowed.In particular, for strong preheating with σ/λ = 10 6 , the allowed range for T reh spans only one order of magnitude.We emphasize that these results are applicable only  for resonant DM production in Starobinsky inflation.Nevertheless, as is discussed in [47], we expect our conclusions to be qualitatively unchanged for plateau-like inflationary potentials.Moreover, for alternative realizations of inflation and dark sector preheating, we have provided here a roadmap to accurately constrain the inflaton couplings.Such constraints on the microphysics of inflation are crucial for narrowing down the UV completions of the Standard Model within which the inflationary paradigm must be embedded.

Figure 1 .
Figure 1.Numerically computed curvature (left) and tensor (right) power spectra for the Starobinsky model (2.2), assuming that the Planck pivot scale k P = 0.05 Mpc −1 leaves the horizon 55 e-folds before the end of inflation.The inflaton mass is chosen to match the measured value of A S .The vertical dashed lines show for comparison the location of the Planck and WMAP (k W = 0.002 Mpc −1 ) pivot scales.

Figure 2 .
Figure 2. Scalar tilt n s (left) and tensor-to-scalar ratio r (right), as functions of the number of e-folds after horizon crossing, N * , for the Starobinsky model (2.2).The continuous black lines corresponds to the numerical evaluation of Eqs.(2.23) and (2.27).The gray dotted lines are the slow-roll approximations (2.11) and (2.12), evaluated at horizon crossing.The red dashed curves correspond to the slow-roll approximation, with the substitution N * → N * + 2.5.

Figure 3 .
Figure 3. Floquet charts for the Mathieu equation (3.18) showing the flow lines for selected comoving momenta for the Starobinsky potential (2.2) (left), together with the corresponding evolution of the energy densities of the inflaton and DM during reheating (right), in the absence of backreaction (top) and in the presence of backreaction and fragmentation (bottom).

Figure 5 .Figure 6 .
Figure 5. Top, left: Evolution of the total inflaton energy density (blue), the energy density of the inhomogeneous inflaton component (red), the energy density of the homogeneous inflaton component (green), and the DM energy density (yellow), for σ/λ = 10 4 .Right: Comoving number densities for the inflaton quanta (blue) and the DM quanta (yellow).Bottom, left: the oscillation-averaged total equation-of-state parameter.

Figure 7 .
Figure 7. Inflaton and radiation energy densities after the end of inflation, in the presence of dark sector preheating (continuous) and in its absence (dashed).For the continuous curves, the full inflaton energy density is shown as the blue line, while the contribution of non-relativistic modes is shown in purple.Their difference is highlighted in the inset.The horizontal dotted line corresponds to the energy density at the end of reheating, ρ reh = ρ ϕ (a reh ) = ρ R (a reh ).

Figure 8 .
Figure 8. Number of e-folds at the Planck pivot scale (left) and reheating temperature (right), as a functions of the inflaton-matter Yukawa coupling y, for a selection of inflaton-DM effective couplings σ/λ.

Figure 9 .
Figure 9. Parameter space for dark matter production from strong preheating, for a selection of effective couplings σ/λ.Along each line the measured dark matter relic abundance is Ω χ h 2 = 0.12.For the dashed curves we ignore the change in the reheating duration due to DM backreaction.Along the continuous curves we use the full approximation (4.17) with ∆ ̸ = 1.The shaded region is excluded by the Lyman-α measurement of the matter power spectrum.

Figure 10 .
Figure 10.Scalar tilt (left) and tensor-to-scalar ratio (right), as a functions of the inflaton-matter Yukawa coupling y, for a selection of inflaton-DM effective couplings σ/λ.The left panel shows the 68% and 95% CL regions from Planck in combination with BK18+BAO[5].In the right panel these contours are not shown, since all curves lie inside the 1σ region.

Table 1 .
Lower bounds on the inflaton-SM Yukawa coupling y and the reheating temperature T reh , as functions of the effective coupling σ/λ, at the 68% and 95% CL.These constraints are particular to the Starobinsky model (2.2), derived from the marginalized Planck+BK18+BAO bounds on n s at r 0.002 ≃ 0.004.

Table 2 .
Summary table for the upper and lower bounds on the reheating temperature and dark matter mass for a selection of effective inflaton-DM couplings σ/λ, from structure formation constraints (upper bound) and the 68% and 95% CL Planck+BK18 bounds on n s for Starobinsky inflation (lower bounds).