Nonlinear Dynamics of Preheating after Multifield Inflation with Nonminimal Couplings

We study the postinflation dynamics of multifield models involving nonminimal couplings using lattice simulations to capture significant nonlinear effects like backreaction and rescattering. We measure the effective equation of state and typical timescales for the onset of thermalization, which could affect the usual mapping between predictions for primordial perturbation spectra and measurements of anisotropies in the cosmic microwave background radiation. For large values of the nonminimal coupling constants, we find efficient particle production that gives rise to nearly instantaneous preheating. Moreover, the strong single-field attractor behavior that was previously identified persists until the end of preheating, thereby suppressing typical signatures of multifield models. We therefore find that predictions for primordial observables in this class of models retain a close match to the latest observations.

Introduction.-Postinflation reheating plays a critical role in our understanding of the very early Universe (see Ref. [1] for a recent review). By the end of the reheating phase-and before big bang nucleosynthesis (BBN) can commence [2]-the Universe must achieve a radiationdominated equation of state and become filled with (at least) a thermal bath of standard model particles at an appropriately high temperature. Although the earliest stages of reheating can be studied within a linearized approximation, some of the most critical processes arise from nonlinear physics, including backreaction and rescattering among the produced particles.
In addition to setting appropriate conditions for BBN, the reheating phase plays a critical role in comparisons between inflationary predictions and recent high-precision measurements of the cosmic microwave background (CMB). In particular, if there were a prolonged period after inflation before the Universe attained a radiationdominated equation of state (EOS), that would impact the mapping between perturbations on observationally relevant length scales and when those scales first crossed outside the Hubble radius during inflation [3][4][5][6]. Residual uncertainty on the duration of reheating, N reh , is now comparable to statistical uncertainties in measurements of CMB spectral observables. Hence understanding the timescale N reh is critical for evaluating observable predictions from inflationary models.
In this Letter we study the nonlinear dynamics of the early preheating phase of reheating in a well-motivated class of models. These models include multiple scalar fields, as typically found in realistic models of high-energy physics [7,8], and each scalar field ϕ has a nonminimal coupling to the spacetime Ricci curvature scalar R of the form ξϕ 2 R. Such nonminimal couplings are quite generic: they are induced by quantum corrections for any selfinteracting scalar field in curved spacetime, and they are required for renormalization [9,10]. Moreover, the dimensionless coupling constants ξ grow with energy scale under renormalization-group flow, with no UV fixed point [11]. Hence they can attain large values at inflationary energy scales. Upon transforming to the Einstein frame, such models feature curved field-space manifolds [12].
Multifield models with nonminimal couplings naturally yield a plateaulike phase of inflation at large field values, of the sort most favored by recent observations [13]. During inflation the fields generically evolve within a single-field attractor, thereby suppressing typical multifield effects that could spoil agreement with observations, such as large primordial non-Gaussianities and isocurvature perturbations [14][15][16].
Previous work, which studied the onset of preheating in this class of models semianalytically, identified three regimes that yielded qualitatively distinct behavior: ξ ≲ Oð1Þ, ∼Oð10Þ, and ≳Oð10 2 Þ [17][18][19]. In this Letter we significantly expand this work, employing lattice simulations to study the complete preheating phase, deep into the nonlinear regime. We restrict attention to coupled scalar fields, and neglect the production of standard model particles such as fermions or gauge fields [20][21][22][23][24][25][26][27][28][29][30][31][32][33]. Nonetheless, we are able to analyze the typical timescales required for the Universe to achieve a radiation-dominated EOS, for the produced particles to backreact on the inflaton condensate, ultimately draining away its energy, and for rescattering among the particles to yield a thermal spectrum. For large couplings, ξ ≳ 10 2 , of the sort encountered in Higgs inflation [34], we find very efficient preheating, typically completing within the first two e-folds after the end of inflation, thereby protecting the close match between predictions for primordial observables and the latest CMB measurements.
Model.-In the Jordan frame, the nonminimal coupling between the N scalar fields and the spacetime Ricci scalar R remains explicit in the action through the term fðϕ I ÞR. Upon rescalingg μν ðxÞ → g μν ðxÞ ¼ Ω 2 ðxÞg μν ðxÞ, with pl , we transform the action into the Einstein frame. (Here M pl ≡ 1= ffiffiffiffiffiffiffiffiffi 8πG p ¼ 2.43 × 10 18 GeV is the reduced Planck mass.) The Einstein-frame potential is stretched by the conformal factor, Vðϕ I Þ ¼Ṽðϕ I Þ=Ω 4 , compared to the Jordan-frame potentialṼðϕ I Þ. Taking canonical scalar fields in the Jordan frame, the nonminimal couplings induce a curved field-space manifold in the Einstein frame, with field-space metric given by G IJ ðϕ K Þ ¼ ½M 2 pl =ð2fÞfδ IJ þ 3f ;I f ;J =fg [12]. The equation of motion for the fields in the Einstein frame is then where Γ I JK ðϕ L Þ is the Christoffel symbol constructed from G IJ . We consider an unperturbed, spatially flat Friedmann-Lemaître-Robertson-Walker (FLRW) spacetime metric, so the Einstein field equations yield H 2 ðtÞ ¼ ρ total =ð3M 2 pl Þ, where ρ total is the total energy density of the system, HðtÞ ≡ _ a=a, and overdots denote derivatives with respect to cosmic time.
We consider two-field models, ϕ I ¼ fϕ; χg, with The topography of the Einstein-frame potential generically includes "ridges" and "valleys" along certain directions χ=ϕ ¼ const. For non-fine-tuned parameters, the fields quickly fall to a local minimum (valley) of the potential, and the background dynamics obey a strong "single-field attractor" [15][16][17]. For symmetric couplings, with ξ ϕ ¼ ξ χ and λ ϕ ¼ g ¼ λ χ , any initial angular motion within field space damps out within a few e-folds after the start of inflation, and the system flows toward the minimum of the potential along a single-field trajectory [35]. Within a single-field attractor, the predictions for the spectral index n s , the tensor-to-scalar ratio r, the running α ¼ dn s =d ln k, primoridal non-Gaussianities f NL , and isocurvature perturbations β iso remain consistent with the latest observations across large regions of phase space and parameter space [15][16][17]. Field fluctuations in these models are sensitive to the curvature of the field-space manifold, which is greatest near the origin. During preheating, as the inflaton condensate oscillates through zero, the effective mass for the fluctuations δχ receives quasiperiodic "spikes" proportional to a component of the field-space Riemann tensor. In the limit ξ I ≫ 1, these scale as R χ ϕϕχ ∝ ξ ϕ . These large spikes lead to sharp violations of the adiabatic condition for those modes, driving efficient particle production [17][18][19]36].
Results.-We employ a modified version of Grid and Bubble Evolver (GABE) [37] to evolve the fields and the background, according to Eq. (1) and the Friedmann equation. Whereas the original software was used to simulate nonminimally coupled degrees of freedom (d.o.f.) [38], we have modified the code significantly to allow for a curved field-space metric in both the dynamics of the fields as well as the initial conditions. We start the simulations when inflation ends, defined by ϵðt init Þ ¼ 1 where ϵ ≡ − _ H=H 2 ; the Hubble scale at this time is H end . We use a grid with N ¼ 256 3 points and a comoving box size L ¼ π=H end so that the longest wavelength in our spectra corresponds to k ¼ H end =2. We match the twopoint correlation functions of ϕðt init ; xÞ and χðt init ; xÞ to corresponding distributions for quantized field fluctuations. Fourier modes of the quantized fluctuations evolving during inflation within the single-field attractor may be parametrized as δϕ I k ¼ ffiffiffiffiffiffi G II p v I k ðτÞ=aðτÞ (no sum on I), where dτ ≡ dt=aðtÞ is conformal time [17]. Near the end of inflation, we use the Wentzel-Kramers-Brillouin (WKB) approximation to estimate amplitudes jv I k ðτ init Þj ¼ ½2Ω ðIÞ ðk; τ init Þ −1=2 , where Ω 2 ðIÞ ðτÞ ¼ k 2 þ a 2 ðτÞm 2 eff;I ðτÞ. The effective masses m 2 eff;I include distinct contributions from the curvature of the potential and from the curvature PHYSICAL REVIEW LETTERS 123, 171301 (2019) of the field-space manifold, and are analyzed in detail in Refs. [17][18][19]. (Here we neglect contributions from coupled metric perturbations). The initial spectra of the fields are subject to a window function that suppresses highmomentum modes above some UV suppression scale, k UV ¼ 50H end . Figures 1 and 2 show results for case A with ξ ϕ ¼ 10, 100. In Fig. 1, we plot the evolution of the inflaton condensate after the end of inflation as calculated in a linearized treatment (akin to Ref. [19]), and as calculated from the spatial average hϕi on the lattice. Backreaction of produced particles-which is absent in linearized analysesbecomes significant beginning around 2.7 e-folds after the end of inflation for ξ ϕ ¼ 10. For ξ ϕ ¼ 100 backreaction is strong enough to completely drain the inflaton condensate within the first 2 e-folds. Figure 2 shows the evolution of the peak values of the spatial averages hϕi and hχi as well as the growth of fluctuations, characterized by ϕ rms ≡ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi hϕ 2 i − hϕi 2 p and χ rms ≡ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi hχ 2 i − hχi 2 p . (Growth of field fluctuations corresponds to particle production [1].) We have confirmed that the early growth of δϕ and δχ fluctuations in our lattice simulations closely matches the behavior calculated via Floquet analysis in Ref. [18]. Beginning around 2.6 e-folds, nonlinear rescattering among the δχ fluctuations drives rapid growth of the δϕ fluctuations for ξ ϕ ¼ 10. For ξ ϕ ¼ 100 the same effect occurs within the first e-fold. Backreaction and rescattering generally become significant at distinct times as one varies couplings [39].
The dynamics of the δϕ and δχ fluctuations vary with coupling ξ ϕ , as shown in Fig. 3. For ξ ϕ ¼ 1, 10 parametric resonance due to the contribution from the potential to m 2 eff;χ leads to a slow growth of δχ fluctuations; these eventually rescatter, leading to the growth of δϕ fluctuations and lowering the χ rms =ϕ rms ratio. For ξ ϕ ≥ 40 the "Ricci spike" [17,36] leads to a fast growth of δχ fluctuations. This is seen in Fig. 3 as an early rise of the χ rms =ϕ rms ratio. When χ rms grows enough it rescatters with δϕ fluctuations, eventually leading to χ rms =ϕ rms ∼ 1. The case of ξ ϕ ¼ 25 is the most interesting, since it displays several distinct phases. The initial growth occurs due to adiabaticity violation caused by the Ricci spike. After 1.5 e-folds the height of the Ricci spike has redshifted, making it comparable to the potential contribution to the effective mass, thereby shutting off particle production [17]. When the Ricci spike redshifts even more, around 2.5 e-folds, a second stage of parametric resonance commences, due to the potential term alone. Subsequently, rescattering enhances the δϕ fluctuations, lowering the χ rms =ϕ rms ratio. The situation is qualitatively similar for the symmetric case (B) [39].
The rapid growth of fluctuations yields an efficient transfer of energy from the inflaton condensate into radiative d.o.f. Within the single-field attractor, we may approximate the energy density in the inflaton condensate as [17] where we evaluate G ϕϕ with ϕ → hϕi and χ ∼ 0. Figure 4 shows that across cases A and B the fraction of energy density in the inflaton condensate falls sharply within the first few e-folds after the end of inflation; for ξ ϕ ≥ 100, virtually all of the energy density has been transferred out of the inflaton condensate within the first N ¼ 1.5 e-folds.
The rapid transfer of energy to radiative d.o.f. is similarly reflected in Fig. 5, which shows the evolution of the EOS, w ¼ p total =ρ total , where ρ total and p total are the total energy density and pressure for the system, respectively. In this case, the system approaches w ¼ 1=3 rapidly for small couplings ξ ϕ ∼ Oð1Þ, because in that regime the Einsteinframe potential for the inflaton approximates a quartic form, so that even the condensate's oscillations correspond to w ≃ 1=3 [17]. As ξ ϕ increases, the Einstein-frame potential for ϕ approaches a quadratic form, for which the condensate's oscillations behave like w ≃ 0 [17], but in that case, the stronger coupling yields more efficient particle production, so that the system eventually becomes dominated by radiative d.o.f. For ξ ϕ ¼ 100, we find a transient phase with a stiff EOS, w > 1=3, which likely arises because typical momenta for the fluctuations are comparable to m eff;I , and the contributions to ρ total and p total from kinetic and spatial-gradient terms are weighted by components of G IJ , which are significant for ξ ϕ ≫ 1. At later times, as m eff;I → 0, the system relaxes to a gas of massless particles with w ¼ 1=3. Across a wide range of couplings for this family of models, we therefore find that the Universe rapidly achieves a radiation-dominated EOS within N rad ∼ 2-2.5 e-folds after the end of inflation. Preheating in α-attractor models with α ¼ Oð1Þ, in contrast, can lead to a prolonged period with w ≃ 0 [40], shifting the pivot scale accordingly and thereby offering a means to empirically distinguish between such models and the family we consider here.
The strong rescattering among fluctuations yields an efficient start to the process of thermalization, by transferring power between particles of different momenta. In Fig. 6 we show the spectra in field fluctuations δϕ and δχ for case A with ξ ϕ ¼ 10. Although the spectra are dominated at early times by increased power in distinct resonance bands, by later times rescattering has flattened out the distributions for both δϕ and δχ. By N therm ¼ 2.8 e-folds after the end of inflation, both fields have attained a spectrum consistent with a thermal distribution, jδϕ I k j 2 ∝ ½kðexp½k=T − 1Þ −1 , at a temperature T reh ∼ OðH end Þ. We find comparable behavior across cases A and B for ξ ϕ ≥ 1 [39].
The rapid thermalization means that the system reaches the adiabatic limit soon after the end of inflation. We denote N ad ¼ min½N bg ; N therm , where N bg is the time by which  super-Hubble coherence of the inflaton condensate is lost, indicated by ϕ rms > hϕi. Any significant turning of the system within the field space between the end of inflation and N ad could amplify non-Gaussianities and isocurvature perturbations, thereby threatening the close agreement between predictions in these models and measurements of the CMB [41][42][43][44]. In Fig. 7, we plot ω=H across cases of interest, where ω ¼ jω I j is the covariant turn rate [45]. Even as the Hubble rate falls over time, we nonetheless find ω=H < 0.1 through N ad , indicating minimal turning of the system within field space.
Our late-time results were unchanged as we varied the initial UV suppression scale k UV ¼ bH end between b ¼ 25, 50, and 100, and the number of grid points between 128 3 , 256 3 and 512 3 . We discuss this and related numerical convergence tests in Ref. [39].
Conclusions.-Multifield models of inflation with nonminimal couplings generically yield predictions for primordial observables in close agreement with the latest observations, deriving from the strong single-field attractor behavior of these models [15][16][17]. Throughout the cases we have examined and across parameter space, we find that this single-field attractor behavior remains robust until the system reaches the adiabatic limit after inflation, with no significant turning in field space even in the midst of strongly nonlinear dynamics.
Preheating in this class of models is efficient, draining the energy density from the inflaton condensate within N bg ≲ 1.5 e-folds in the limit of strong couplings, ξ I ∼ 100. The system typically reaches a radiation-dominated equation of state within N rad ≲ 2.5, while rescattering yields a rapid onset of thermalization within N therm ≲ 3, thereby fulfilling several of the most critical requirements of the reheating phase. We defer to future work such questions as possible impact of coupled metric perturbations on the fully nonlinear preheating dynamics, and the coupling of the scalar fields ϕ and χ to standard model particles.