Parametric Resonance in the Early Universe – A Fitting Analysis

Parametric Resonance in the Early Universe – A Fitting Analysis

[    [

Particle production via parametric resonance in the early Universe, is a non-perturbative, non-linear and out-of-equilibrium phenomenon. Although it is a well studied topic, whenever a new scenario exhibits parametric resonance, a full re-analysis is normally required. To avoid this tedious task, many works present often only a simplified linear treatment of the problem. In order to surpass this circumstance in the future, we provide a fitting analysis of parametric resonance through all its relevant stages: initial linear growth, non-linear evolution, and relaxation towards equilibrium. Using lattice simulations in an expanding grid in dimensions, we parametrize the dynamics’ outcome scanning over the relevant ingredients: role of the oscillatory field, particle coupling strength, initial conditions, and background expansion rate. We emphasize the inaccuracy of the linear calculation of the decay time of the oscillatory field, and propose a more appropriate definition of this scale based on the subsequent non-linear dynamics. We provide simple fits to the relevant time scales and particle energy fractions at each stage. Our fits can be applied to post-inflationary preheating scenarios, where the oscillatory field is the inflaton, or to spectator-field scenarios, where the oscillatory field can be e.g. a curvaton, or the Standard Model Higgs.

a]Daniel G. Figueroa, b]Francisco Torrentí Prepared for submission to JCAP

Parametric Resonance in the Early Universe – A Fitting Analysis

  • Theoretical Physics Department, CERN, Geneva, Switzerland

  • Instituto de Física Teórica IFT-UAM/CSIC, Universidad Autónoma de Madrid, Cantoblanco 28049 Madrid, Spain.




1 Introduction

Compelling evidence supports the idea of an inflationary phase in the early Universe [1]. The specific particle physics realization of inflation is however uncertain, so the inflationary period is typically parametrized in terms of a scalar field, the inflaton, with a vacuum-like potential. After inflation, the reheating stage follows, converting all inflationary energy into different particle species, which represent all the matter and radiation in the Universe. Eventually, the created particles dominate the total energy budget and ’thermalize’, signaling the onset of the ’hot Big Bang’ thermal era.

In this paper we consider inflaton potentials with simple monomial shapes, as this gives rise to one of the most important particle creation phenomena in the early universe: parametric resonance. This is the case of chaotic inflation models, where the inflaton rolls down a monomial potential during the whole inflationary period. Although these scenarios are under tension with cosmological data [1], the simple addition of a small non-minimal gravitational coupling reconcile them with the observations [2]. Some scenarios which fit perfectly well the observational data, e.g. Higgs-Inflation [3, 4] and Starobinsky inflation [5], also exhibit a monomial potential with a single minimum, but only during the stages following inflation.

In all the scenarios we consider, soon after the end of inflation, the inflaton is in the form of a homogeneous condensate, and starts oscillating around the minimum of its potential. Each time the inflaton crosses zero, all particle species sufficiently strongly coupled to the inflaton, are created in energetic bursts. In the case of bosonic species, the production of particles is resonant, and the energy transferred grows exponentially within few oscillations of the inflaton [6, 7, 8, 9, 10, 11, 12, 13]. In the case of fermionic species, there is also a significant transfer of energy [14, 15, 16, 17], but Pauli blocking prevents resonance from developing. The production of particles in this way, either of fermions or bosons, represents the archetypical example of what is meant by an initial ’preheating’ stage of reheating.

Inflationary preheating is however not the only case where parametric resonance takes place in the early Universe. If a light spectator field is present during inflation, this field forms a homogeneous condensate during the inflationary period, and oscillates around the minimum of its potential afterwards. This is the case e.g. of the curvaton scenario [18, 19, 20, 21]. The curvaton may decay after inflation via parametric resonance, transferring abruptly all its energy to the particle species coupled to it [22, 23, 24, 25]. Another example of a spectator field, naturally decaying through parametric resonance after inflation, is the Standard Model Higgs field. If the Higgs is weakly coupled to the inflationary sector, the Higgs is always excited either during inflation [26, 27, 28], or towards the end of it [29, 30]. The Higgs is then ’forced’ to decay into the rest of the SM species after inflation111Note that in the case of Higgs-Inflation [3, 4], the Higgs also decays after inflation via parametric resonance, into the rest of the SM fields [31, 32, 33]. In this scenario, the Higgs plays however the role of the inflaton. Therefore, the Higgs decay in the case of Higgs-inflation scenarios [31, 32, 33, 34], should rather be categorized within the context of preheating scenarios., via parametric resonance [27, 35, 36, 37, 38, 39, 30].

In this paper, independently of the context, we will often refer to the oscillatory field as the ’mother’ field, and to the created species as the ’daughter’ fields. Particle production of daughter fields via parametric resonance, corresponds to a non-perturbative effect, which cannot be captured by perturbative coupling expansions, not even if the couplings involved are small [10]. During the initial stage of parametric resonance, the system is linear, and analytical methods can be applied. As the particle production is exponential for bosonic species, the daughter field(s) eventually ’backreact’ onto the mother field, making the system non-linear. In order to fully capture the non-linearities of the system, we need to study this phenomenon in the lattice. The approach of classical field theory real-time lattice simulations can be considered valid as long as the occupation number of the different species is much larger than one, and hence their quantum nature can be ignored [40, 41]. Lattice simulations have been, in fact, successfully carried out for different preheating scenarios during the last years, see e.g. [42, 43] and references therein. However, each time a new scenario exhibits parametric resonance, a new re-analysis is often required.

