Ricci reheating reloaded

A Hubble-induced phase transition is a natural spontaneous symmetry breaking mechanism allowing for explosive particle production in non-oscillatory models of inflation involving non-minimally coupled spectator fields. In this work, we perform a comprehensive characterisation of this type of transitions as a tachyonic Ricci-heating mechanism, significantly extending previous results in the literature. By performing 𝒪 (100) 3+1-dimensional classical lattice simulations, we explore the parameter space of two exemplary scenarios, numerically determining the main timescales in the process. Based on these results, we formulate a set of parametric equations that offer a practical approach for determining the efficiency of the heating process, the temperature at the onset of radiation domination, and the minimum number of e-folds of inflation needed to resolve the flatness and horizon problems in specific quintessential inflation scenarios. These parametric equations eliminate the need for additional lattice simulations, providing a convenient and efficient method for evaluating these key quantities.


Introduction
The success of inflation is based on its ability to solve the main shortcomings of the hot Big Bang model while providing a mechanism able to generate the primordial density perturbations seeding structure formation.Although numerous theoretical models of inflation have been already excluded by observations, current data sets are still insufficient to identify a particular scenario within the whole inflationary paradigm.Furthermore, there is still a rather limited understanding of how a specific inflationary model should be embedded into a more fundamental particle physics theory.The situation is expected to improve, however, in the next decade, when missions such as the upgraded Simons Observatory [1], the LiteBird satellite [2,3] or the ground-based CMB Stage 4 program [4] will significantly reduce the uncertainties in the tensor-to-scalar ratio [5].In this context, understanding the details of the heating stage following the end of inflation [6][7][8] becomes of the uttermost importance for providing accurate theoretical predictions to be compared with observations.
A minimal Hubble-induced heating mechanism based on the time-dependence of the Ricci scalar was recently put forward in the context on non-oscillatory quintessential inflation scenarios, where the absence of a potential minimum gives rise, rather generically, to a kineticdominated post-inflationary era [9,10] (cf.also Refs.[11][12][13] and Ref. [14] for a comprehensive review of quintessential inflation).In this setting, any self-interacting spectator field nonminimally coupled to gravity may potentially undergo a second-order phase transition when the Universe evolves from inflation to kination and the Ricci scalar turns negative.Generally speaking, the transition to the new minima of the potential appearing at large field values involves a tachyonic instability that comes to an end only when the effective mass of the spectator field becomes again globally positive.This symmetry restoration can either take place at the eventual onset of radiation domination or be effectively induced by the non-linear dynamics following particle production.In particular, if the energy lost by the spectator in each oscillation is smaller than the one associated to the motion of the new Hubble-induced minima as the Universe expands, the field distribution will be able to cross back the origin of the potential [9], giving rise to a bimodal distribution that disappears after a few oscillations [10].
The main purpose of this work is to perform a full characterisation of Hubble-induced phase transitions as a tachyonic-heating mechanism, significantly expanding the results in Refs.[10,15].By performing a large number ≃ O(100) of 3+1-dimensional classical lattice simulations we will explore a wide portion of the parameter space of the theory.In particular, the outputs of our simulations will allow us to numerically characterise the main timescales in the process, namely the initial tachyonic instability, the onset of backreaction effects and the virialisation phase leading to a radiation-like equation-of-state for the spectator field.These numerical results constitute the starting point for obtaining a set of ready-to-use parametric formulas characterising the total energy-density at each timescale, the efficiency of the heating process and the radiation temperature at the onset of radiation domination for all the model parameters in the scan, in the spirit of Ref. [16].This line of work comes with several benefits.First, the non-linear phenomenology of Hubble-induced phase transitions is investigated comprehensively, without resorting on ad hoc homogeneity assumptions and with a strong emphasis on the interplay between the different model parameters.Second, the simple relations given by the fitting formulas allow us to accurately determine the minimal number of e-folds of inflation needed to solve the flatness and horizon problems in specific quintessential inflation scenarios, without resorting to additional time-consuming lattice simulations.Third, the high efficiency of the tachyonic process enable us to directly apply our results to a plethora of variations of the model including, for instance, quartic couplings among species leading potentially to energy transfer effects.
This work is organised as follows.In Section 2, we review the basic working assumptions of Hubble-induced phase transitions, focusing for concreteness on a Z 2 -symmetric model.After deriving analytical solutions for the tachyonic amplification of fluctuations in the linear regime, we proceed to study the non-linear dynamics of the system using 3 + 1 classical lattice simulations.The parametric formulas for the heating efficiency and the radiation temperature following from this procedure are presented in Section 3, where we describe also their impact on the main inflationary observables and the potential effects of non-gravitational interactions with other beyond the Standard Model sectors.Finally, our conclusions are presented in Section 4.

