Computation of Gravitational Particle Production Using Adiabatic Invariants

Analytic and numerical techniques are presented for computing gravitational production of scalar particles in the limit that the inflaton mass is much larger than the Hubble expansion rate at the end of inflation. These techniques rely upon adiabatic invariants and time modeling of a typical inflaton field which has slow and fast time variation components. A faster computation time for numerical integration is achieved via subtraction of slowly varying components that are ultimately exponentially suppressed. The fast oscillatory remnant results in production of scalar particles with a mass larger than the inflationary Hubble expansion rate through a mechanism analogous to perturbative particle scattering. An improved effective Boltzmann collision equation description of this particle production mechanism is developed. This model allows computation of the spectrum using only adiabatic invariants, avoiding the need to explicitly solve the inflaton equations of motion.

Intuitively, the fast oscillatory component in an otherwise slow Hubble expansion rate H after the quasi-dS era leads to the production of particles. In contrast with preheating related particle production (e.g. [35][36][37][38][39][40]), this scenario requires mediation by the massless gravitational sector, which cannot be integrated out. Also, unlike the scenarios of refs. [27,[41][42][43], the phase space structure of the initial inflaton degrees of freedom is fixed by the Bunch-Davies/adiabatic vacuum prescription during inflation. It is in this sense that the particle production considered in ref. [8] might naively be thought of as part of the gravitational particle production resulting from the transition out of the quasi-dS era.
In ref. [44], an analytic formalism for computing the particle production in this scenario was constructed based on a double expansion in H/m φ 1 and H/m χ 1. In using the formalism, a Gaussian spectral model was introduced to obtain a compact formula. It was based on the assumptions that the Fourier spectrum of the fast time variation of the inflaton field is approximately Gaussian, and that most of the particle production takes place during the first few oscillations. A production rate was then estimated using a naive matching to the Fermi's golden rule, and the Boltzmann equation was integrated to obtain the final number density.
However, there were at least two problems with this. Firstly, the matching to Fermi's golden rule was done in a non-systematic manner, and neglected the long time behavior of the dynamics, although it certainly is correct to O(1) accuracy. Secondly, the high k part of the particle production spectrum requires the fast oscillation component of H to be tracked for longer than a few oscillation periods. This is because the time t * at which the time-integral in the formalism picks up the dominant contribution to the particle production amplitude is when corresponding to the on-shell condition of eq. (1). One of the main results of this paper is to correct these shortcomings. More specifically, using the technique of adiabatic invariants, we construct an accurate time model instead of a spectral model. We will see that these will lead to a 2+O(1) factor correction to the final spectrum compared to that of ref. [44].
It is important to note that the adiabatic invariant formalism of this paper is distinct from the adiabatic vacuum of refs. [45][46][47]. In particular, the latter is essentially the WKB approximation, and applies only to linear wave equations, while the former requires oscillatory motion, and can be applied to non-linear Hamiltonian dynamics. Furthermore, in the context of gravitational particle production, the WKB approach is used to solve the χ k mode equation, and is controlled by an expansion in H/m χ . In contrast, the adiabatic invariant formalism of this paper is used to analyze the dynamics of the inflaton φ and the scale factor a, and is controlled by H/m φ . This then determines the χ mode dispersion, which controls its particle production.
A second main result of this paper demonstrates that one can use the formalism of ref. [44] to compute the spectrum numerically more efficiently than the brute force methods used in the literature. We note that what is expensive to compute numerically is the slow time-varying component weighted integration over fast oscillating functions which converge slowly to a negligible contribution. 3 Instead, the formalism of ref. [44] already isolates the fast time variation such that its numerical integration can be more than O(1000) faster in convergence. The relative gain in efficiency depends on the amplitude of the spectrum.
It is important to note that although the original motivation for this scenario is the production of dark matter, the computational tools presented in this work are useful for constraining top down constructions as well. For example, in contrast with moduli problem scenarios such as those of refs. [48,49], moduli much heavier than the inflationary expansion rate can be gravitationally produced as well during the first few oscillations out of inflation [50] (where the largest number is produced in the first oscillation). Furthermore, even if the particles decay too rapidly to be dark matter, they can have observable consequences as long as they are sufficiently long lived, as discussed for example in refs. [2,5,[51][52][53][54][55][56][57][58][59][60]. 4 The order of presentation will be as follows. In section II, we summarize the previous work of separating the H time scale from the m φ time scale and explain its shortcomings [44]. Here, we also establish the notation for this paper. In section III, we explain how the adiabatic invariant can be constructed and how to use it to compute the long time behavior of slowly varying fields. This tool will play a crucial role in section IV, where we present a more accurate approximation of the Boltzmann equation, which is one of the main goals of this paper. We parameterize the associated errors based on a coarse graining over short time scales. In section V, we introduce an explicit time model of the inflaton dynamics that separates the short and long time scale dependences.
In section VI, we derive an analytic formula for the coarse grained particle production amplitude that includes {1 → 2, 2 → 2, 3 → 2} resonance contributions. The formulae of this section are analogous to the S-matrix amplitude, which will need to be supplemented with an analog of phase space integrals to express a physically measurable number density. The coarse graining time scale is computed in section VII by minimizing the competing errors associated with quantum interference and general sum to integral approximations. We calculate the spectrum and number density estimates in section VIII by utilizing the explicit solutions obtained from the adiabatic invariant tool. In section IX, we present our second main result: numerical integration advantage of the "fast only" formula of eq. (162) over the brute force Bogoliubov formulation of eq. (8). In that section, we also generalize the Gaussian spectral model of ref. [44] by incorporating the adiabatic invariant formalism, and compare it with other computational methods.
We summarize our results in section X. The appendices include a review of the adiabatic invariant formalism used in this paper and a derivation of an explicit formula for the time dependence of the slowly varying Hubble expansion rate after the end of inflation for a generic inflationary potential.

II. ASPECTS OF PREVIOUS WORK ON THE TOPIC
Here we review our conventions and previous work on this topic. Consider the action in the metric signature (1,-1,-1,-1): where M P = 1/ √ 8πG is the reduced Planck mass, R is the Ricci scalar, φ denotes the inflaton field with potential V (φ), and χ denotes a real scalar field. We assume that χ interacts only through gravity controlled by a non-minimal coupling ξ. Pure Einstein gravity (i.e. minimal coupling) corresponds to ξ = 0 while conformal coupling corresponds to ξ = 1/6. The standard background cosmological metric is written as where the scalar factor a is a function of either coordinate time t or conformal time τ . The Hubble parameter H and Ricci scalar R have the conventions whereȧ = ∂ t a and a = ∂ τ a. The energy-momentum tensor is assumed to arise primarily from a homogeneous inflaton field φ, and therefore Einstein's equation can be written as It is well known that the number density n χ can be written as where the Bogoliubov coefficient expression for gravitational production of χ particles with mass m χ is with ω 2 k ≡ k 2 + a 2 m 2 χ + 1 6 (1 − 6ξ)a 2 R as long as |β k | 1 [38,62]. 5 After separating out the time scale much faster than H −1 , the expression can be written as where t end is the time when inflation ends (when the slow-roll parameter is unity), with and the subscript UV denotes the quantities with the low frequency components filtered out [44]. 6 In the work of ref. [44], a Gaussian Fourier spectrum function with a high frequency peak was fit to the properties of the potential to effectively define these UV components. Although this did allow one to analytically compute the dominant contribution to the particle production in a clean formula, there were some shortcomings: 1. Although the largest amplitude of inflaton oscillations occurs at t end , particle production remains significant at later times when the amplitude is appreciably smaller. Indeed, the final number density was computed in that reference by converting the first few oscillations computation into a rate using (where σ is related to the Gaussian spectral model), setting t end → t, and using a Boltzmann equation to compute the production rate. However, no formalism was provided on how to solve for the long time scale dependence such that the replacement t end → t can be made for an inflaton potential that is not purely quadratic. One expects this long time scale to be particularly important for the k modes where since for these cases, β k settles to its asymptotic final value only for t t end . 5 We will restrict ξ to values such that ω 2 k > 0 for all times. 6 More explicitly, low frequency components refer to time dependent quantities that vary on a time scale of order H −1 end or longer. Also, by filtered out, we mean subtracted out from the total.
2. The choice of the coarse graining time scale identification (∆t) −1 = σ/ √ 2π in equation 6.11 of ref. [44] was an estimate for the coarse graining time scale based on the naive expectation of the Gaussian frequency model: However, since the actual Gaussian model comes from the square ofφ 2 whereas eq. (16) came from the φ spectral representation itself (see equation 5.23 of ref. [44]), there is an ambiguity in this identification.
Both of these shortcomings will be addressed in this paper through a novel time dependent model based on an adiabatic invariant construction replacing the Gaussian Fourier spectral model of ref. [44]. The use of the adiabatic invariant also makes more precise exactly whichφ 2 component is being kept when evaluating β k contributions from "fast varying" (UV) modes.

III. CONSTRUCTING THE ADIABATIC INVARIANT
One of the issues outlined near eq. (14) is the long time scale behavior of the inflaton oscillation amplitude necessary for the accurate determination of the amplitude for the high k modes. To compute this, we use the adiabatic invariant technique reviewed in appendix A. This technique allows one to construct an approximate invariant Q because of the combination of the properties that Hamiltonian equations of motion preserving phase space and energy being conserved in the limit that an external time-translation breaking function becomes a constant.
To use the adiabatic invariant to isolate the slowly varying part from the fast varying part of T µν cleanly, we decompose the FLRW scale factor as a(t) = a slow (t) + a fast (t). Let the slowly varying source function for the adiabatic invariant construction be λ(t) = a slow (t) and let the coordinate of the phase space q (in the notation of the appendix A) be the inflaton field φ. An interesting part of the procedure we are going to follow is that a slow (t) does not need to be specified to obtain the adiabatic invariant.
We start by carrying out the integral in eq. (A25). We choose C to be the φ path that begins and ends at φ(v). We can set the initial condition φ(v) at which π(v) = 0. This can typically be done since the motion is oscillatory by assumption and the turning point typically corresponds to π = 0. We can therefore write E(v) from eq. (A4) as where H is the Hamiltonian and φ is any field value in the path C. This results in an adiabatic invariant equation for φ(v): Since this expression is a constant, we have a prediction for how φ(v) varies as a function of a slow (v).
Let us now obtain a more explicit expression for φ(v) using the conservation of energy per integration cycle¸C dφ. First, canonically normalize the fields such that the Lagrangian is a 3 If there are non-minimal gravitational couplings, such terms can be absorbed into V . The momentum conjugate to φ is π = a 3 slowφ , and thus the Hamiltonian is H(φ, π, a slow ) = 1 2 Solving for Π using eqs. (17) and (19), the adiabatic invariant can now be written as where φ C,− (v) and φ C,+ (v) are respectively the smallest and largest values φ takes during an oscillation cycle beginning at time v. The turning points are determined by where φ C,− (v) can then be viewed as a function of φ C,+ (v).
The mass dimension of Q is 3 which is the same as that of the phase space density of canonically normalized scalar fields. This adiabatic invariant is the phase space density multiplied by the time-dependent volume scaling factor a 3 slow (v), which is intuitively the reason why it is constant. The interesting nontrivial aspect is that the phase space density here is not defined with respect to free particles but with interactions turned on. Starting with an initial condition φ C,+ (v i ), one can replace Q with this initial data and thereby solve for the slowly varying φ C,+ (v) usinĝ where v i is the initial time (which we will often take to be the time at the end of inflation in practice). This can also be rewritten as where For some choice of potentials V (φ), I(φ, φ C,+ (v)) will be simple enough to compute such that φ + (v) will be easy to obtain from eq. (23). We will give some examples below. An alternative method is to compute V m ≡ V (φ C,± ) as a function of time using the formalism of appendix B, and solve for φ C,± (v) through where the coefficients are defined in appendix B, and V m is the potential value at the peak of the oscillations. Let's turn to some examples.
A. Quadratic potentials In this notation, we know for the quadratic potential V (φ) = m 2 φ 2 /2 such that eq. (24) becomes and using eq. (21), we find that which gives This is one of the most important generic examples since most potentials have the quadratic mass terms dominating as the amplitude of the oscillations decreasing. That is why we separated this example from the next set of examples. By taking A = 0 in eq. (B51) of appendix B, one can see the time scale of a slow /ȧ slow is determined entirely by H −1 end , which is much larger than the oscillatory period given by 2π/m φ . This is because the solution of eq. (20) is Q = 6πM 2 P a 3 end H 2 end /m φ for a purely quadratic potential.

B. 2n power potentials
In the case of more general even powered potentials, for some integer n ≥ 1 we have with eq. (21) giving φ C,− [φ C,+ ] = −φ C,+ . Now eq. (24) evaluates to where 2 F 1 is the hypergeometric function. Note that which allows eq. (23) to be evaluated as Consider n = 2, which corresponds to a quartic potential. The field φ C,+ (v) scales consistently with the conformal scaling dimension of the scalar field and significantly different from the quadratic power index −3/2 (see eq. (29)). That is because even when the oscillation amplitude decreases, there is no mass term for such potentials.

C. Trigonometric integrability
With the potential eq. (21) gives φ C,− [φ C,+ ] = −φ C,+ . We also find that as |φ| → ∞ the potential is bounded from below. Therefore, eq. (24) evaluates to where E(φ, m) is the elliptic integral of the second kind. Note that which allows eq. (23) to be evaluated as expressing a nontrivial scale factor dependence of φ C,+ (v). One can easily plot the function to see that φ C+ (v) behaves very close to a −3/2 . That is simply because as a slow (v) increases, the field samples the minimum of the potential which becomes quadratic. This example also illustrates that in general, one cannot find an expression for φ C,+ (v) as a function of a slow (v) using elementary functions.

D. Asymmetric potentials
Consider the potential which has asymmetric turning points due to the cubic interaction term. Even though this potential is not bounded from below, as long as one considers the initial oscillation amplitude φ C,+ (v i ) close to the minimum point of the potential φ min = 0, there will be stable classical oscillations. Because of |φ C,− [φ C,+ ]| = |φ C,+ | in the turning points, this is an example that is not covered by techniques such as that presented by ref. [63]. Although eq. (24) can be evaluated giving a result similar to eq. (37), it is more instructive to use eq. (26). Using eqs. (B14), (B32), and (B49), we find the long term field time dependence of where the physical, initial condition dependent charge is as given by eq. (B28). 7 Clearly, eq. (40) shows that the asymmetric time dependence of φ C,± is suppressed parametrically by which is intuitively expected since this is a measure of the cubic potential strength to the quadratic potential strengths: i.e. the deviation ∆φ from a symmetric excursion at a linearization level. It is interesting to see explicitly in a simple manner that field interactions lead to a non-scaling behavior of the energy density, in contrast with other interaction results such as eq. (33) (or equivalently eq. (B7)).

IV. APPROXIMATION OF THE BOLTZMANN EQUATION
The Boltzmann collision equation allows number density computation of a particle production process. The collision term includes both a production rate and an annihilation rate. As noted in eq. (13), an effective production rate was estimated in ref. [44], but as explained in section II, the effective production rate can be improved by better time domain computation of the inflaton field.
We present here an approximation of the quantum production rate induced by Bogoliubov transform, which models the χ particle number density spectrum f χ (k, t) = |β k (t)| 2 using a non-negative effective production rate spectrum. This is nontrivial as the actual time evolution of f χ (k, t) is not a monotonically growing function unlike a production rate of a Boltzmann equation. 8 However, when coarse-grained over a short quantum fluctuation time scale, 9 the production rate is positive definite. Ultimately, this is because the initial state is a vacuum, which by construction has no particles, and the final state at time t cannot have negative number of particles. The main goal of this section is to set up the parameterization of such a positive definite quantity approximation and compute the error incurred by neglecting the fluctuations and coarse-graining.
Given the short time scale separation eq. (9) of the Bogoliubov coefficient, the Boltzmann collision equation will be approximated as where we definedb and ∆t is a characteristic coarse graining time that encompasses many oscillations and will be chosen such that the error of this approximation is minimized. After accounting for all errors, an expression for ∆t is given by eq. (119). The approximation eq. (45) of the Boltzmann rate equation is equivalent to the spectrum density being modeled as where the errors are estimated by with v m in a discrete set of times whose Let us now justify this result. First, to make uniform the kinematic description of the fast time dynamics, we partition the real time variable t into a pair of variables {v m , s} for m ∈ X ≡ {0, 1, 2, ...}: where and where the finite time interval ∆t(v m ) will be fixed through an approximation scheme that minimizes the error terms of eq. (47). The estimation of the phase interference error eq. (50) will be deferred to section VII due to the prerequisite results from sections V and VI. It will be shown below that the other errors associated with eqs. (48), (49), (51), and (52) have upper bounds that approximately scale with positive powers of H(v)∆t(v), which is assumed to be a small expansion parameter. This subsection will justifying neglecting the term eq. (51) by showing that the relative sum to integral error, defined as is suppressed by H(v)∆t(v). With the Euler-Maclaurin formula, the sum can be approximated with an integral as where P 1 (x) is the first Bernoulli polynomial and we define its argument as Slowly varying functions such asb k (v) and ∆t(v) derive their time dependence from a slow (t), and therefore the derivative of such functions are typically on the order of the function multiplied by the Hubble parameter. Thus the error terms of the right hand side of eq. (59) are estimated as respectively to leading order in H(v m )∆t(v m ). Therefore the relative error is of the form with an upper bound given by where c m are order one numbers. This is suppressed if H(v m )∆t(v m ) 1.

C. Remainder errors
We will now estimate the relative remainder errors of eq. (47). These errors will vanishes as t → ∞. For finite time, they are suppressed by positive powers of H(v Nt )∆t(v Nt ).
The relative error associated with eq. (48) is written as Using the inequality we obtain an approximate upper bound of As will be shown in section VI, the fast varying component of the Bogoliubov integrand can be estimated as and therefore an approximate upper bound on this remainder error is The relative error arising from eq. (52) is given by The upper bound of eq. (72) of the Bogoliubov integrand yields and therefore an approximate upper bound on this remainder error is Eqs. (73) and (76) are both suppressed by positive powers of H/m φ and H∆t.

V. A TIME MODEL OF INFLATON DYNAMICS
Although Einstein equations give directly the Hubble rate H in eq. (12), it is useful to compute its time derivative since that will tend to suppress the slow frequency components such that the fast frequency components become manifest. The time derivative of H giveṡ Similarly, the Einstein equation for the Ricci scalar R gives which shows up in eq. (12). To extract the high frequency time variation of this term, it is also useful to take a derivativeṘ where one notes that the first term obtained a large contribution from the potential term. In this form, the second term in eq. (78) is subdominant to the first term because H is much smaller than m φ , which is the inverse time scale of the time derivative. Remarkably we have reduced the determination of the Bogoliubov integrand B k (t) to modelingφ 2 . In this section, we will develop a model ofφ using an adiabatic invariant and a fast-slow time decomposition. Since φ C,± are left-right bounds on the inflaton field φ(t), we can parameterize the field time dependence in terms of a phase Ξ φ (v, s) such that whereφ and the trajectory is only over a time period s ∈ [0, ∆t] . To finish defining the approximate model, we must give a map between t and (v, s). As we will see, we will define this map through v only having ∆t resolution such that the set of points {v ∈ (t 1 , t 2 ), s ∈ [0, ∆t]} approximately covers the same domain as t ∈ (t 1 , t 2 ). This equivalence class definition of v is possible because all the quantities with v dependence will be slowly varying such that they are approximately constant in the time interval ∆t.
We define φ C,± (v) to be the amplitude derived from the adiabatic invariant. This automatically gives v a resolution of the period of the adiabatic invariant. In this case, we see Ξ φ is real since φ C,± are bounds of periodic motion in the adiabatic invariant approximation. As a definition of the model, we require . Now, we attribute the fast time behavior of φ(t) ≈ φ(v, s) in eq. (79) through approximate energy conservation. We can express ∂ s φ as where can be obtained using either φ C,± (v). Eqs. (83) and (84) lead to a quadrature integral that determines the phase Ξ φ by In this paper, we will consider only potentials where there is a single dynamically relevant minimum φ min at the end of inflation. If we keep only the leading asymmetric term of the potential, we have the expansion and if we define this leads to the expansion which makes manifest how α 3 controls the asymmetric nature of the maximum and minimum value of the field excursion represented by φ C,± . Putting this into eq. (86) gives where Using eqs. (93) and (92), the time dependence of [∂ s φ] 2 is approximated (to leading order in α 3 ) as Note the squaring produced a slowly varying term even though ∆φ(v) cos Ξ φ naively looked to be fast varying only. The fast varying part is obtained by throwing out the terms that depend on v only. We will call this UV:

VI. AMPLITUDE OF THE PRODUCTION RATE
In the time interval of s ∈ [0, ∆t(v)], the relevant quantities for evaluation of the Bogoliubov coefficient integrand are eqs. (12), (77), and (78). To leading order in H/m φ , we therefore have which is essentially determined by the inflaton fast velocity squared function [∂ s φ(v, s)] 2 UV . For convenience of interpretation later, this can be evaluated using eq. (96) rewritten in a complexified form as where we substituted V m = 3M 2 P H 2 slow . The key quantity eq. ( (53)) can therefore be written as where we defined and the order of terms in the brackets in eq. (99) is in the order of typical importance in magnitude.
We now arrive at one of our main analytic results of this paper. Using the Boltzmann equation approximation of section IV, we write the spectrum as For brevity of notation F (E k ) has its v dependences that exist through The fact that odd → even processes are mediated by the m φ α 3 (φ − φ min ) 3 ≡ m φ α 3 δφ 3 interaction vertex is clear since without the cubic vertex, the symmetry δφ → −δφ would forbid all odd δφ number changing processes.
To obtain an intuition for a 3 consider the example potential of The inflaton mass will be m φ = √ 2nM 2 /v φ . We choose v φ = 0.5M P . The leading asymmetry term will be α 3 = n − 1 and thus which is an O(1) number. At later times, a 3 (t) will decay as H slow (t) and therefore become a small expansion parameter at large resonance times.

VII. ESTIMATION OF COARSE GRAINING TIME ∆t
Until now, we still have not specified the time width ∆t(v). Its specification dominantly affects the computation in several ways. First, it fixes the errors of eq. (47). Second, it enters in the approximation eq. (79) and the trajectory interval s ∈ [0, ∆t(v)] specifying the adiabatic invariant based time model. As far as most of the errors of eq. (47) are concerned, the smallest error is obtained when ∆t is the smallest. On the other hand, there is a preferred ∆t > 0 that makes the quantum interference error of eq. (50) negligible. The adiabatic invariant interpretation also sets a minimum ∆t. We will see that they can be made commensurate. We will find below that the quantum interference error considerations lead to a time width at an intermediate scale between those set by the inflaton mass and Hubble expansion rate, i.e. m −1 . This scale is larger than the one set by the adiabatic invariant minimum of a single oscillation.
We will first consider the quantum interference error of eq. (50). The dominant interference term is where v p is the nearest time to resonancev k , and therefore dominates the denominator sum. Neglecting the cubic coupling a 3 , and for large frequencies E k ∼ ω * m χ , eqs. (105) and (106) lead to where sinc(x) = sin x x , and A 2 is given by eq. (100). As such, we have for times v p nearest the central peak ofb k . The largest off-diagonal contributions will be where we made use of 1 to neglect sub-leading terms. Therefore, the quantity is an estimate of the relative error.
To minimize this error, we require the arguments of the sine-cardinal functions in the numerator to lie outside the central peak: where the other sign equation is automatically satisfied if this condition is satisfied. The above inequality will be satisfied if for all v p that satisfy eq. (113). The coarse graining time ∆t(v p ) is therefore set to the lower bound: where we have assumed E k ≈ m φ . For general times v away from the peak, we take H(v p ) → H(v) to ensure H(v)∆t(v) remains a small expansion parameter. As such, we set as the coarse graining time width of the Boltzmann rate approximation. This time interval is longer than the single period of the adiabatic invariant construction oscillations which is of order ω −1 * ∼ m −1 φ since H(v p ) m φ for these scenarios. However, we see that H(v)∆t(v) → 0 as v → ∞, which allows for an accurate fast-slow decomposition. A slow varying function can be treated as approximately constant over the time interval (v, v + ∆t(v)) because the Taylor expansion terms are suppressed by positive powers of H(v)∆t(v).

VIII. PREDICTED SPECTRUM AND NUMBER DENSITY
Using the time model described by eqs. (104), (105) and (106), the particle density spectrum 10 f χ (k, t) and the number density n χ (t) can be computed using slow varying quantities derived from the adiabatic invariant formalism. We begin with the integration approximation with F (E k ) given in eq. (106). Here F (−E k ) and the cubic vertex amplitudes eqs. (102) and (101) were treated as subdominant to F (E k ) and A 2 of eq. (100) respectively. It is evident that the particle production peaks at the resonance timev k when the energy matching condition E k (v k ) = ω * (v k ) is satisfied. As the inverse time width is much smaller than the energy scale set by the inflaton mass, one can treat the limit ∆t −1 → 0 as a good approximation. The spectrum of the particle production rate then follows Fermi's golden rule: with A 2 /2 acting as the transition matrix element, and the delta function as the energy density of states. In analogy to refs. [41,42], this can be interpreted as the particle scattering process φφ → χχ of two inflatons with total energy 2ω * (v) under going gravitational annihilation to produce two χ particles with total energy 2E k (v) in the center of mass frame. After integration of eq. (122), the spectrum and number density limit to respectively, wherev k is the time such that The number density was obtained by resolving the delta function through wave-vector k integration. Integration of the spectrum eq. (123) with bounds k ∈ k (v k ),k(t) yields an alternative formula for the number density. Using the explicit form of eq. (100) gives where we treated and ω * = m φ respectively, and this yields the estimates where H end = H slow (t end ). The number density scales as n χ ∝ a −3 (t) at large times, as expected for nonrelativistic massive particles. Note that the spectrum scales as f χ ∝ k −9/2 , matching the results of ref. [8]. These estimates of the number density are larger by a factor of 2 from the results of the spectral model of ref. [44]. The spectral model suffered from an ambiguity of setting the correct time scale of Fermi's golden rule. The time model has no such ambiguity, with eq. (122) emerging as a result without being put in by hand.
For a more general potential which is quadratic near its minimum, we show in appendix B that a perturbative expansion in the adiabatic invariant charge Q allows us to write the Hubble expansion rate as where we defineQ and are the relevant coefficients. The leading u −3/2 term of eq. (131) quickly dominates at later times. The spectrum eq. (126) can therefore be estimated as for large k. Let us now neglect the time dependence of the frequency ω * as H 2 /m 2 φ relative corrections to m φ . The integral for the number density in eq. (127) is then proportional tô and therefore the number density is estimated as at late times for m φ > m χ . Note the differences between this number density and eq. (130) for a purely quadratic potential. The effect of the non-quadratic terms of the inflaton potential is summarized by the factors F andH end /H end , which are derived from the adiabatic invariant Q and the α n coefficients. In contrast with ref. [63], this formalism is applicable for potentials that are asymmetric about their minimum.
As an example, consider the potential where the large numbers stem partly from the large M p /m φ ratio (see e.g. eq. (B49) and the comment below (B51) for further explanation). When the inflaton potential is quadratic, the Hubble rate has the time behavior of eq. (128). For general potentials, the Hubble rate is more accurately modeled as eq. (131), with the leading order u −3/2 term dominating at late times. The ratio of these two estimates will deviate from unity for non-quadratic potentials, and this will change the magnitude of predicted spectrum and number density. For the potential eq. (141), this ratio isH Overall, the predictions given by eqs. (130) and (139) differ by a factor of around 1.6 for this example. As shown in figure 1, numerical results demonstrate this for the example potential eq. (141), with eq. (136) being a better numerical fit of the spectrum than eq. (129) at large k.

IX. COMPARISON OF DIFFERENT COMPUTATIONAL APPROACHES
In this section, we compare some different approaches to computing the super-Hubble mass scale particle production. The second main result of this paper presented in section IX B is to note the O(m φ /H) (which can easily be a factor of 100) enhancement in numerical integration efficiency is achieved for a formulation obtained by subtracting out the slowly varying component of H. All of our numerical illustrations will be done with a prototype inflationary model presented in subsection IX A. In subsection IX C, we compare the particle production results of the Gaussian spectral model of ref. [44] with a generalized Gaussian spectral model based on an adiabatic invariant and an analogous time-model of V. In subsection IX D, we illustrate the time model in the context of potentials where the maximum displacement from the minimum of the potential is asymmetric on two sides of the minimum.

A. Inflationary model used for illustration
As a prototype of an inflaton potential for low-scale inflation, we test the time model developed in this paper by considering the hilltop model : where n > 2 is an integer and v φ > 0 is the vacuum expectation value of the inflaton at the minimum of its potential [31,[64][65][66][67]. The mass of the inflaton near the minimum is and there is an effective cubic interaction of field displacements from the minimum (as well as higher order interactions). The Hubble expansion rate at the end of inflation is estimated as H end √ 3M P −1 M 2 . By calculating the large-scale curvature perturbations, the scalar spectral index and the tensor-to-scalar ratio is found to be where the e-folding number N of the CMB lies between 50 and 60 [68]. As observed by the Planck satellite, the overall normalization of the curvature perturbation implies which relates M and v φ . Hence there is one free parameter left, which can be taken as v φ [64]. The measured range of the spectral index (n s = 0.968(6) at 1σ level [69]) can be made consistent with n ≥ 6 [64]. We take n = 6 in the potential eq. (148) and v φ = 0.5M P to compare to the numerical results of ref. [8]. This is a somewhat tuned model that could be destabilized by loop generated Planck suppressed operators, but it serves as an algebraically simple demonstration of the adiabatic invariant formalism. 11 It is interesting that the cosmological data driven phenomenology favoring n ≥ 6 also enforces m φ /H 1 since which is equivalent to the condition m φ ≈ 30H end with n = 6. As we saw in eqs. (143) and (144), the expansion in non-quadratic parameters α n ∝ ∂ n φ V defined in eq. (B13) allows us to capture many different models [15,70,71] even though we focused for numerical illustrations on the model of eq. (148).

B. Exact versus fast component numerical integration
In this section, we will summarize our numerical procedure to compute β k using the brute force exact integration eq. (8) and the fast only component integration eq. (9). We will then discuss figures that illustrate the differences between these methods, and explain observed features.
Both methods first required solving for the solution φ ode (t) to the inflaton equation of motion, given by the non-linear ordinary differential equation where we assumed a background homogeneous inflaton field that dominates the energy density. Here the initial conditions and parameters were chosen such that at least N ≈ 55 e-folds in the scale factor occurred before the end of inflation. We used the initial conditions φ(t i ) = 0.17φ end andφ(t i ) = 0, and found that t i − t end ≈ −63H −1 end . 12 It was found that the relevant portions of the ODE solution were nearly independent of the slow roll initial conditions. The time t end when the quasi-dS era ends is given by the solution to following the criterion used in ref. [8]. For the sake of convenience, we shifted time such that t end = 0. We also set our time scale such that M = 5297H end and M P = 1.565 × 10 7 H end , where we have defined which is approximately the usual Hubble expansion rate at the end of inflation. The Hubble rate and Ricci scalar were determined by the Einstein equations 3M 2 P H 2 ode = 1 2φ 2 ode + V (φ ode ), and M 2 P R ode =φ 2 ode − 4V (φ ode ), respectively. The scalar factor was then computed by solvingȧ ode = H ode a ode . With these quantities in hand, we evaluated the exact β k integration as where t i was chosen sufficiently in the past compared to t end to approximate the adiabatic initial conditions of the Bunch-Davies vacuum. In our computations, we chose H end (t i − t end ) = −5.
To obtain the necessary components for the fast only integration, we solved the adiabatic invariant equation to obtain the slow time behavior of the Hubble expansion parameter. We first defined and computed it for various values of V m starting with V (φ end ) and proceeding to smaller values as needed.
The step function ensured that finding the turning points was numerically unnecessary, and the smooth monotonic nature of the integration ensured that only a few sampling points of V m were needed to obtain a good interpolation. The monotonic nature of J(V m ) also guarantees the existence of its inverse. The adiabatic invariant equation Q = a 3 slow J(V m ) can then be written as which determines V m in terms of a slow . We computed the adiabatic invariant using the initial conditions as Q = a 3 end J (V (φ end )). The time dependence of a slow was then obtained by integratinġ with the initial condition being a slow (t end ) = a end . For convenience, we chose for our numerical work the normalization a end = 1.
With a slow (t) in hand, we arrived at the desired quantities where H fast is defined to be zero before t end . One can compute other slowly varying quantities by taking a slow (t)). For example, ρ slow = V m and R slow = −V m /M 2 P . However, only H slow was necessary for our purposes here asṘ ode is dominated by its oscillatory components.  156) and (162), respectively, are compared for the inflationary model of eq. (148). Although both methods tend towards nearly identical asymptotic values, the fast integration converges 3 orders of magnitude faster in integration time than the exact integration. The long convergence time of the latter is due to the slowly varying component in eq. (8), which damps down only as t −1 for large times. On the right, we see that the exact integration is sensitive to the initial conditions before the end of inflation, while the fast integration is not. Starting the exact integration at t end would lead to a large deviation in the asymptotic limit, which is due to the adiabatic condition not being satisfied at that time.
The fast only integration refers to subtracting out the slow components and sub-leading contributions, such as assuming thatṘ HR and m 2 χ R, as indicated by eq. (9). We also assumedṘ slow Ṙ fast and kept R ode for the sake of computational simplicity and accuracy. For similar reasons, we avoided replacing a ode with a slow . In summary, we evaluated where E k (t) = k 2 /a 2 ode (t) + m 2 χ . The integration time t f should be chosen long enough such that the dominant resonance occurs for a given k mode i.e. t f >v k , wherev k is defined in eq. (125). This is particularly important for particle mass m χ close to the m φ threshold. The scale factor at resonance is a(v k ) ≈ k/ m 2 φ − m 2 χ , and thereforev k → ∞ as m χ → m φ . In practice, one uses a cutoff k c for the wave-vector integration of the number density. The relative error incurred in the number density can be estimated to be as implied by eq. (130). Figure 2 illustrates the efficiency advantage of fast integration formulation of eq. (162) over the exact integration formulation of eq. (8). For this illustration, the parameters were chosen such that the 2 → 2 resonance at time (see eq. (122)) in the m φ /H end ≈ 30 model of section IX A is at t satisfying which explains the spike in the dashed curve at near unity time. The fast integration formulation of eq. (162) occurs because the integration over the sinc function in eq. (121) converges on a time scale as given by eq. (119). The comparatively slow convergence of the "exact" numerical integration stems from which has a fluctuation amplitude ∆β k scaling as t −1 if k/a 1 after the resonance is reached. This requires an integration time t of for |∆β| is the accuracy of the amplitude that is desired. For figure 2, the final |β k | is O(10 −4 ), which means we require an accuracy of |∆β| ∼ |β k | ∼ 10 −4 . This gives a convergence time of for the brute force numerical integration ("exact integration"). As shown by eq. (167), the advantage of the fast only integration method increases for small |β k | which is typical for the scenarios of physical interest. Figure 2 illustrates another advantage of the fast only integration technique over the exact integration for a mode slightly off resonance at the initial time t end . The exact integration is more sensitive to satisfying the adiabatic vacuum condition during the quasi-dS era because imposing the standard adiabatic vacuum condition at a non-adiabatic time of t end is apparently an excited state which contains a sizable amount late time high frequency particle modes. On the other hand the fast only integration started at t end has by construction subtracted out the leading nonadiabaticity governed by O(H end ) dynamical scales, giving rise to an apparently acceptable level of consistency with the implicitly assumed adiabaticity (with respect to high frequency modes) at t end .
To put this another way, assuming an adiabatic boundary condition at a nonadiabatic time with an adiabaticity violating scale of O(H end ) is equivalent to assuming the extra presence of late-time high frequency (i.e. ω k /a end H end ) modes, even though naively O(H end ) nonadiabaticity should not contain any high frequency components. This somewhat surprising result may be due to the fact that imposing the adiabatic vacuum boundary condition at a nonadiabatic time t end implicitly contains a step function in time with a non-negligible amplitude, which necessarily includes non-negligible high frequency components.
The produced χ number density as a function of m χ is shown using various approximations in figure 3. The "Exact integration" and "Fast only integration" use eqs. (156) and (162), respectively, along with eq. (7). The figure clearly shows that the number density reaches a significant threshold at m χ = m φ . This is due to the δφδφ → χχ resonance becoming kinematically suppressed. A similar suppression of the δφ → χχ resonance causes the slope change around the m χ = 1 2 m φ threshold. The largest deviation between the fast only and exact integrations occurs at low m χ , where the assumption of m 2 χ R made in eq, (10) becomes less accurate.   171), and the time model eq. (175). Note that for mχ < m φ , and k > m φ a end , the time model and exact integration match while the spectral model differs from both by a factor of 2. The spectrum in this regime scales like k −9/2 , which appears as linear behavior on this log-log plot.

