Energy distribution and equation of state of the early Universe: matching the end of inflation and the onset of radiation domination

We study the energy distribution and equation of state of the universe between the end of inflation and the onset of radiation domination (RD), considering observationally consistent single-field inflationary scenarios, with a potential 'flattening' at large field values, and a monomial shape $V(\phi) \propto |\phi|^p$ around the origin. We include a quadratic interaction $g^2\phi^2X^2$ between the inflaton $\phi$ and a light scalar 'daughter' field $X$, as a proxy for (p)reheating. We capture the non-perturbative and non-linear nature of the system dynamics with lattice simulations, obtaining that: $i)$ the final energy transferred to $X$ depends only on $p$, not on $g^2$; $ii)$ the final transfer of energy is always negligible for $2 \leq p<4$, and of order $\sim 50\%$ for $p \geq 4$; $iii)$ the system goes at late times to matter-domination for $p = 2$, and always to RD for $p>2$. In the latter case we calculate the number of e-folds until RD, significantly reducing the uncertainty in the inflationary observables $n_s$ and $r$.

Introduction -. Cosmological observations strongly support the idea of an inflationary period in the early universe [1][2][3][4]. Inflation must be followed by a (p)reheating stage, where most of the energy in the universe is transferred into light particle species, with only one strong requisite: the universe must arrive at a radiation dominated (RD) thermal state before the start of Big Bang Nucleosynthesis (BBN), at temperatures T BBN ∼ 1MeV. The state of the Universe at BBN is based on the Standard Model (SM) particle content, which is fairly known. However, the way the Universe arrives at this state from the previous inflation stage is largely unknown, and depends strongly on the underlying particle physics model.
Measurements of the cosmic microwave background (CMB) provide an upper bound on the inflationary Hubble rate, H inf 6.6 × 10 13 GeV [5,6], corresponding to energy scales just below ∼ 10 16 GeV. The energy gap between the end of inflation and the onset of BBN may therefore span up to ∼ 19 orders of magnitude. Characterizing this primordial dark age period is important, as it represents a natural 'cosmological window' to probe beyond the SM (BSM) physics, potentially displaying a very rich phenomenology, see [7][8][9][10] for reviews and references therein. Moreover, the equation of state (EoS) during this period is required for making accurate predictions of inflationary CMB observables [11][12][13].
In the context of slow-roll single-field inflation, a preheating phase emerges when the inflaton φ, the field responsible for inflation, starts oscillating around the minimum of its potential. In this work we consider a broad class of observationally viable scenarios inspired by αattractors [14], with 'flattening' of the inflaton potential at large field values, and monomial behavior V (φ) ∝ |φ| p around the origin, with p ≥ 2 including fractional values. The inflaton is directly coupled to a '(p)reheating sector' represented by a light scalar field X, which will be called the daughter field from now on. We consider a quadratic interaction g 2 φ 2 X 2 , as it does not require the introduction of new mass scales, and serves as a proxy for the leading term in gauge interactions [15]. Under these considerations, the universe goes first through a stage of preheating, in which the initially homogeneous inflaton condensate fragments via non-perturbative particle production effects, see [16][17][18][19][20][21][22] for the pioneering studies and e.g. [23][24][25][26] for recent numerical works. Preheating can happen through two separate phenomena: 1) broad parametric resonance of the daughter field, in which the inflaton transfers to the former a large amount of its energy exponentially fast, and 2) self-resonance of the inflaton, in which the inflaton amplifies its own fluctuations. In both cases, a departure from the (initially homogeneous) inflaton oscillation-averaged EoS is ensued, affecting the following expansion history of the universe.
To investigate (p)reheating we use very long classical lattice simulations, considering a large range of inflatondaughter couplings. We are interested in the number of e-folds ∆N RD from the end of inflation till the onset of RD. Previous works [23,25] have obtained this number in the absence of inflaton-daughter interactions, ∆N RD | g 2 =0 . However, the time scale of the daughter field excitation through broad resonance is faster than via inflaton self-resonance, so ∆N RD | g 2 =0 represents only an upper bound to this quantity. Furthermore, recent criticism to (gravitational) reheating in the absence of inflaton couplings to other species [27,28], reinforces the idea that the universe is most naturally reheated if inflatondaughter interactions are present.
In this Letter we consider a large range of inflatondaughter interactions, exploring both a large coupling regime (which leads to broad resonance of the daughter field), and a small coupling regime (which recovers the coupling-less results from [23,25]). We characterize in detail the energy distribution, EoS, and ∆N RD as a function of p and g 2 , and use this information to reduce drastically the uncertainty in the prediction of the inflationary scalar tilt n s and tensor-to-scalar ratio r. Our analysis goes beyond ad-hoc analytical parametrizations of the post-inflationary EoS [29], and beyond previous numerical works [24,30,31], which only considered the initial preheating and early non-linear stage for specific choices of p. Here we simulate, for the first time, the very long-term evolution of the system, for arbitrary values of p ≥ 2.
Parametric resonance and self-resonance -. Consider a scenario with an inflaton φ and daughter field X, where M and Λ are mass scales, and g is a dimensionless coupling. The first term is the inflaton potential, responsible for slow-roll inflation. The interaction term allows to transfer energy between φ and X.
The inflaton potential features a plateau at φ M and an inflection point at entailing that φ always oscillates in the positive-curvature region of the potential, which can be approximated around the origin by the power-law V inf (φ) µ 4−p |φ| p /p, with µ 4−p ≡ Λ 4 /M p for p = 4, and µ 4−p ≡ λ = 1 for p = 4. After inflation, the inflaton oscillates initially as a homogeneous condensate, with decaying amplitude φ(t) ∝ a(t) −6/(p+2) , and timedependent oscillation frequency . This leads to an EoS [32] where p φ osc and ρ φ osc denote the oscillation-averaged pressure and energy densities of the inflaton. For M m pl , two preheating effects emerge due to the initial homogeneous oscillations: parametric resonance of the daughter field and self-resonance of the inflaton. Choosing a * = 1 and re-defining the variables as dz ≡ p+2 (X/φ * ), the linearized mode equations of the daughter and inflaton fields, during the early oscillatory phase, correspond to oscillator-like equations with time-dependent mass terms is an effective resonance parameter, decreasing in time for p < 4, remaining constant for p = 4, and growing for p > 4. This leads to (p)reheating dynamics described by the interplay of two phenomena: broad resonance of the daughter field, and self-resonance of the inflaton field. Fluctuations of both fields evolve as |δχ k | 2 ∝ e 2µ k z and |δϕ k | 2 ∝ e 2ν k z , where µ k ≡ µ k (κ, q res ; p) and ν k ≡ ν k (κ; p) are the respective Floquet indices. These functions are positive for some bands of momenta, leading to an exponential growth of the field modes.
If parametric resonance is broad (q res 1), the range of amplified δχ k is much wider than the one for δφ k . Thus, if both effects are present, the excitation of δχ k is the dominant one. The daughter field is also excited if the resonance is narrow (q res 1), but this effect is negligible compared to broad resonance, and in any case, due to natural limitations of the lattice, it cannot be captured in our simulations. On the other hand, the momenta excited during broad resonance scale (modulo scale factor powers) as p br ∼ q 1/4 * ω * 10 13 GeV, which justifies neglecting a mass term of the daughter field in (1) for m X p br . Results -. We present now our numerical results, obtained from classical lattice simulations of the EOM f −a −2 ∇ 2 f +3Hḟ = −∂ f V for f = {φ, X} and the Friedmann equation, for different choices of p and q * . We have performed simulations in 2+1 and 3+1 dimensions, and checked that they are almost identical, see supplemental material for a direct comparison. However, simulations in 2+1 dimensions have the advantage of allowing us to investigate a much larger region of parameter space. We have used a number of lattice sites per dimension ranging from N = 2 7 to N = 2 10 , and explored different infrared and ultraviolet lattice cut-off's, ensuring a range of momenta encompassing well the scales excited by the different resonances. We have simulated the cases M = 4m pl − 10m pl , which guarantee that the inflaton oscillations occur in the positive-curvature region.
Energy distribution and equation of state. The different energy density components ρ i characterize the evolution of the system. Different contributions include kinetic ρ k,f and gradient ρ g,f energy components, the inflaton potential ρ pot = V (φ) [first term in (1)] and the interaction term ρ int = 1 2 g 2 φ 2 X 2 . As expected from previous studies, the system virializes very quickly [23][24][25]33], with the fields obeying a relation of the type ḟ 2 = |∇f | 2 + f (∂V /∂f ) , where brackets indicate oscillation and spatial averaging. Introducing energy density ratios as ε i ≡ ρ i / j ρ j (so that j ε j = 1 by construction), the virial relations imply The instantaneous EoS w ≡ p/ρ is sourced by the different energy contributions as This means that whenever ε pot , ε int become negligible (as it happens e.g. at later times for p > 2), then ε k,ϕ + ε k,χ 1/2, which leads to a RD universe with w = 1/3. Furthermore, taking averages on both sides of Eq. (6) leads to the effective EoS during the first inflaton oscillations: initially it holds that ε k,ϕ + ε pot 1, so Eq. (4) implies ε k,ϕ p/(p + 2) and ε pot 2/(p + 2), and from there we recover w hom in Eq. (2).  After the initial homogeneous phase, {ε a } and w evolve very differently depending on the choice of p and q * , see panels in Fig. 1. We discuss now their evolution, in particular their asymptotic late time behaviour: • p = 2 : There is no self-resonance of the inflaton field, but the daughter field energy grows exponentially fast via broad parametric resonance, as long as q res > 1 (top-left panel Fig. 1). Thanks to the inflaton-daughter interaction, a growth of the inflaton gradient energy is also induced. However, as q res decreases in time, parametric resonance eventually becomes narrow. For q * 6 · 10 3 , broad resonance persists long enough for backreaction from the daughter field to break the homogeneous inflaton condensate. In such a case, the effective EoS jumps from w hom 0 to a positive value w max < 1/3 (closer to 1/3 the larger q * ). Gradient energies redshift as ∼ a −4 , whereas the inflaton potential/kinetic energies redshift as ∼ a −3 . Therefore, once q res < 1, the daughter energy fractions become gradually negligible, independently of the daughter-inflaton coupling strength. Similarly, the equation of state asymptotically tends to the homogeneous value w → w hom = 0. We observe this (otherwise expected) result for the first time with very long lattice simulations, as shorter simulations in previous works were only able to observe a transitory stabilization of the EoS around w 0.2 [30,31].
• 2 < p < 4 : The inflaton can now develop fluctuations via self-resonance, but if q res 1, the daughter field energies grow much faster via broad parametric resonance. If q * is large enough, backreaction effects from the daughter field break the initially homogeneous inflaton configuration, making the EoS jump from w = w hom to w = w max < 1/3. In fact, a transitory regime of equipartition can be observed for very large values of q * , equally distributing the energy between the two fields (bottom-left panel of Fig. 1). In any case, the resonance eventually becomes narrow when q res = 1, and the daughter field energy fractions become gradually negligible. However, the inflaton self-resonance is still present, which triggers a slow cascade of the inflaton spectra towards the ultraviolet, as well as a growth of its gradient energy at the expense of its potential. This phenomenon, originally reported in [23,25] in the absence of inflaton-daughter interactions, is observed now remarkably even after the inflaton fragments due to the parametric resonance of the daughter field. As a consequence, the EoS always goes to w 1/3 at sufficiently late times. For 2 < p 3 and certain values of q * , the self-resonance is so weak that a temporary decrease of w towards w hom is observed after the end of broad resonance (see top-middle panel), though the system always goes to w 1/3 at later times.
• p ≥ 4 : The resonance parameter q res remains constant for p = 4, or grows in time for p > 4. In the latter case, the system always ends up in broad resonance, even if q * < 1. As inflaton self-resonance effects are also present, the system never ceases to exchange energy between the two fields at late times. Due to this, it achieves an equilibrium state, in which the energy is evenly distributed: 50% of the energy is stored in the daughter field, and the other 50% in the inflaton. The equation of state also goes to w → 1/3 at late times.
Inflationary observables and discussion -. To compute the inflationary scalar tilt n s and tensor-toscalar ratio r, we need to determine the number of efolds N CMB before the end of inflation, when the pivot scale k CMB = 0.05Mpc −1 exited the horizon. For this, we need to know the exact evolution of the universe after inflation. In particular, we need [34,35] with V CMB denoting the potential energy when k CMB leaves the horizon, ρ br and ρ RD the energy densities when backreaction becomes noticeable and when the universe becomes RD, respectively, ∆N br the e-folds between the end of inflation and backreaction (see Fig. 2), andw the mean EoS between backreaction and the onset of RD. For p 3, the transition from backreaction to RD is actually almost instantaneous, independently of g 2 (see Fig. 1), making the last term in Eq. (7) negligible. In Fig. 2 we show ∆N br for different choices of p and g, extracted from simulations. The inflaton fragments due to self-resonance for p > 2 even if g = 0, which provides the value for ∆N br found in [23,25]. However, the presence of an interaction reduces significantly this quantity, as long as backreaction from the daughter field fragments the inflaton condensate. This requires e.g. p 3.4, 2.6 for g = 10 −5 , 10 −4 respectively.
For p = 2, the system never achieves a RD state in our set-up, so we cannot determine N CMB . However, according to our results, the difference in N CMB compared to a case in which the inflaton remains homogeneous until RD is δN CMB 1 for q * < 10 6 . For example, we get δN CMB ≈ 0.4 for the case depicted in the top-left panel of Fig. 1. For p > 2 we can compute N CMB exactly, provided we note that V CMB depends also on N CMB , making (7) a non-linear equation. For p = 3 − 6, we find the narrow range N CMB 56.1 − 56.9 for M = 4m pl and N CMB 56.7 − 57.6 for M = 10m pl . Using this, we obtain very precise values for the inflationary observables: n s 0.9643 − 0.9647 and r 0.01 − 0.009 for M = 4m pl , and n s 0.9633−0.9622 and r 0.047−0.05 for M = 10m pl , see Fig. 3. This represents a drastic reduction in the uncertainty of these quantities, compared to the traditional bounds obtained from N CMB = 50−60.
Discussion. We have characterized in detail the evolution of the energy distribution and effective EoS of the universe from the end of inflation till the onset of RD, considering an inflaton with monomial potential during the (p)reheating stage, and a quadratic coupling to a light daughter field. Remarkable facts emerge: i) Broad parametric resonance dominates over inflaton self-resonance, and backreaction from the daughter field is responsible for breaking the initial homogeneity of the inflaton. However, broad resonance eventually ends for 2 ≤ p < 4. For p = 2 the system goes back to a higher degree of homogeneity, while the EoS approaches gradually the homogeneous value (w = 0). For p > 2, inflaton fluctuations are also created via self-resonance, remarkably even after the breaking of the inflaton homogeneous condensate. Due to this, the system always goes eventually to RD in this case, either in the presence or absence of interactions with a daughter field species.
ii) The final amount of energy transferred into the daughter field is essentially independent of the coupling strength between the two fields, and depends only on the power law exponent p: it becomes (eventually) negligible for 2 ≤ p < 4, and of order ∼ 50% for p ≥ 4. Therefore, in order to achieve a complete decay of the inflaton in these scenarios, some new ingredient is needed.
iii) Viable models of inflation with p > 2 allow for a precise calculation of N CMB , with accuracy δN CMB 1, leading to very precise predictions for n s and r. This highlights the relevance of characterizing the postinflationary stage in detail.
To conclude, we mention some of the limitations of our analysis, which can provide interesting avenues for future studies. For example, our study could be generalized to e.g. trilinear interactions or higher order operators [36], or to an initial excitation via tachyonic preheating [37][38][39][40]. Oscillons can also form whenever the inflaton oscillates around flatter-than-quadratic regions of the potential, via self-resonance effects [41] or tachyonic oscillations [42]. This would push the EoS towards w 0 [23,25,43] during their lifetime. While there are various effects that can change the early stages of preheating, we expect that the late-time energy distribution and EoS will depend mainly on the inflaton potential around its minimum, and on the type of inflaton-daughter coupling. It would be also interesting to generalize our analysis to multi-field inflation scenarios [44][45][46][47][48]. Finally, if the inflaton is coupled to several light scalar fields, preliminary lattice simulations (for quadratic couplings) show that the energy transferred to the preheat sector can be enhanced up to N f /(N f + 1)%, with N f being the number of different light daughter fields. We plan to explore some of these topics in the future.