As lattice simulations are computationally expensive and time consuming, and not everybody has the expertise on the appropriate numerical packages [44, 45, 46, 47, 48], many studies often resort to over-simplified analytical analysis, which capture only the initial linear stage. A systematic study of parametric resonance, fitting the dynamics through all the relevant stages, from the initial linear growth till the relaxation towards equilibrium, passing through an intermediate non-linear stage, is missing in the literature. In this work, we fill in this gap. We have used massively parallelized lattice simulations to charaterize the dynamics of parametric resonance through all its stages. We have parametrized the dynamics by scanning over the relevant circumstances and parameters: role of the oscillating field, particle coupling, initial conditions, and background rate of expansion. We have obtained in this way simple fits to the most significant quantities, like the characteristic time scales and energy fractions of the different particle species. Our fitted formulas can be applied to the study of parametric resonance in scenarios where the mother field dominates the energy budget of the universe (i.e. preheating), or in scenarios where the mother field represents only a sub-dominant component (e.g. inflationary spectator fields).

As parametric resonance in the context of the early Universe has been well studied in the past, let us emphasize here that with our present work, we simply aim to aid in the analysis of future scenarios exhibiting parametric resonance. The advantage of using our fitted formulas will be twofold: on the one hand skipping the tedious task of running new simulations, and on the other hand preventing the use of over-simplified linear analysis of the problem.

The structure of this work is as follows. In Section 2 we describe general aspects of parametric resonance, while we derive an analytical estimate of the decay time of the mother field, based on a linear calculation. In Section 3 we preset the numerical results from our lattice simulations. We describe our results for preheating with a quartic potential in Section 3.1, and for preheating with a quadratic potential in Section 3.2. We compare these results against the analytical estimations from Section 2. In Section 3.3 we present the analogous numerical study for scenarios where the mother field represents only a sub-dominant energy component of the Universe. In Section 4 we list all fitted formulas together from all scenarios considered. In Section 5 we discuss the context where our results can be useful. In the appendices we present details on the lattice formulation we have used, and discuss briefly the evolution of the field spectra in some of the scenarios considered.

From now on we consider units, and represent the reduced Planck mass by GeV. We take a flat background with Friedman-Robertson-Walker (FRW) metric , where is the scale factor, and the cosmic time.

2 Parametric Resonance: Analytical Calculation

Before we move into the specific scenarios of Sections 3.1, 3.2, and 3.3, let us discuss some general aspects of parametric resonance, while we derive an analytical estimation of the decay time of the mother field. In the following Sections we will compare this analytical estimation with the results obtained from lattice simulations.

Let us begin by considering a scalar field with a quartic potential , coupled to another scalar field through an interaction , with a dimensionless coupling constant. The equations of motion (EOM) of the system read


where is the Hubble rate. We will consider the field to be initially homogeneous with some initial amplitude , and null initial velocity , whilst the field is not excited initially, . If we neglect for the time being the interaction term, the equation for the homogeneous part of the field, corresponds to an anharmonic oscillator in the presence of a friction term. As has a single minimum at , the field will start rolling down towards the minimum. If the friction term dominates over the potential term, the system is overdamped and the field rolls down very slowly, in the so called slow-roll regime. Eventually, as the Hubble rate diminishes due to the expansion of the Universe, there will be a time when the system becomes underdamped. This time signals the onset of the mother field oscillations around the minimum of its potential. More specifically, we will define an initial time as the moment when the Hubble rate just becomes smaller than the effective frequency of oscillation. In light of Eq. (2.1), the period of oscillation is , so we can determine from the condition , with and .

After a convenient conformal transformation of the time and field variables


the EOM read


where , , and is the so called resonance parameter,


Neglecting the interaction term, the EOM of the homogeneous part of reduces to


In the case when the mother field dominates the energy budget of the universe (e.g. in preheating), the energy density scales as radiation dominated (RD) [49], so the scale factor behaves as . In this case, the term on the of Eq. (2.6) simply vanishes, . If the field does not dominate the energy budget of the universe, the behavior of the scale factor depends on the equation of state of the dominant energy component of the Universe. For fixed , one can find , which either dies away as if , or vanishes directly for . We will therefore set , for the simplicity of the discussion. The solution of Eq. (2.6) with initial conditions , , is the Elliptic function222In reality, the initial conditions should be , , with some value propagated from the past when the field was deep in the slow-roll condition . Taking into account this does not change the essence of the oscillatory regime once the field enters into the underdamped regime. Hence, for the easiness of the discussion, we will simply stick here to the solution .


The equation for the Fourier modes of the field (assuming RD) can be written as


In this form the equation for the fluctuations of the field does not depend on the expansion of the universe, and it is completely reduced to a problem in Minkowski space-time333This is of course only a special feature of the conformally invariant theory .. Given the behavior of in Eq. (2.7), Eq. (2.8) corresponds to the class of the Lamé equations, which has a well-understood structure of resonances. Whenever , with (i.e. [1, 3], [6, 10], …), there is an infrared band of modes , for which the modes can be exponentially amplified as , with a parameter known as the Floquet index [11]. Considering the mode frequency , we can speak of adiabatic modes if the condition is fulfilled. The set of unstable modes correspond to the modes that violate the adiabaticity condition each time crosses around zero, verifying the opposite condition, . The instability of the resonant modes is naturally interpreted as a strong particle creation of the field, as the occupation number grows as .