Spectator field dynamics
The simplest realisation of a Hubble-induced phase transition is described by an action containing an energetically-dominant quintessential scalar field ϕ [14,17,18] and a subdominant Z 2 -symmetric spectator field χ non-minimally coupled to gravity, with R the Ricci scalar and M P = (8πG) 1/2 = 2.4 × 10 18 GeV the reduced Planck mass.The precise form of the inflaton Lagrangian density L ϕ in this expression will not play a relevant role in the following discussion, being the associated scalar field only responsible for driving the background evolution of the Universe at early times and, most importantly, for inducing a transition from inflation to a kinetic-dominated epoch.Regarding the spectator field counterpart, it includes a dimensionless coupling constant ξ, only restricted a priori by the self-consistency of the procedure.In particular, the field-dependent contribution to the overall Ricci scalar prefactor is required to be subdominant as compared to the usual Planck mass counterpart, namely ξχ 2 ≪ M 2 P .In this limit, the graviton propagator is well-defined for all field values of interest, while coinciding approximately with the graviton propagator in General Relativity.Moreover, in order to avoid the formation of large isocurvature perturbations, the strength of the non-minimal coupling is restricted to exceed 1/12, guaranteeing the spectator field to be sufficiently massive as not to undergo large displacements during inflation [9,10,19].As it will become clear in Sections 2.1 and 2.2, a more restrictive lower bound comes from minimizing the effects caused by reentering modes close to the Hubble horizon, which enforces ξ ≳ 15.In this parameter range, inflationary isocurvature perturbations are safely under control.Finally, note that, although the phenomenology of Hubble-induced phase transitions can be easily extended to arbitrary monomial potentials [14,15], we have intentionally restricted ourselves to a marginal quartic self-interaction with dimensionless coupling constant λ.This avoids the introduction of additional energy scales beyond the constant Planck mass and the dynamical scalar curvature entering the Einstein-Hilbert term.
In a Friedmann-Lemaître-Robertson-Walker background g µν = diag(−1, a 2 (t)δ ij ) with scale factor a(t), the Klein-Gordon evolution equation following from the variation of Eq. (2.1) with respect to χ takes the form where dots denote derivatives with respect to the coordinate t.The energetic subdominance of the spectator field allows us to interpret the Ricci scalar in this expression as an effective timedependent mass for χ, whose dynamics in a Friedmann-Lemaître-Robertson-Walker metric is completely determined by the evolution of the inflaton field, namely with H = ȧ/a the Hubble rate and w ϕ the effective inflaton equation-of-state parameter.In these expressions, the Ricci scalar constitutes a cosmic timekeeping device, changing sign precisely at the transition between inflation (w ϕ = −1) and kination (w ϕ = 1) and triggering with it the spontaneous breaking of the Z 2 symmetry in the χ sector of the theory.Given the time-dependence of the resulting symmetry-breaking potential, the true vacua appearing at large field values migrate monotonically towards zero as the Hubble function decreases with time.
In a generic quintessential-inflation scenario, the transition from inflation to kination takes place on a timescale dictated by the detailed shape of the inflationary potential at the end of inflation.Since we are mainly interested in the production of particles upon spontaneous symmetry breaking, we will assume the crossover stage between these two cosmological epochs to be essentially instantaneous, neglecting with it any particle creation effects prior to the onset of kinetic domination.This assumption allows us to parametrise the evolution of the scale factor and the Hubble rate as with a kin = a(t kin ) and H kin = H(t kin ) the associated values at the beginning of kination.In order to remove the Hubble-friction term in Eq. (2.2), it is also convenient to transform the field and spacetime coordinates to their conformal version where τ is the usual conformal time variable (dτ = dt/a) and χ kin = √ 6ξH kin is a typical field value associated to the initial curvature of the effective potential around the origin.Taking into account all these changes, the final version of the Klein-Gordon equation for the scalar field Y takes the form with an effective time-dependent mass term and From Eq. (2.1), we can additionally compute the spectator energy-momentum tensor and, in particular, its energy density and pressure [10] ) (2.10) Note that, due to the non-minimal coupling to gravity, these quantities turn out to depend non-quadratically on the spectator field value and its velocity, being therefore not guaranteed to be positive definite at all times [20][21][22][23].Although this well-known fact leads to the violation of the weak energy condition within the spectator field sector, the total energy density of the Universe, i.e. that including the dominant inflaton contribution, is still positive definite at all times.