C. Time model versus Gaussian model
In this subsection, we generalize the Gaussian spectral model of ref. [44] and compare its computational accuracy with that of the symmetric time model of eq. (121). In both models, we compute the spectrum using the Boltzmann equation approximation form of eq. (47): where ∆t(v) andb k (v) will be defined differently based on the model: i.e. {∆t → ∆t (GS) ,b k →b As discussed in section II, the Gaussian model had several shortcomings partly because of the ambiguities associated with mapping to Boltzmann equations as well as not having a long time understanding of oscillating field amplitudes. The Gaussian spectral model with the adiabatic invariant evolved fields resolving these deficiencies is defined by generalizing the Gaussian model of ref. [44] as where the width σ(v) is defined as with as a small positive parameter to ensure computational convergence. This model still has the limitation that the spectrum has a limited shape, but it will be able to reproduce the high k part of the spectrum as can be seen in figure 4. This analytic fit to the high k part of the spectrum was used as the primary guidance in the generalization of the Gaussian model.
In the time model, we have instead from eq. (121) the expressions which are similar to eq. (171) in having a peaked function, but different in how fast that peaked function falls off with E k − ω * . The sinc function falls off much slower and has multiple peaks. Furthermore, the prefactors dependent on the non-minimal coupling parameter ξ are different and can become significant for large ω * /E k ∼ m φ /m χ applicable for small k/(a end m φ ). 13 These features can be seen in figure 4. Shown in figure 3 is also the integrated number density of the "Time model" (eq. (175)) and the "Gaussian spectral model" (eq. (171)). The "Time model" plot does not fall off with large m χ /m φ because eq. (175) grows as m 2 χ and therefore the rate is boosted by an m 4 χ factor. Even though the minimally coupled case plot seems to indicate that the large m χ /m φ region has a vanishing asymptote, this is an illusion associated with the vertical resolution of the figure. This mismatch between the "Time model" approximation and the better approximation of the "Fast only integration" is due to the loss of validity in the kinematic region in which eq. (2) does not have a solution (e.g. see eq.(118)).
In fact, for the symmetric time model, one can interpret the nonzero a 3 n χ for the m χ > m φ parametric region as an error on the particle production computation for any finite v. Using eqs. (7), (47) and (175), we can estimate this error as which with which is suppressed by H i /m φ and the small prefactor. The Gaussian model in contrast has an exponentially suppressed error as m χ /m φ becomes larger than unity.

D. Asymmetric time model
The time model in subsection IX C did not account for the δφδφδφ → χχ and δφ → χχ scattering. Using eq. (104), we can make the substitution {∆t → ∆t ASTM ,b k →b ASTM k } in eq. (169) where and A n are amplitudes given by eqs. (102), (100), and (101). Notice that because ω * ≈ m φ while E k redshifts, one can see from the denominators of eq. (181) that δφδφδφ → χχ amplitude A 3 will be reached first before the δφδφ → χχ amplitude as a function of time for 2m χ < 3m φ The time evolution of the particle spectrum for an example parametric point is shown in figure 5. In the following argument, we will set t end = 0 for convenience. Using the resonance condition 3ω * = 2E k , we would predict a step in particle spectrum time evolution at which is consistent with the step feature in the figure near H end t ≈ 3. The step at H end t ≈ 5.5 is due to the dominant 2ω * = 2E k resonance.

X. SUMMARY
In this paper, we developed analytic and numerical techniques to compute the Bogoliubov coefficient β k more accurately and efficiently. Using the adiabatic invariant equation, we computed the slowly varying components of the β k integrand necessary to implement the fast-slow decomposition formalism of ref. [44]. In the process, a numerical integration technique was created that is O(1000) faster than a brute force computation.
To compute the particle production rate and number density, an approximation of the Boltzmann equation was presented, and a coarse graining time ∆t was introduced. We demonstrated that ∆t minimizes the error of this approximation at an intermediate scale between the inflaton mass and Hubble rate, and that the overall error is suppressed by powers of H/m φ . This simplified the integration by linearizing the phase of the Bogoliubov integrand, yielding a form reminiscent of a Fourier transform.
We derived a time model of inflaton dynamics and used it simplify the β k integrand to a form amenable to exact analytical integration. This was done using the adiabatic turning points to create an envelope of the oscillatory motion. A differential equation for the fast oscillatory phase was then obtained using approximate energy conservation. By treating the slowly varying components as constant, and using a perturbative expansion of the inflaton potential, the phase was obtained as a function of adiabatic invariants.
A closed form production rate proportional to a sine-cardinal function was derived, with a resonant peak when the inflaton oscillatory frequency ω * equals the time derivative of the Bogoliubov phase E k . We then integrated the production rate to obtain the spectrum and number density. For m χ < m φ and large k, these results matched the exact Bogoliubov integration.
For k > m φ a end , the rate limits to the form of Fermi's golden rule for the tree-level graviton-mediated process φφ → χχ in a frame where both inflatons are at rest. In this case, ω * and E k can be interpreted as the energies of the φ and χ particles, respectively. The delta function in the production rate allowed an analytical evaluation of the time integral, which yielded simple closed form approximations of the spectrum and number density. These made O(1) corrections to the predictions of the Gaussian spectral model of ref. [44].
This formalism demonstrates a correspondence between the Bogoliubov and scattering methods of computing gravitational particle production. The predicted spectrum and number density are the same, up to factors of 2, as the results of papers that considered a gas of inflatons at rest undergoing gravitational annihilation [41,42]. These factors of 2 are not currently understood, but we suspect that the inflaton gas must be treated as an appropriate coherent state in the scattering picture to obtain an exact correspondence.
The noisy behavior of the numerical spectrum as a function of k (see for example numerical results of ref. [8] and figure 4 of the present work) remains unexplained and is a subject of ongoing research by our group. Preliminary results suggest this is due to quantum interference between different resonances of the Bogoliubov integral, and have interesting implications for the particle scattering method of computing gravitational particle production. It would also be of great interest to find cosmological observables that are sensitive to the details in the momentum spectrum of the gravitationally produced particles, which is one of the main results of this paper.

ACKNOWLEDGMENTS
This work was supported in part by the Ray MacDonald Fund at UW-Madison.

Appendix A: Adiabatic invariant
Here we give a clarification of the adiabatic invariant construction presented in ref. [73]. In addition to the derivation of the conserved adiabatic invariant expression of eq. (A25), the error incurred over a long time period T / (where T is the period of fast periodic motion and is a small adiabaticity parameter) is explicitly evaluated in eq. (A34) using eq. (A26). We also compare it to the canonical transformation approach of ref. [74].
Consider an external influence on the system represented by the function λ(t). The Hamiltonian of the system is denoted as H( q, π, λ(t)), where π is the conjugate momentum to q. Suppose this function is slowly (adiabatically) changing such that on a time scale T . The energy E(t) ≡ H( q(t), π(t), λ(t)) is not conserved because of the time dependence of λ(t). We write the time average of dE/dt at time v over a time period T as where the leading order error ∆ 1 is estimated as through a Taylor expansion. Next, restrict to a 1-dimensional motion which becomes periodic in t when λ(t) becomes independent of t as λ(t) = λ(v). Define Π(q, E, λ) to be the solution to as an algebraic identity. Using this, we can define the time period to be where C T (v) represents a chosen 1-dimensional path where q C (t) returns to itself periodically with λ(t) = λ(v) fixed, i.e. q C (t) is a solution to the equation of motion with λ(t) = λ(v). From eq. (A4), we conclude as E and λ are independent variables. One can use this to write where we made use ofq = ∂ π H (q, Π (q, E, λ) , λ) just as in eq. (A6). Since there are 2 values of t for a given value of q in the time interval [v, v + T ], we will call one branch of eq. (A2) integral "branch 1" and the other "branch 2". Putting eq. (A9) into eq. (A2) gives wheret(q ) is the implicit solution of eq. (A9) such thatt(q(t)) = t.
As the external parameter is fixed at λ = λ(v) without time variation for which q C (t) is defined below eq. (A7), a small error is incurred by putting in approximate solutions with λ fixed. Therefore, we define the error ∆ 4 as where which is still a function of time. Putting this in to eq. (A10) gives which implies where we replaced T (v) in the second line using eq. (A5). Taking a derivative of eq. (A4) with respect to E yields is conserved when ∆ n are all neglected. To be more explicit, the identity is when errors ∆ n are neglected. The correction to Q(t) is to leading order in ∆ n . Let us compute the change in Q over a long time period over which λ changes significantly, denoted as to leading order in 1. 14 For a generic non-adiabatic-invariant quantity O NAI (t) that derives time dependence from λ(t), the change over the long time period ∆t(v) is where the time average has been estimated as Ȯ NAI v ∼ O NAI /T (v) which exactly vanishes in the limit that λ(t) becomes time independent. Note regardless of how small the time variation is, the change in O NAI over the long time period ∆t is large. This is in contrast with the case of the adiabatic invariant for which the change is where we have assumed, for example,Ė(t) ∼λ(t)E/λ,λ(t)/λ ∼ O( 2 )/T 2 , Q ∼ ET , and ∂ 2 λ H( q, π, λ(v)) ∼ H/λ 2 , based on smoothness of these functions. Simplifying the large piece ∆ 3 , we find and thus to be the variation in the adiabatic invariant charge Q over a long time period ∆t ∼ T / . This makes manifest the importance of the periodicity of q motion. If the motion were not periodic, then we would conclude leading to a result analogous to eq. (A28). However, for periodic motion relevant to the adiabatic invariant construction, the first term in eq. (A31) turns into resulting in which shows that Q is conserved to leading order even on long time scales over which λ changes. This matches the more subtle analysis of ref. [74] where the key idea is to first derive the leading order approximation , Ω(t), λ(0))] q=q(Ω,J,λ(0)) (A35) where F 1 (q, Ω, λ(0)) is the generating function of the canonical transformation to a special set of variables: action angle variables (Ω, Q) for the time independent problem except with the constant parameter λ(0) lifted to a function of time λ(t). 15 More explicitly, one first solves the phase space symplectic form preserving map (q, p; λ(0)) ↔ (Ω, Q; λ(0)) for the time independent problem, 16 use this to find p(q, Ω; a(0)) and Q (q, Ω; a(0)), make the replacement a(0) → a(t), and solve the differential equations p(q, Ω; a(t)) = ∂ q F 1 (q, Ω, a(t)), Q (q, Ω; a(t)) = −∂ Ω F 1 (q, Ω, a(t)), to construct F 1 (q, Ω, a(t)). Because of the action-angle variable choice, one can show that , Ω(t), λ(0))] q=q(Ω,J,λ(0)) and ∂ ∂Ω [∂ λ F 1 (q(t), Ω(t), λ(0))] q=q(Ω,J,λ(0)) , are periodic with period T such that eq. (A35) vanishes.
Appendix B: A perturbative expansion of the Hubble rate time dependence In this section we will develop a formalism to estimate the time dependence of the adiabatic energy V m = V (φ C,± ) by using an expansion of the inflation potential from leading order behavior at its minimum. As the time dependence is given by the first few terms of this expansion, this formalism is broadly applicable to many models of inflation.