Figure 1: Left: We show the stability/instability chart of the Lamé equation (2.8). Coloured bands indicate the regions of the (,) parameter space in which the real part of the Floquet index is a positive number and hence the solution of the Lamé equation is exponential. The darker the colour, the greater the index, up to a maximum of for black areas. White areas are the regions in which . Right: Some examples of the Floquet index derived numerically from the Lamé equation for resonance parameters ranging between and . In each panel, we plot the corresponding Floquet index as a function of the momentum . We have divided the different ’s in two groups: those inside one of the resonance bands , , , which excite modes down to (blue solid lines), and those which are in between resonance bands (red dashed lines), which only excite modes down to some minimum momentum .

If the resonance parameter is not within one of the resonant bands, but lies in between two adjacent bands, then there is still a resonance of the type , but within a shorter range of momenta , and hence with a smaller Floquet index . There is a theoretical maximum value for the Floquet index given by  [11], so that any is always constrained as for . For resonant parameters , is typically of order , see Fig. 1.

For simplicity, in the remaining of this Section we will consider the resonance parameter to be within one of the resonant bands, , , . The growth of the fluctuations in the initial stages of resonance is described by the linear Eq. (2.8). Even if the amplitude of the fluctuations grows exponentially, Eq. (2.8) is expected to represent a good description of the field excitation during the initial stages. Of course, one is ignoring in this way the backreaction of the bosons into . This is a good approximation for as long as the energy tranferred into the field is only a marginal fraction of the energy available in the mother field . In our numerical analysis of Section 3.1, we will quantify exactly when the linear approximation breaks down. For the time being, to continue with our analytical approach, we will just consider valid the linear regime all the time through.

The energy density of the created particles due to the resonance, is given by


where we have introduced an oscillation-averaged effective mass for the field,


with the oscillation period of  [11]. From the violation of the adiabaticity condition for , i.e. , we can determine an estimation of the maximum (comoving) momentum possibly excited in broad resonance,


where we have identified . From Eqs. (2.10),(2.11), we conclude that


In other words, in broad resonance , the decay products are always non-relativistic. Correspondingly we can approximate the effective mode frequency as , where . If is within a resonant band, then all modes with momenta are excited with some Floquet index varying within . This corresponds to the cases with blue solid lines in Fig. 1. We can therefore model the occupation number of the excited modes simply as a step function , with a mean Floquet index. It follows that444Notice that the scaling is characteristic of relativistic species, despite the fact that we stated that the decay products are non-relativistic. This is because the energy density of the daughter fields is given by , with the number density and their mass, as it corresponds to any non-relativistic species. However, while , the effective mass is also time dependent, , and hence the total energy density scales as radiation .


This is how the energy density of the daughter fields (those fully within a resonant band) will grow, at least as long as their backreaction into the mother field remains negligible. Using this linear approximation we can estimate the moment at which an efficient transfer of energy has taken place from into the bosons, characterized by . This will be just a crude estimate of the time scale of the mother field decay, since by then backreaction and rescattering effects will have become important, invalidating the linear approach. However, the nonlinear effects due to backreaction of the decay products, simply tend to shut off the resonance. Hence, the calculation in the linear regime should provide at least, in principle, a reasonable estimate of the time scale for when the energy has been efficiently transferred into the daughter fields. Whether is also a good estimate of the decay time of the mother field, will be contrasted against our lattice simulations in the next Section.

The energy of the oscillating field, since the onset of the oscillations, decays as [38]


where in the second equality we have used . We can now find by simply equating Eqs. (2.13) and (2.14),


so that


For instance, for chaotic inflation with quartic potential, , and hence . Looking at Fig. 1, we see that the Floquet index of the modes for which is within a resonant band (blue solid lines in the figure), can be approximated, as said, by a simple step function , with a mean Floquet index . Taking this into account, for chaotic inflation we find


It is clear that the larger the , the shorter it takes for the mother field to transfer energy efficiently into the daughter fields. This is expected, as the stronger the interaction is, the faster the decay should be. We see that the decay time, however, according to the above calculation, is always some value of the order . Therefore, contrary to ’popular wisdom’ about parametric resonance, the time scale , identified with the decay of the oscillatory field in the linear approximation, is in practice mostly independent of . Though it is certainly true that the larger the the shorter the decay, the dependence is only logarithmic, see Eq. (2.16), so the time scale does not change appreciably. For instance, increasing in 10 orders of magnitude, only speeds up the decay time in a factor . In the following Section we will check the validity of these estimation by comparing it with the numerical outcome obtained directly from lattice simulations.