Linear regime
At this stage, one could be tempted to treat the spectator field χ as an homogeneous component akin to the runaway inflaton field ϕ, reducing the analysis of Hubble-induced phase transitions to the mere integration of the homogeneous equations of motion, as done for instance in Refs.[9,11,13].Despite its simplicity, this naive approach is based on the false premise that the only relevant mode for the spectator field's evolution is the zero mode.On the contrary, as argued in Refs.[15,19,24] and shown explicitly in Ref. [10], the existence of spatial gradients plays a crucial role in the formation of topological defects, in the overall energy budget and in determining the effective equation of state of the spectator field (cf.also Refs.[25,26] and [27,28] for previous analyses supporting this view).Moreover, the homogeneous-field approximation does not respect the fundamental Z 2 symmetry of the model under consideration, requiring ad hoc initial conditions that cannot be sourced by stochastic quantum fluctuations after inflation comes to an end.
Having in mind the above arguments, let us proceed to the analysis of the spectator field fluctuations.To this end, we notice that immediately after symmetry breaking the nonlinear term in the evolution equation (2.6) plays a rather subdominant role, allowing us to approximate this expression by that for a free scalar field with time dependent mass M (z), making therefore the problem exactly solvable.As described in detail in Refs.[10,14,15], the associated spectator field quantisation proceeds according to the standard procedure in non-trivial backgrounds [29,30].In particular, each Fourier mode κ = k/(a kin χ kin ) of the scalar field Y is promoted to a quantum operator that satisfies the equal-time commutation relations [ Ŷκ (z), Πκ ′ (z)] = iδ(κ + κ ′ ), where Πκ = d Ŷκ /dz, âκ and â † κ are the standard annihilation and creation operators and f κ is a set of isotropic mode functions fulfilling the Klein-Gordon equation To solve this initial-value problem, we will assume the ν parameter within the effective square mass M 2 (z) to be large enough as to sufficiently separate the scales corresponding to the Hubble horizon and the typical momenta amplified by the tachyonic instability, thus reducing potential horizon-reentry effects [15].Given this requirement and assuming vacuum initial conditions , the Klein-Gordon equation (2.12) admits a solution [15] with J ν and Y ν the Bessel functions of the first kind and the quantities α κ , β κ , δ defined as This solution consists of a rapidly growing mode (J ) and a decaying mode (Y).By considering only the former, one can easily obtain the following approximated expression for the amplitude of the spectator field mode functions in the large-ν limit [15], which, for highly-symmetric spherical configurations, leads to a two-point correlation function with root mean-square perturbation (rms) and shape function (2.20) Here erfi stands for the imaginary error function and κ * (z) The analysis of the system in its linear regime gives us the time-evolution of the amplified scalar modes during the initial tachyonic instability.Their exponential growth leads to a rapid enhancement of the mode occupation number and to the classicalisation of the system [10,14,15].We will use this result for two different purposes.Firstly, it allows us to estimate the typical momenta at play, namely the lowest and highest momenta in the tachyonic band (see Eq. (2.13)) and the typical momentum κ * .This information is crucial when we perform numerical lattice simulations of the system in its full non-linear glory.Secondly, we can take advantage of the solution in Eq. (2.19) to characterise the moment at which the nonlinear backreaction effects become so large as to end the linear regime and the tachyonic amplification.We will be able to define a univocal criterion that can be applied to the analytical solution as well as to the numerical results, thus allowing a direct comparison.

Non-linear regime
As quantum fluctuations grow during the tachyonic phase, the previously-neglected quartic self-interaction term in Eq. (2.6) becomes increasingly important, up to the point of completely halting the tachyonic growth of perturbations and marking the beginning of the oscillations of the spectator field.The assumptions at the core of the linear analysis in Section 2.1 cease then to be valid.In order to study this highly non-linear regime, we will resort in what follows to fully-fledged classical lattice simulations in 3 + 1 dimensions.To this end, we will make use of a modified version of the publicly-available CosmoLattice code [16,31], designed to evolve interacting scalar fields in an expanding background.In particular, we modified the default lattice implementation of the Klein-Gordon equation in this numerical tool to account for the Hubble-dependent mass term in Eq. (2.2), introducing also related definitions of the energy density and pressure in Eqs.(2.9) and (2.10).
In all our simulations, the lattice parameters are chosen to ensure the stability and reliability of the output.In particular: 1.The comoving grid size L and the number of lattice points per dimension N are selected in such a way that all relevant modes are always well within the associated infrared (IR) and ultraviolet (UV) resolution in momentum space, In particular, these quantities are set to properly cover the tachyonic band found in Eq. (2.13) as well as modes enhanced by subsequent rescattering effects.Therefore, following the results of the linear analysis, we identify the smallest momentum in the tachyonic band with κ IR = H, while the largest amplified momentum is set to be smaller than the lattice's UV momentum, i.e. √ 4ν 2 − 1H ≪ κ UV .The last condition implies a constraint on the minimum number of lattice sites N > 2 √ 4ν 2 − 1/ √ 3 which is always fulfilled in our simulations, as much larger lattices are needed to cover the late-time dynamics of the system.Since the mode excitation depends mainly on the non-minimal coupling parameter ν, we set N = 128 for ν < 10, N = 256 for 10 ≤ ν < 25 and N = 512 for ν ≥ 25.In order to verify that our results do not depend on the lattice size, we perform when possible a few test simulations with increased lattice size.
2. The time-step variable is chosen according to the stability criterion δt/δx ≪ 1/ √ d [31], with d = 3 the number of spatial dimensions and the length of the side of a lattice cell.More specifically, we set δt = 0.1 for ν ≥ 10 and δt = 0.01 for ν < 10.
Given the fixed power-law background expansion in (2.4), a symplectic 4th order Velocity-Verlet evolver guarantees stability and precision of the numerical solutions when the conservation of energy cannot be explicitly checked.In agreement with the overall picture developed in the previous section, the initial conditions for our lattice simulations are set as χ(0) = χ ′ (0) = 0, with fluctuations over this homogeneous background included as Gaussian random fields, as done customarily for systems with short classicalisation times like the one under consideration [10].The resulting evolution is therefore deterministic up to a randomlygenerated initial seed.A robust (and time-consuming) approach would require to perform several realisations of the same simulation, averaging subsequently over the final output to make sure that the random initial conditions do not influence the dynamics.However, because of the computational resources at our disposal, we choose to keep the base seed constant in all our simulations, making them exactly comparable.Since the system looses memory of the initial conditions soon after the development of the tachyonic instability, this choice does not influence the overall macroscopic evolution, as we have explicitly checked by performing a few simulations with different base seeds.