Integral equation
If we want to solve for V m as a function of scale factor a slow , we want to solve obtained from eq. (20) for the potential at the oscillation maximum V m where 15 Here Q is the conjugate momentum to the angle variable Ω in contrast with the usual notation of calling Q as J. This notation change was natural since Q is often used to denote a conserved charge. 16 Symplectic form preservation is a sufficient condition for the construction of canonical transformations for a time-dependent Hamiltonian, even if the form-preserving transformation does not lead to any special simplification in the general timedependent case: i.e. when one uses the time-independent problem derived action-angle canonical transformations on the time-dependent system, the action variable Q is not constant, but the transformation is still canonical. and we have shifted the potential such that minimum of the potential is at V (φ min ) = 0. For potentials V (φ) where the inverse function V −1 (V ) is simple enough for the integration to be executed in a closed form (such as monomial potentials), eq. (B1) is useful. Since Q is a constant, we can take ratios to write .
(B4) An interesting information offered by this expression is that even though V m will generically not scale as a power law with a slow , there is always some function of V m only that will scale as a 3 slow /a 3 end . As an example of using eq. (B1), consider the potential of eq. (30). We can easily work out making eq. (B1) evaluate to which allows us to write by taking the ratio as in eq. (B4) a simple expression matching the result of eq. (33). In the limit n → ∞, the potential is infinitely flat near the minimum, corresponding to an effectively massless scalar energy density behavior. Indeed, it is well known that for kination dominated fields, the equation of state is such that the energy density scales as a −6 consistently with eq. (B7). As a second example, consider eq. (141) for which and which can be inserted into eq. (B1) to evaluate V m implicitly in terms of hypergeometric functions. However, writing the formal expressions are not as illuminating as the original integral expression of eq. (B4) with the substitutions of Eqs. (B8) and (B9). This in turn is not as illuminating as the Taylor expansion derived expression of eq. (131) whose derivation we turn to next.

Small √ Vm expansion
The equation (20) for the adiabatic invariant can be written as where we defined the dimensionless parameters The inversion of this result gives Using the relation H slow = V m /3M 2 P yields a perturbative expansion of the Hubble rate as We can summarize in terms of the potential coefficients as To gain intuition for this formalism, consider a toy potential of eq. (38) with m φ = m = M giving and φ min = 0. The potential Taylor expansion coefficients are which gives consistently with an easily obtainable exact solution. 17 Note that in situations where Aφ 3 vanish, α 3 and 2V m /m 3 would be absent. This also illustrates that α n is not necessarily small, partly due to the hierarchy between M P and the non-gravitational dynamics scales m and A. As long as the potential is analytic, the Taylor expansion of eq. (B12) can be exact, which means that α n has no requirement of being small. Eq. (B43) becomes where one sees that the corrections to the u −3/2 coming from the cubic interactions do not depend on the M P /m hierarchy of the intermediate steps coming from α n . Although the growing pure numbers such as 3045/128 seem alarming, the extra power of A 2 Q/(2πm 5 u 3 ) (which can be orders of magnitude smaller than 10 −1 in practice) leads to a suppression, particularly at large times when u 3 grows large. For a cosmologically relevant example, see eq. (147).