Before we move into the numerical results, let us note that a similar computation can be carried out for a mother field with a quadratic potential . The details are more cumbersome in this case, because contrary to the quartic case previously described, in the quadratic case (when the expansion of the Universe cannot be ignored), the Floquet index is not fixed for a given mode. This is because there is now a new mass scale, , which breaks the conformal invariance, making impossible to reduce the problem into a Minkowski analogue 555Of course if there was no expansion of the Universe, the problem is directly formulated in Minkowski, so the structure of the resonance bands is fixed. In such a case, there is a well defined Floquet index for each mode. However, whenever the expansion of the universe cannot be ignored, as it is the case in preheating, each mode scans several resonance bands, and therefore one cannot ascribe a given Floquet index to a given mode., as it happened in the quartic case. In the quadratic case the resonance of a given mode is such that each mode scans several resonance bands, and the evolution of a resonant mode function is in fact stochastic, see [7] for a detailed explanation on this. Without entering into further details, as the linear computation in the quadratic case was carried out in [10], we do not repeat it here. We just quote their result, adapting it to our notation. They find that the maximum momentum excited during parametric resonance in a quadratic potential is approximately


Taking as a reasonable averaged value of the stochastic Floquet index , for chaotic inflation with , Eq. (112) of [10] is equivalent to


with . As in the quartic case, we see that one expects this scale to be always of the order of , changing only logarithmically with resonance parameter.

3 Parametric Resonance: Lattice Simulations

As mentioned before, parametric resonance in the early Universe can be realized in two main different circumstances: when the mother field dominates the energy budget of the Universe, and when the mother field is only a sub-dominant energy component of the Universe. In this Section we will perform lattice simulations of both situations:

  • Inflaton Preheating. In this case we identify the field with the field responsible for inflation, the inflaton. We consider single-field slow-roll scenarios where the inflaton has a monomial potential . Short after inflation ends, when the slow-roll parameters become approximately of order unity, the Hubble rate just becomes smaller than the inflaton mass. As the inflaton has a very large vacuum expectation value (VEV), the inflaton amplitude starts then oscillating around the minimum of its potential. This induces a strong creation of all particles coupled to it, if the coupling strength is sufficiently large. The creation of these particles represents possibly the most important particle creation stage in the history of the Universe: as the inflaton and its decay products are the dominant energy component of the Universe, this stage represent the creation of (most of) the matter in the universe. This adds an extra difficulty, as the time-evolution of the scale factor must be obtained by solving self-consistently the fields EOM together with the Friedmann equations. We consider the two paradigmatic particular models of chaotic inflation, where the inflaton has either a quartic potential (Section 3.1) or a quadratic potential (Section 3.2):


    The strength of the parameters and is fixed by the amplitude of the observed CMB anisotropies. In the quartic model, the energy density of the inflaton scales (after averaging over oscillations) as in a RD background, with , with the scale factor evolving correspondingly as . In the quadratic model the energy density of the inflaton (again after oscillations-averaging) evolves as in a MD background, with , and the scale factor evolving correspondingly as . Both scenarios of inflation are in fact challenged by recent CMB measurements [1] (the quartic case more severely), but in reality, the simple addition of an non-minimal gravitational coupling to the inflaton can easily reconcile these scenarios with the observations [2].

  • Inflationary Spectator Fields. In this second type of scenarios, we consider the field to be just a spectator field during inflation, hence representing a very subdominant component of the energy budget. This does not prevent however the amplitude of these fields to be rather large at the end of inflation (though not as large, in principle, as in single field chaotic inflation scenarios). When inflation ends and the Hubble rate becomes smaller than the effective mass of the spectator field, the amplitude of the field starts oscillating around the minimum of its potential. The expansion rate of the universe after inflation is determined by the inflationary sector, which we will not model explicitly. It is in fact only the evolution of the scale factor that we really need to introduce in the simulations. For instance, for matter-dominated (MD), radiation-dominated (RD) and kination-dominated (KD) universes, the scale factor behaves as , , and , respectively. The most obvious case of a spectator-field is a curvaton, which is normally described with a quadratic potential666Other polynomial potentials have been considered, but the realization of the curvaton mechanism seems much more contrived in those cases [50]. of the type in the context of a RD background [18, 19, 20]. We will restrict our numerical analysis to this case (Section 3.3), taking as a free parameter varied over a certain range. A relevant case of a spectator-field with a quartic potential, although not a curvaton, is the Standard Model (SM) Higgs in the weak coupling limit [26, 27, 28, 29, 30]. The study of the Higgs dynamics after inflation has triggered recently an intense activity [27, 35, 36, 38, 39, 30]. In particular, in [38], the outcome of the dynamics was parametrized in a similar fashion to what we will do here in Section 3.2 for the RD quadratic curvaton case. Therefore, we will not repeat the details of the quartic case here, though we will include a summary of those results in Section 4, where we collect the fits from all the cases studied (inflaton or spectator field cases, with quadratic or quartic potential).

In all scenarios, we will always consider a symmetric interaction between the mother field and the daughter field . This interaction is scale free, with a dimensionless coupling constant. This is particularly convenient from the point of view of the lattice, since any other form of interaction would require the introduction of a new mass scale. Besides, this interaction has been often assumed in the context of preheating, and it is the leading interaction term in the context of gauged spectator fields, as demonstrated in [38] for the case of the SM Higgs. It is also interesting to note that this interaction does not lead to a tree level decay of the mother field into the daughter species, so all the transfer of energy from into will be due only to the non-perturbative effects characteristic of parametric resonance.

3.1 Lattice Simulations of preheating with quartic potential

We consider in this section preheating in the case of a massless self-interacting inflaton with potential