The timescales of heating
In this section, we present the main results of our extensive parameter scanning in the (ν, λ) plane, obtained by performing repeated 3+1 dimensional classical lattice simulations with λ ∈ [10 −6 , 10 −1 ] and ν ∈ [5,35].The domain of ν is set to avoid sizeable horizon-reentry effects (ν ≥ 5) while ensuring the covering of all relevant modes throughout the evolution (ν ≤ 35).As will become clear in what follows, the main timescales in the evolution of the system for ν > 35 can be inferred from simulations with smaller non-minimal couplings.The range of λ is chosen to ensure perturbativity (λ < 1) while partially covering the typical values appearing in Standard Model settings (10 −6 ≲ λ ≲ 10 −1 ), namely those of the Higgs self-coupling at energies around 10 10 − 10 16 GeV [32][33][34].In principle, the requirements of energetic subdominance of the spectator field and of sub-Planckian contribution to the reduced Planck mass set boundaries on the parameter space as well.However, being these boundaries dependent on the scale of kination H kin , we can proceed independently from its exact value and work with dimensionless quantities such as χ(z)/H kin and ρ(z)/H 4 kin , as we did in the previous section through χ kin .The discussion on the correction to the Planck mass and to the energetic subdominance of the spectator field will be made a posteriori.
A typical output of our simulations is presented in Fig. 1, where we display the evolution of the spectator-field energy density, separated in its kinetic (K), gradient (G), potential (V ) and interaction (I) contributions, With these definitions, the interaction term I contains a total-divergence piece that averages to zero once integrated over the full lattice volume with periodic boundary conditions [35,36].Note that both the interaction term I and the total energy ρ χ are displayed as absolute values in Fig. 1, since the non-minimal coupling (tachyonic) contribution is not positive definite.However, the sum of the inflaton and spectator fields energy densities is always positive definite [20][21][22][23].Note also that in Fig. 1 each term has been multiplied by a a 4 factor, in order to highlight the radiation-like evolution of the kinetic and gradient terms at late times.
The presence of the non-minimal coupling to gravity leads to a rapid growth of fluctuations that comes to an end when the non-linearities associated to the quartic self-interaction term become relevant.The backreaction time z br can be identified with the first time at which the second time-derivative of the spectator field vanishes, Y ′′ (z) = 0, which is itself related to the time at which the effective frequency of the field's oscillations vanishes, (2.24) Contrary to other choices in the literature, this definition can be conveniently applied to both analytical and numerical results, enabling their direct comparison. 1 In particular, due to the stochastic nature of the lattice approach, we choose to consider the lattice-averaged amplitude of the fluctuations as a good indicator of the overall dynamical evolution, comparing the  In the upper panel, dashed coloured lines correspond to the analytical curve obtained from Eq. (2.25).In the middle and lower panels, solid lines are derived from the numerical output of lattice simulations while the black dashed lines correspond to the fitting formulas (2.26) and (2.27).
time at which the root-mean-squared ⟨Y (z) ′ ⟩ rms reaches its maximum value to the analytical definition of z br following from with κ * = 2 √ 2ν + 1H the typical momentum scale entering the analytical two-point correlation function in Eq. (2.18).
The behaviour of the backreaction time z br following from the analytical expression (2.25) is displayed in the upper panel of Fig. 2 as a function of the model parameters in the scan.We observe that, due to the fast-growth of tachyonic modes upon symmetry breaking, this quantity is roughly comparable with the duration of the first semi-oscillation, z br ≈ O (10).The explosive enhancement of fluctuations reaches, however, a saturation point at large ν, with different lines approaching an almost constant value in one-to-one correspondence with the choice of λ.Consequently, for non-minimal couplings ν ≳ 20, the difference in backreaction times at constant λ becomes completely negligible.We also recall here that our analytical estimate of the backreaction scale relies on the large-ν limit, a factor that can contribute to the small discrepancies in the interval 5 ≤ ν ≤ 10.
From the data points in Fig. 2 we can additionally derive fitting formulas for the amplitude of the spectator field at backreaction time and the total energy of the spectator field at that moment, which are all simple first-order polynomials in n = − log(λ).These formulas are displayed as black dashed lines in the middle and lower panels of Fig. 2. Note that the mild dependence on the quartic self-coupling allows us to obtain a reliable fit on a logarithmic interval for λ that spans six orders of magnitude.
As mentioned previously, the self-consistency of our procedure requires the explosive tachyonic enhancement of spectator field fluctuations to fulfil the condition on negligible contributions to the reduced Planck mass.The parametric formula (2.26) allows us to estimate the peak amplification of the scalar field and with it the effective change in M 2 P at the level of the model's Lagrangian. Figure 3 shows with solid lines the 10% threshold for different choices of the energy scale at the beginning of kination.Given the self-similar evolution of the lines representing ⟨χ 2 br ⟩ rms in Fig. 2, we can extrapolate the maximal amplitude even for λ < 10 −6 , thanks to the parametric dependence in (2.26).The red rectangle indicates the actual region explored in our scanning, which confirms that our setup is fully consistent with a subdominant scalar field for H kin ≤ 10 12 GeV.
For z > z br , the system enters a highly-oscillatory regime driven initially by the Hubbleinduced contributions (cf.Fig. 1).The impact of these time-dependent terms becomes, however, completely negligible at later times, where the steady flow of energy from the IR to the UV part of the spectrum effectively virializes the system [10], thus approaching an asymptotic energy configuration with ⟨K⟩ ≈ ⟨G⟩ + 2⟨V ⟩.This condition provides a natural lattice definition for the radiation time z rad at which the spectator field starts behaving as a radiation fluid, namely with the subscript t denoting a time-average over several oscillations of the scalar field.Note that this indirect definition of z rad holds as long as the contribution of the non-minimal coupling terms is negligible, coinciding in that limit with the usual condition w χ (z rad ) = 1/3, as can be easily inferred from Eqs. (2.9) and (2.10).Fig. 1 shows that the condition (2.30) is The dashed red rectangle delimits the range explored in our extensive scanning of the parameter space.first fulfilled at z rad ≈ O(100), when the energy stored in the non-minimal interaction is at least one order of magnitude smaller than the kinetic, gradient and potential counterparts.
The specific dependence of the radiation time z rad on the model parameters in the scan is displayed in the top panel of Fig. 4. In deriving these results we have checked that the contribution of the non-minimal coupling becomes subleading as the original symmetry of the theory is restored and the Hubble function decreases with time.This translates into the condition I < 0.1 × K, G, V which is satisfied by all the data points in the figure.For small values of ν, the virialisation process is completed soon after the non-linear regime sets in.At that stage, the spectator field is still undergoing large oscillations, making the precise determination of z rad more difficult and explaining the irregular features at ν < 15.We note also that, due to a stronger rescattering of modes, z rad is generally shorter for higher values of λ.
The above results allow us to derive fitting formulas for the timescale z rad as function of the model parameters.For all values of λ, the behaviour of z rad deduced from Fig. 4   second-order polynomials in n = − log(λ).We can additionally derive a fitting formula for the total energy of the spectator field at z rad , namely This is displayed in the bottom panel of Fig. 4 as black dashed lines.By comparing the bottom panels of Fig. 2 and Fig. 4, we notice that the total energy available at the end of the virialisation process is roughly a constant fraction of the energy produced during the tachyonic phase.