The time for the onset of the oscillatory regime is defined through the condition . This constitutes the initial time of our lattice simulations. We will write all quantities evaluated at time with a sub-index , so this condition can be simply written as . From a simple numerical calculation of the homogeneous Klein-Gordon and Friedmann equations, , , we obtain and . The equations of motion (EOM) of the inflaton and the daughter field can be easily derived, but for convenience, let us first define new field and space-time variables, similarly as in Sec. 2,


where are the old cosmic time and comoving coordinate variables. We denote this set of field and spacetime variables as the ’natural’ variables of the problem. We indicate differentiation with respect cosmic/natural time with a dot/prima respectively, so and . Spatial derivatives should equally be understood, from now on, as taken with respect natural variables. In these variables, the EOM are




is the resonance parameter. These equations are of course the same as Eqs. (2.4) from Sec. 2. However, whereas before, in order to gain some insight on the dynamics of parametric resonance, we used the homogeneous part of the equation for and the Fourier transformed equation of , now we will be rather solving the (lattice version) of the full Eqs. (3.4) in real space.

We take in this Section, as this is fixed by the observed amplitude of the CMB anisotropies [2]. The strength of the coupling is in principle arbitrary. However, in order not to spoil inflation, radiative corrections in the effective inflaton potential must be under control. This sets a constraint  [51]. Unfortunately, in practice we are not capable of simulating resonance parameters outside of the range . Since , this means that we can only simulate couplings . The lower limit is due to the natural limitations of the lattice to simulate fields with narrow resonance bands, as we cannot resolve well the relevant dynamical range of momenta with an appropriate number of modes. The upper limit emerges because the required simulation time and number of lattice points grow with . This is discussed in more detail in Appendix B, so we refer there to the interested reader. Fortunately, as we shall see, the results for the ’s simulated are well described by simple power-law fits, allowing in principle to extrapolate the outcome to larger ’s.

3.1.1 Onset of non-linearities, energy evolution and decay time

Let us briefly recall first the properties of the system from our discussion in Sect. 2. As the mode function of the daughter field follows the Lamé equation [Eq. (2.8)], there are unstable solutions of the type , with the -dependent Floquet index. For certain values of , , causing an exponential growth of the given field mode, and hence of the occupation number. When , the growth of is much stronger than for other values, as can be seen in the pattern of resonance depicted in Fig. 1.

Figure 2: We show the initial oscillations of the volume-averaged conformal amplitude of the inflaton field . We show the cases , , , and for the preheating scenario with quartic potential. We use notation of Eq. (3.3). The dashed vertical red line indicates the time , when backreaction of the daughter fields become relevant, triggering the decay of the inflaton amplitude and energy density (see also Fig. 4).

Let us move now into the results from the lattice simulations. In Fig. 2 we plot the conformal amplitude of the inflaton field for the resonance parameters and . It is clearly appreciated that during a certain number of oscillations, the conformal amplitude of the inflaton remains just constant, like if it was not coupled to the daughter field(s). However, there is a time (which differs for the different ’s) when the amplitude of the conformal inflaton starts decreasing significantly. This is the initial moment when the inflaton starts decaying due to the backreaction from the daughter fields. We shall refer to that time as (the br subindex meaning backreaction)777Let us note that our definition of backreaction differs from the standard condition labeled as ’backreaction’ in the seminal paper [10], which corresponds to the moment when becomes equal to the effective inflaton mass. The latter is a condition that determines the onset of the modulation of the inflaton’s frequency of oscillation. However, we prefer to define the moment of backreaction as the onset of the decay of the (conformal) amplitude of the inflaton, because it is then when the presence of the excited field becomes truly noticeable, and hence the inflaton energy start decreasing significantly.. During the time , the daugther fields have been experiencing parametric resonance, so their energy density has been growing exponentially from initially small quantum fluctuations888See Appendix B for a discussion about the introduction of initial field fluctuations in the lattice.. As the energy flows from the mother field into the daughter fields, at the amount of energy transferred onto the bosons is not anymore a negligible fraction of energy stored in the mother field. Therefore, from then onwards, the (conformal) inflaton amplitude starts to decrease noticeable, see Fig. 2. The time corresponds, in order words, to the onset of the inflaton decay, when the backreaction effects from excited daughter fields become non-negligible. In practice, we have determined as the moment when the (conformal) energy of the mother field drops with respect its initial amplitude.

Figure 3: We depict as a function of for the range . Each point corresponds to the value obtained directly from a lattice simulation, and we have joined the different points with straight lines. Yellow vertical bands indicate the position of the resonance bands of the Lamé equation . The dashed, purple, lower line indicate the estimate [Eq. (3.6)] for values within resonance bands, while the upper one indicates the fit Eq. (3.7) for the relative maxima.

In Fig. 3 we have plotted the different ’s obtained from our simulations, for several resonance parameters in the range . We observe that follows a clear oscillatory pattern, in clear correspondence with the particular structure of resonance bands shown in Fig. 1. In general, the wider the resonance band in the Lamé equation for a given , the shorter is. For those values of emplaced within resonance bands, we find in fact an almost constant value


On the other hand, the behavior of for values outside the resonance bands, i.e. for , is quite different. For values that are in the left extreme of these intervals, i.e. , takes its maximum value, as this corresponds to the right end of a resonance band at , see Fig. 1. We provide the following phenomenological fit to these relative maxima (excluding the particular case ), which we also plot in the Figure,


As increases inside one of the intervals outside the resonance bands, decreases until hitting at the center (more or less) of the nearest resonance band, see Fig. 3. In conclusion, we observe a direct translation of the resonance structure of Fig. 1 into the lattice simulations. This happens because for , the backreaction effects of onto the is negligible, and hence the Lamé equation (2.8) is really at work.

Let us compare now this result with the analytical calculation from Sect. 2. There, using the linear regime, we derived the time scale in Eq. (2.16), and identified it with the decay time of the mother field. However, we see now that this identification is misleading, as rather corresponds to a rough indication of the time scale when the transfer of energy from the mother field to its decay products becomes significant. In other words, it corresponds to the onset of backreaction, which as explained, it only determines the initial moment when the inflaton starts decaying, see Fig. 2. For the range of values shown in Fig. 3, , so the analytical prediction only overestimates in a factor the actual number , found in the simulations at the onset of backreaction. Failing in a factor is not surprising, as the estimation of in Eq. (2.16) involved in fact many approximations. However, the relevant observation to make here is not that can be considered as an order of magnitude estimation of . Rather, the relevant point, is than should not be identified with a decay time, as it rather signals the moment of backreaction, when the linear approximation breaks down. The time scale for determining the end of the transfer of energy from the mother field into the decay products, which we shall identified as the truly ’decay time’ scale of the inflaton, will be referred to as . As we will explain shortly, it corresponds in fact to a much longer time scale, , which cannot be estimated analytically, as the dynamics at become non-linear.

To follow the post-inflationary dynamics in the non-linear regime, it is useful to see how the different contributions to the total energy of the system evolve as a function of time. The total energy can be written as a sum of its different contributions as




where and are the kinetic and gradient energy of the fields , and and are the interaction and potential energies, all written in terms of the natural variables of Eq. (3.3) (i.e. in terms of the field variables and derivatives of these with respect ).

Figure 4: Evolution of the different energy components of the system as a function of time, see Eq. (3.8), for the inflationary scenario , where . Left: We plot for the initial stages of the inflaton decay, and we have indicated with a vertical dashed red line. Right: We plot the same case for later times. To see better how the equipartition regime holds, we have removed the oscillations by taking the oscillation average of the different functions. We have added two new lines that indicate the sums and , see Eq. (3.11).

In the left panel of Fig. 4 we show the evolution of the volume-averaged amplitude of the different energy components of the system. There we can clearly observe how, at first, the inflaton energy dominates the energy budget of the system, alternating between kinetic and potential energies as the oscillations go on. Short after the onset of the simulation, the rest of energies start growing (including the inflaton gradient energy, which indicates the formation of inhomogeneities), becoming very soon an important part of the total energy. At time , these energies have grown enough so that they start backreacting onto the inflaton condensate, inducing its decay (i.e. the decrease of the inflaton kinetic and potential energies). This can also be appreciated in Fig. 2, where from the (conformal) inflaton amplitude starts decreasing significantly.

Let us note that, although the energy fractions at show some scattered dependence on , in reality they are quite independent of the resonance parameter. From the numerical outcome we find

Energy Fractions at :

with the errors , simply reflecting the scattering of energies with . We see from this that at , most of the energy remains yet in the inflaton. However, we also learn that only when of the total energy is already transferred into the daughter field(s), does backreaction really becomes noticeable, making the inflaton amplitude to initiate its decay. The other energy components remain always at sub-percentage levels during , independently of .

At times , the energy components evolve substantially from the given values in Eq. (3.1.1). The energies evolve towards an ’equiparted’ distribution among components, until the system eventually reaches a stationary regime, where the energy components do not change appreciably. This is observed in the bottom panel of Fig. 4, where we have removed the oscillations by taking the oscillation average of the different energies. We observe different equipartition identities for the and fields respectively,


As it can be appreciated in Fig. 4, the second identity holds almost exactly for all times, while the first one only holds for late times (though it is not a bad approximation at earlier times).

From the analysis of the energies we see that a new time scale, much longer than , can be naturally identified with the decay time of the mother field. This scale can be defined by how long it takes the system to relax from into the stationary regime. We shall call the moment when the stationary regime is onset as . It is this time, and not , that signals the true end of the inflaton decay, because it is at that there is no (appreciable) transfer of energy anymore from the inflaton into the daughter field(s). Although the exact definition of is more arbitrary than , we find appropriate to provide an operative definition based on the level of accuracy of equipartition. In particular, at the moment when the inflaton equipartition energy holds at a better level than , i.e. , the inflaton kinetic and gradient energies are stabilized and do not evolve appreciably further, see Fig. 4. The stabilization of the inflaton energy components when equipartition is set to a level is in fact independent of . This is very relevant, as this makes defined in this way, a good indicator of the decay time of the mother field.

Figure 5: Points show the different obtained for different lattice simulations with different values of , for preheating with quartic potential. The dashed line indicates the best fit (3.12).

The relevant property of is that it grows with the resonance parameter , following a simple power-law fit. We show in Fig. 5 the value of as a function of , as extracted from our lattice simulations with different ’s. We obtain the following fit