is linear
In particular, the fitting formulas (2.27) and (2.33) tell us that ρ χ rad (λ, ν) ∼ 10 −4 × ρ χ br (λ, ν).At this stage, it is instructive to compare our finding with the results of an homogeneous treatment, where the system dynamics gives rise to an effective equation-of-state parameter w χ = 1/3 soon after the onset of the oscillatory regime [9,10].As apparent in Fig. 1, the presence of spatial inhomogeneities is certainly not a negligible effect, with the energy density into gradients accounting to a sizeable fraction of the total energy density and generically exceeding the potential energy contribution.This is confirmed in Fig. 5, where we display the gradient-to-total energy ratio G/ρ χ computed at the radiation time z rad for all the model parameters in the scan.Since the initial tachyonic production is more effective at sourcing spatial inhomogeneities, this ratio is generically larger for larger values of ν, leading to gradient contributions of up to 40% of the total energy density.The kink feature visible at ν = 25 is a mild numerical effect due to the change of lattice size, which slightly modifies the relative importance of gradients.

Heating efficiency and temperature
After a relativistic equation of state at z rad the spectator field's energy density evolves proportionally to a −4 , with the associated power spectrum approaching slowly a thermal distribution [39][40][41].From there on, the evolution of the Universe's energy budget can be studied analytically via the so-called heating efficiency introduced in Ref. [42].
The most straightforward definition of this quantity is based on the ratio of the energy density of the spectator field at the onset of kination to that of the inflaton component, Θ ht = ρ ψ (a kin )/ρ ϕ (a kin ), which is valid in case of instantaneous production of relativistic daughter fields ψ at the end of inflation [14,42].However, in the present model the tachyonic production of χ particles is completed in about one e-fold after the onset of kination, taking another e-fold for the spectator field to reach a radiation-like equation of state.Therefore, it is more convenient to define the heating efficiency at the onset of the radiation-like behaviour namely with ρ ϕ (a) = 3M 2 P H 2 kin (a/a kin ) −6 the energy density of the inflaton field in the kinetic dominated regime.Note that, since Θ ht | rad = Θ ht | kin (a rad /a kin ) 2 , it is always possible to switch between the two alternative definitions.By inserting the parametric formulas (2.31) and (2.33) in Eq. (3.1) we gain also a parametric expression for Θ ht , showing explicitly its dependence on the model parameters.As displayed in Fig. 6 for the specific case of H kin = 10 11 GeV, the computed values fulfil the big bang nucleosynthesis bound Θ ht ≳ 10 −16 derived in Ref. [14].Higher values of this quantity correspond to longer and stronger tachyonic phases from large-ν and small-λ regimes, as demonstrated by the tight correlation between Θ ht (λ, ν), ρ rad (λ, ν) and ρ br (λ, ν) in Figs. 2, 4 and 6.
The expression for the heating efficiency in Eq. (3.2), together with Eqs.(2.26), (2.27), (2.31) and (2.33), constitutes the central result of this work.The timescales corresponding to backreaction time and to radiation time mark the interval where the complicated non-linear dynamics has to be solved with numerical simulations.The initial tachyonic phase is well captured by the analytical solution to the mode equation in (2.14).The parametric formulas for ρ χ rad and z rad summarise the macroscopic characteristics of the complicated self-interacting system that we have numerically simulated in more than one hundred runs.The quality of the fitting procedure has been checked for all of our fitting formulas via the R 2 method, yielding always R 2 > 0.99.Given the simplicity and reliability of our results, the parametric formulas deliver a ready-to-use description of the main phases of heating, eliminating the need for further simulations over (and beyond) the parameter space we have explored.
Although the spectator scalar field achieves a radiation-like equation of state at z rad , a fully thermalised state is not reached until much later, once the non-linear rescattering processes has redistributed the available energy-density within a large range of modes [10,[39][40][41].Nonetheless, it is convenient to define an instantaneous temperature scale, with g ht * = 106.75 the Standard Model number of relativistic degrees of freedom at energies above O(100) GeV and ρ χ ht = ρ χ (z ht ) = ρ ϕ (z ht ) the total energy-density of the scalar field at the end of the heating phase.This temperature should be understood as a mere estimate of the typical energy of a particle at the onset on radiation domination, coinciding with the most common concept of (re)heating temperature in the fast thermalisation limit.Independently of its thermal character, this scale plays an important role in the production of dark matter relics via UV-freeze-in [43,44] and many mechanisms aiming to explain the observed matterantimatter asymmetry of the Universe [45].Assuming that the overall equation of state of the Universe is w ≈ 1 until the end of heating, we obtain with the heating efficiency Θ ht containing an implicit dependence on the scale of kination, cf.Eq. (3.2).As expected, due to their long-lasting tachyonic instabilities, simulations with smaller values of λ are characterised by a significantly higher heating temperature, up to T ht ∼ 10 12 GeV for H kin = 10 11 GeV.However, even the lowest temperature in our simulations, T ht ∼ 10 5 GeV, falls way above the Planck and big bang nucleosynthesis bounds on the heating temperature, T ht ≥ 5 MeV [46,47].
The parametric formulas for ρ χ rad and z rad allow also for the computation of the minimal number of e-folds of inflation N needed to solve the flatness and horizon problems in scenarios where the Hubble-induced phase transition is the only particle production channel.As computed in Ref. [42], this is given by with the subscript hc referring to quantities evaluated at the horizon crossing of the pivot scale k hc = 0.002 Mpc −1 , kin, rad, mat denoting respectively quantities evaluated at the beginning of the kination, radiation and matter-domination and 0 indicating their present-day values.Numerically this corresponds to [42] where the parametric formulas (2.31) and (3.2) can be directly included.

Impact on inflationary observables
The heating process summarised in the parametric formulas (3.2) and (3.4) plays an important role when computing the total duration of inflation, with the contribution of heating efficiency in Eq. (3.6) accounting for up to ∼ 4 e-folds of expansion history for H kin = 10 11 GeV.Given the capability of future Cosmic Microwave Background experiments to detect a tensor-toscalar ratio at the level of r = O(10 −3 − 10 −4 ) while marginally improving the current constraints on the spectral tilt n s [1,3,4], it might become eventually possible to set bounds on the physics of heating [48,49].Having this in mind, let us make use of the lattice-based relations (3.6) and (3.2) to determine the precise value of the inflationary observables in a specific quintessential-inflation scenario involving only gravitational particle production.Among the many inflationary models in the literature, one can consider, for instance, the simple α-attractor family, encoding, among others, several incarnations of the Higgs inflation paradigm [50,51] and covering a big portion of the (n s , r) space [52] currently favoured by the Planck/BICEP2/Keck collaborations [53] using a single parameter α.The presence of the pole restricts the field's movement in field space, imposing a bound on its value.By transforming to a canonical field, the pole can be effectively shifted to infinity while stretching the potential for the canonical variable near the pole, resulting in the emergence of plateaus [54][55][56].
Although initially formulated for oscillatory models, the α-attractor idea can be easily adapted to the quintessential inflation paradigm by considering non-symmetric potentials [57][58][59].The simplest realisation involves a linear term, adjusted to satisfy the Planck/COBE normalisation constraint, and a small cosmological constant Λ responsible for the late-time acceleration of the Universe, namely V (ϕ) = γϕ + Λ with γ/( √ αM 3 P ) ∼ 2 × 10 −11 in order to satisfy the normalisation of the primordial power spectrum [57].When written in terms of a canonically-normalised field φ = √ 6αM P arctanh ϕ/( √ 6αM P ), this potential becomes effectively stretched, coinciding in the large field limit φ ≫ √ 6αM P with the T-model potential and sharing therefore the very same predictions for the spectral tilt n s , the tensor-to-scalar ratio r and spectral amplitude A s , (3.9) Combining these expressions with Eqs.(3.2) and (3.6), and neglecting the variation of the Hubble rate during inflation, V 0 ≃ 3M 2 P H 2 kin , we can link the inflationary observables to the heating efficiency following from our detailed lattice simulations.The result is displayed in Figure 7, where we have computed r and n s assuming specific values for the Hubble scale H kin .The vertical grey lines correspond to the predictions for the maximum and minimum number of inflationary e-folds (N max = 63.8,N min = 58.0)obtained from the minimum and maximum heating efficiency respectively (cf. Figure 6).The region between the black lines is therefore the full area spanned by the parameter space (λ, ν) of our simulations.

Beyond the single field approximation
The prototypical case of a Z 2 -symmetric scalar field considered in the previous section provides us with an extensive understanding of the tachyonic phase (both from the analytical and numerical point of view) and of the timescales of the heating process.This is to be understood as a basic scenario whose phenomenology can be extended to more elaborated field content.
Let us consider, for instance, a scale-free interaction 1 2 g 2 χ 2 φ 2 that respects the global Z 2 symmetry of the spectator field, with φ some scalar field in a beyond the Standard Model sector and g 2 a dimensionless coupling constant.This term, customary in the literature of heating [38,60], opens up the possibility of depleting the spectator field energy density via perturbative and non-perturbative effects.Given the non-linear nature of the self-scattering interaction and the coupling between fields, let us turn once more to lattice simulations.Due to the fast population of the daughter field's UV modes, a larger lattice is generically needed to cover the evolution of the field's fluctuations, especially in case of large g and prolonged tachyonic phases for χ (i.e.small λ).In order to maintain a good coverage of all momenta within a lattice box of N = 256, we will restrict ourselves to model parameters in the fiducial intervals g ∈ [10 −4 , 10 −3 ], ν ∈ [10,20] and λ ∈ [10 −5 , 10 −1 ], considering only particular benchmark points due to the longer computational time needed for each simulation.
The time-evolution of the different terms contributing to the total energy-density in a typical simulation is displayed in Fig. 8.As noted already in Section 2.3, the overall evolution of the spectator field is characterised by three main phases: the initial tachyonic amplification of the field's modes, the evolution towards virialisation and the thermalisation of the system.Additionally, we can observe one additional timescale related to the daughter field.The equipartition timescale marks the moment when the total energy density is evenly distributed between the two scalar fields.A numerical computation of this timescale is a challenging task, since it can require a very long time to be reached, depending on the production mechanisms and on their efficiency.We simply define it through the condition ρ χ (z eq ) = 1  2 ρ tot (z eq ) for those simulations in which equipartition is reached within their total duration.
The parametric dependence of the timescales z br , z rad and of the energy densities ρ χ br , ρ χ rad following from our simulations is presented in Fig. 9.We observe that, for the values of g considered in our simulations, the end of the tachyonic instability is caused primarily by the self-interactions of the spectator scalar field, thus matching closely the results obtained in the single field case.After the peak production of the spectator field modes, the associated kinetic and gradient energies evolve as radiation, with both fields contributing to the overall virialisation of the system.The timescale z rad is computed numerically by extending the condition in (2.30) to a 2-component interacting fluid, i.e. by including the additional kinetic, gradient and interaction terms.Again, we notice a good matching between z rad for the single-and two-fields scenarios.The interval between z rad and z eq has a variable duration, depending on many factors, including the intensity of the production processes, the appearance of parametric resonances and the magnitude of the coupling constant g.For the cases (g = 10 −3 , λ = 10 −3 ) and (g = 10 −4 , λ = 10 −4 ), the two fields are reaching equipartition of the total energy already during the virialisation phase.This phenomenon induces stronger oscillations on the spectator field which makes the estimation of z rad somewhat less reliable, see the middle panel of Fig. 9. On the other hand, a short equipartition timescale allows us to compute it numerically for the special cases (g = 10 −3 , λ = 10 −3 ) and (g = 10 −4 , λ = 10 −4 ), finding it to be in the range 250 < z eq < 400.After rescattering between fields has redistributed the total energy density, the system evolves towards the end of the heating phase.We notice once more that the results in Fig. 9 are in excellent agreement with the previous results for the a single spectator field.Lastly, Fig. 10 displays the heating efficiency for the two-field setup with g = 10 −3 and g = 10 −4 compared with the parametric formula obtained from Eq. (3.2).
It is interesting to notice that for sufficiently large couplings g ∼ O(0.1) and small selfinteractions λ, the backreaction is expected to be produced by the daughter field itself.In this scenario, the rapid growth of φ fluctuations could dynamically generate an effective mass for the spectator field at large times, even if the quartic self-coupling is completely negligible (λ → 0).Unfortunately, the energy transfer to UV modes in this case is fast and affects a Figure 9: Dependence of the timescales z br (first panels), z rad (third panels) and of the energy densities ρ χ br (second panels) and ρ χ rad (fourth panels) on the model parameters ν and λ.The daughter coupling as been set to g = 10 −3 for all simulations.Dots represent the numerically-computed timescales in the two-fields scenario.Dashed lines refer to the fitting formulas obtained from Eq. (2.27) and (2.33) while solid lines are derived from the lattice simulations of the single-spectator-field scenario.Lighter lines correspond to smaller quartic self-couplings (down to λ = 10 −5 ) while darker lines to larger couplings, up to λ = 10 −1 .large band of momenta that cannot be properly covered by a lattice with N = 512, the largest one we are currently able to simulate.Therefore, we leave the testing of this enticing scenario to a future work.

Conclusions and outlook
In the present article, we have carried out an extensive analysis of the Hubble-induced heating mechanism appearing in quintessential inflation scenarios involving spectator fields nonminimally coupled to gravity.By performing a large number of 3+1-dimensional classical lattice simulations, we have probed a wide range of the parameter space of the theory, numerically determining the main timescales in the process, namely the initial tachyonic insta-  bility, the onset of backreaction effects and the virialisation phase leading to a radiation-like equation-of-state for the spectator field.These findings allowed us to derive a set of parametric equations that can be conveniently applied to evaluate the efficiency of the heating process, the temperature at the onset of radiation domination and the minimum number of e-folds of inflation required to solve the flatness and horizon problems in specific quintessential scenarios, without relying on further lattice simulations.In this sense, our results constitute a major step forward in addressing the question of heating in non-oscillatory models of inflation [42,58,[61][62][63], within a minimal, natural and non-perturbative setup.In particular, since the gravitational coupling considered in this work is required by the renormalisation of the energy-momentum tensor in curved space-time [29,64], findings provide a lower bound on the heating efficiency and the associated radiation temperature, allowing generically for successful big bang nucleosynthesis even in the absence of direct couplings between the inflaton and the matter sector.
For the sake of completeness, we also explored a prototypical yet significant extension of the single-field scenario.By introducing an explicit coupling with a secondary scalar sector beyond the standard model, we have assessed the potential depletion of the spectator field fluctuations prior to the onset of non-linearities.The two-field scenario constitutes an excellent testing ground for the robustness of the single-scalar-field fitting formulas.Indeed, our analysis shows that, for moderate values of the coupling between fields, the results of the single-field model constitute an accurate description of the two-fields model as well.Note also that these results can be directly applied to scenarios involving Abelian gauge bosons as daughter fields, as first shown in Ref. [37] in the context of parametric resonance and explicitly verified here.
Our approach could be easily adapted to scenarios far beyond the toy-model cases con-sidered in this work.The most direct application would be the Standard Model itself, with a Higgs field stable up to the Planck scale [34] playing the role of the spectator field. 2 In this context, the analysis of secondary interactions becomes extremely important, since the tachyonic productions of Higgs perturbations is expected to compete with the secondary amplification of gauge bosons and their potential decay into fermions, along the lines of the Combined Preheating scenario developed in Refs.[35,66,67].Interestingly enough, the inclusion of these dissipative decay channels could strongly modify the single-field picture presented in this paper, opening the door to stabilising the defects generated during the transition and modifying with it the associated gravitational waves spectrum [24,68].Further applications of our setting could include models with a larger field content, involving, for instance, additional dark matter or beyond-the-standard sectors, as well as higher-dimensional operators.In this context, it would be interesting to explore the early nonthermal production of dark matter relics [69][70][71] and the generation of the matter-antimatter asymmetry in the Universe [19].Finally, directly coupled gauge fields, involving for instance Chern-Simons interactions, could also be considered as interesting secondary fields with a rich phenomenology, e.g. in the context of primordial magnetogenesis.

Figure 1 :
Figure1: Evolution of the volume-averaged kinetic (K), gradient (G), potential (V ) and interaction (I) contributions in(2.23)  to the total energy density ρ χ for exemplary values ν = 10 and λ = 10 −4 .Each energy component is multiplied by a 4 to highlight asymptotic radiation-like behaviours.In order to facilitate the comparison, we display the absolute values of the interaction and the total-energy terms, since they are not positive definite at all times.The timescales identifying the backreaction time z br and the virialisation time z rad for this specific case are indicated at the top of the image.

Figure 2 :
Figure 2: Parametric dependence of the backreaction timescale z br (upper panel), of the amplitude of the scalar field fluctuations ⟨χ 2 br ⟩ rms (middle panel) and of the total energy at backreaction ρ χ br (lower panel) on the model parameters ν and λ, as indicated in the figure.In the upper panel, dashed coloured lines correspond to the analytical curve obtained from Eq. (2.25).In the middle and lower panels, solid lines are derived from the numerical output of lattice simulations while the black dashed lines correspond to the fitting formulas (2.26) and (2.27).

Figure 3 :
Figure 3: Correction to the squared Planck mass given by the peak mode enhancement at backreaction ⟨χ 2 br ⟩ rms as a function of the model parameters.The coloured lines set the contours for a 10% correction at different values of the Hubble scale at the onset of kination.The grey shading indicates the duration of the tachyonic phase in e-folds N br = 1/2 ln(1 + z br /ν).The dashed red rectangle delimits the range explored in our extensive scanning of the parameter space.

Figure 4 :
Figure 4: Parametric dependence of the backreaction timescale z rad (upper panel), the total energy at completion of virialisation ρ χrad (lower panel) on the model parameters ν and λ, as indicated in the figure.Solid coloured lines correspond to the numerical results obtained from the full set of lattice simulations.In the lower panel, black dashed lines correspond to the fitting formula(2.33)   .

Figure 5 :
Figure 5: Ratio of the energy density stored in gradients and the total energy density of the spectator field at z rad , as a function of the model parameters ν and λ.

Figure 6 :
Figure 6: Parametric dependence of the heating efficiency Θ ht on the model parameters λ and ν for a fiducial energy scale of kination H kin = 10 11 GeV.Solid lines correspond to the values computed from the simulations scanning of the parameter space, while the dashed lines correspond to the fitting formula (3.2).

Figure 7 :
Figure 7: Predicted values of tensor-to-scalar ratio and spectral tilt in an α-attractor quintessential-inflation scenario involving only gravitational particle production.Lines of different colours correspond to different values of the H kin parameter, where the end points correspond to the largest (left) and smallest (right) heating efficiencies achieved in our simulations.The vertical grey lines delimit the region allowed by the numerically-obtained heating efficiency.The results are shown in comparison with the constraints from Planck + Bicep/Keck at one and two sigma.

Figure 8 :
Figure 8: Energy components of the model as defined in (2.23).Each component has been multiplied by a 4 in order to highlight the radiation-like behaviour.In this figure, the model parameters have been set to ν = 20, λ = 10 −3 , g = 10 −3 .

Figure 10 :
Figure 10: Parametric dependence of the heating efficiency Θ ht on the two model parameters λ and ν.Dashed lines correspond to the fitting formulas obtained for the single-scalar-field case.The colour coding matches the one defined for the previous figures.