which we also show in Fig. 5. Note that for , the scales and are not particularly separated, with . This explains why these point must be fitted with a different power law. Note that the inflaton decay takes longer the greater the resonance parameter (i.e. the larger the mother-daughter coupling), which is in principle counter-intuitive. Following the standard logic of the linear calculation, the larger the the shorter the decay time should be. However, once we have learned that ought not identified with the decay time, but with the onset of back-reaction , then the linear logic does not prevail anymore. The reason as to why the truly decay time follows the opposite trend, increasing with , lies on the fact that for the system has become non-linear. Although one would tend to think that the stronger the coupling the faster the stationary regime should be achieved, our lattice simulations – fully capturing the non-linear dynamics – clearly prove the opposite. This was in fact, also noticed already in [38].

As mentioned, we can only obtain our fits for resonance parameters up to due to the limitations of the lattice approach. However there is nothing specially different in the physics of parametric resonance for . Therefore, there is no impediment, in principle, to extrapolate the scaling law Eq. (3.12) to higher ’s.

Let us note that the energy fractions at do not change appreciably any more in our simulations. Some small change should be expected nonetheless, as the system approaches equilibrium. However this is not captured in our simulations. The energy from the end of the inflaton decay onwards are actually rather independent of , given by the fractions

Energy Fractions at :


again with the errors reflecting some (rather random) scattering of the energies with . We see from this that at , the energy is almost ’democratically’ split between the mother and the daughter field(s), though with some more energy stored in the latter, with , , and . At these moments it is also verified the approximate equipartion and .

3.2 Lattice Simulations of preheating with quadratic potential

Let us now consider preheating after chaotic inflation with an inflaton quadratic potential


In this case, we define the onset of the oscillatory regime when the condition holds, which we take as the initial time of our lattice simulations. From a numerical calculation using the homogeneous Klein-Gordon and Friedman equations, , , we find and . Let us define again a set of ’natural’ variables as


where are the old cosmic time and comoving coordinates. As before, we indicate differentiation with respect cosmic/natural time with a dot/prima respectively, and . Spatial derivatives should be understood as taken with respect natural variables, and corresponding momenta will be referred as . The fields’ EOM in these variables are


where the resonance parameter is defined this time as


We take , as this is fixed by the observed amplitude of CMB anisotropies [2].

Let us focus first on the case of a non-expanding universe, so we set and in the equations above. As before, during some time, the particles are very sub-dominant with respect to the inflaton condensate, and hence the effect of their backreaction onto the inflaton can be neglected. During this regime, the mode equation of the daughter fields , corresponds to the so called Mathieu equation [10], which similarly to the set of Lamé equation, is characterized by a well-known structure of resonance bands. More specifically, for some regions in the plane (with ), there is a solution of the type with . One can distinguish two different regimes in the preheating process, depending on the particular value of . If , the narrow resonance regime holds. In this case, the size of the resonance bands is so small that they cannot be well captured in the lattice. On the other hand, if , the system is in a broad resonance regime, and the bands are large enough so that lattice simulations can be applied in this case.

When the expansion of the universe is introduced, the scale factor affects the EOM of in a non-trivial way: even if the system starts in broad resonance with , as the Universe expands, the system rapidly redshifts towards neighboring bands of lower resonance parameter. This is due to the term in Eq. (3.17), which makes the effective resonance parameter to decrease as time goes by. The system does not remain therefore in a single resonance band, but redshifts due to the expansion of the universe. As a consequence, even if the system starts in a broad resonance regime, it can only be maintained as such for some finite time, until it ends up in a narrow resonance regime. For a detailed analysis of the behavior of the mode functions obeying the Mathieu equation both in Minkowski and in an expanding Universe, we recommend to read the seminal work [10]. In our present work we will just focus mostly, from now on, on the outcome from lattice simulations.

Let us note that, in principle, the coupling can be arbitrarily small, so that we could be in the regime of narrow resonance from the very beginning of the oscillations. As we cannot simulate in the lattice narrow resonance, we certainly want to avoid such cases. Furthermore, even if we start in broad resonance with , we need to be sufficiently large, so that does not turn smaller than unity before the backreaction effects from the daughter field(s) are noticed. Taking into account that the scale factor behaves as in this scenario, a transition of broad-to-narrow resonance takes place whenever , i.e. in a time from the start of the simulation. Therefore, we want this time to be larger than the back-reaction time . In practice, we cannot simulate cases for , because for these , and hence we would enter into narrow resonance before backreaction matters. We have simulated cases in the interval . Let us notice that the upper bound on the coupling to prevent radiative corrections, , corresponds to . Of course, in supersymmetric theories, radiative corrections from bosons and fermions tend to cancel each other. In such theories the coupling constant can be in principle much greater than . As in this work we want to be as generic as possible, we will allow ourselves to consider higher couplings. However we will only reach up to , as this corresponds to the largest resonance parameter we are capable of simulating. See Appendix B for an extended discussion about this.

Figure 6: We plot the different times obtained from the lattice simulations of the inflationary model with different resonance parameters. We have joined the points with a straight line, and the orange band corresponds to the values of Eq. (3.19).

3.2.1 Onset of non-linearities, energy evolution and decay time

In Fig. 6 we show the backreaction time , obtained from our lattice simulations. We define again as the moment when the inflaton conformal amplitude starts decreasing abruptly, due to the back-reaction of the excited fields. We show as a function of the resonance parameter . For all simulations, we see that


We do not observe a clear pattern for as a function of , as we saw in the case. This is however expected. The reason is that, in the present case, we cannot differentiate whether a mode is placed in the middle of a resonance band or not. Now each mode experiences a rapid scanning of bands due to the expansion of the Universe. Actually, as described in [10], the resonance in this system is stochastic, precisely due to the scanning over the resonance bands. In [10], it was well appreciated that when solving the Mathieu equation for different modes , for the same initial resonance parameter , the Floquet index oscillates constantly around zero as we go surveying the various modes999The Floquet index alternates between positive and negative values inside a certain envelope curve. The specific form of this envelope is however irrelevant for us now, so we just refer to the interested reader to check Fig. 10 and Eq. (81) of [10].. While for a given mode the Floquet index can be positive , at a neighboring mode it may become negative, , even if . The occurrence of positive and negative ’s is of course not symmetric, but in a proportion 3:1, so that overall there is always a net effect of particle creation [10]. The excitation of a given mode goes receiving alternating positive and negative ’kicks’ in a proportion 3:1, so that in some moments grows, and in others it decreases, but in the overall there is always a net growth. The ’wiggly’ pattern of as a function of is, therefore, just a reflection of the stochastic nature of the resonance in this system. To our knowledge, the pattern depicted in Fig. 6, has never been shown before. Due to the stochastic nature of the resonance, one cannot predict exactly for a specific initial resonance parameter .

Looking at Fig. 6, we appreciate that the onset of the backreaction, and hence the start of the inflaton decay, happens always in a time . Similarly to the analytical calculation presented in Sect. 2 for the quartic case, one can also derive an estimation, based on the linear regime, of the time it takes for an efficient transfer of energy into the daughter field(s), for the quadratic case101010Given the stochastic nature of the resonance, this calculation is perhaps less transparent, but it is expected to capture well, in principle, the order of magnitude value.. As such computation was presented in [10], we just quoted their result (adapted to our notation) in our Eq. (2.19). Taking as a reasonable averaged value of the stochastic Floquet index , then . For , then . As in the quartic case, we see that is a good estimation of the back-reaction time (ignoring of course the stochastic pattern seen in Fig. 6). It is not, however, a good estimation of the decay time of the inflaton, which we estimate next.

We can understand better the post-inflationary dynamics at if we analyze again how the different energy contributions evolve as a function of time. The total energy can be written as a sum of its components as




where and are the kinetic and gradient energy of the fields ( labeling their conformal amplitude), and and are the interaction and potential energies.

Figure 7: Left: We show for the quadratic preheating case and , the evolution of the different energy components of the system as a function of time, see Eq. (3.22). We normalize them to the total energy at initial times, . The gray, red, and blue vertical dashed lines indicate the times , and . Right: We show the times (red circles) and (blue squares) as a function of obtained from lattice simulations.

In Fig. 7 we show the evolution of the energy contributions as a function of time for a particular resonance parameter. We take, as before, the oscillation average of the different functions. One of the most interesting properties of this system is that the equipartition identities


hold for all times. This can be observed in Fig. 7.

Let us begin by noting that, despite the spiky patter of exhibited in Fig. 6, the dominant energy fractions at show however, much less scattering with than in the case of . The energy fractions are mostly independent of the resonance parameter, and are given by

Energy Fractions at :

The errors simply reflect the (random) scattering of energies with . We see again that at , almost all of the energy remains yet in the inflaton. When the tresshold of of energy transferred is surpassed, backreaction then becomes noticeable, and the inflaton amplitude starts decaying. The other energy components, , , remain always at less than levels during , independently of .

We can define again a time scale characterizing the moment when the system enters into a stationary regime. As equipartition holds all the time, we cannot determine now a specific moment when equipartition is verified to better than a certain degree (as we did in the inflationary case). However, we can define at the onset of the stationary regime, understanding the latter now as the regime when the inflaton kinetic and potential energies do not evolve appreciably anymore within one inflaton oscillation period. In practice, we define at the moment when these energies do not change more than within one oscillation. This threshold is not as arbitrary as it seems: at the earlier times , the relative change of the dominant energies not only is bigger than , but also changes in time. However, at times , with defined as just said, the relative change simply remains always below the threshold. Let us note, however, that this does not mean that these energies do not evolve in time at . Actually they evolve smoothly, but the relative change (within an oscillation time scale) is simply very small. Extracting that way from our lattice simulations, we find the data to be very well fitted (see right panel of Fig. 7) by,


Once again, we see that the larger the resonance parameter , the longer it takes the flow of energy from the inflaton to the daughter fields to cease. At this time, the dominant energy components are actually rather independent of the resonance parameter for . Their relative fractions are given by

Dominant Energy Fractions at ():


again with the errors reflecting some scattering of the energies with . The interaction energy is a very sub-dominant component which remains also almost constant after . The inflaton gradient energy and the potential energy density , also sub-dominant components, show however some trend of energy exchange: as increases, grows and decreases. We provide the following estimations based on fits obtained within the range ,

Sub-dominant Energy Fractions at ():


For , we observe that the potential energy becomes marginal, with , while the inflaton gradient energy seems to saturate to a fraction , which still remains subdominant as compared to