Quantum interference in gravitational particle production

Previous numerical investigations of gravitational particle production during the coherent oscillation period of inflation displayed unexplained fluctuations in the spectral density of the produced particles. We argue that these features are due to the quantum interference of the coherent scattering reactions that produce the particles. We provide accurate analytic formulae to compute the particle production amplitude for a conformally-coupled scalar field, including the interference effect in the kinematic region where the production can be interpreted as inflaton scattering into scalar final states via graviton exchange.


Introduction
During the period of coherent oscillations of the inflaton field following the quasi-de Sitter (quasi-dS) phase of inflation [1][2][3][4][5][6][7], particles (including dark matter candidates) may be produced via gravitationally-mediated nonthermal scattering processes in addition to the comparatively well-studied inflaton-decay and thermal-scattering processes . Numerical investigations of gravitational particle production (GPP) employing the Bogoliubov approach have displayed unexplained oscillations as a function of the wavenumber k in the final phase-space distribution f χ (k) of the produced particles [61][62][63]. For example, figure 1 of ref. [62] showing the final phase-space density of dark-matter particle production in a hilltop inflationary model displays large oscillations that resemble numerical noise. Similar large oscillations in the final phase-space density can be seen in the right-hand panel of figure 1 in ref. [63] for GPP of the helicity-1/2 component of a spin-3/2 Rarita-Schwinger field.
In this paper we explain these oscillatory features as the result of a quantum effect arising from an interference of different amplitudes, which are analogous to gravitationallymediated nonthermal scattering processes, 1 nφ → 2χ for n ≥ 1. Typically the 2φ process dominates nonthermal scattering production, but it has recently been pointed out [64] that the nφ processes with n = 2 may also be important. Most of the effect comes from interference of 2φ with the next leading amplitude, which is 3φ if cubic interactions exist and 4φ otherwise. We compute analytically the scattering contribution to the Bogoliubov amplitude including the interference terms, and find the results compare well with numerical computations. We also give a less technical semi-quantitative estimate of this interference amplitude based on a coherent scattering picture of a modified Boltzmann evolution. In this latter picture, the interference arises because the initial macroscopic inflaton scattering state can be viewed as a cold coherent superposition of nφ states, e.g., c 1 |φφ +c 2 |φφφ , such that the interference arises from |c 1 χχ|U |φφ + c 2 χχ|U |φφφ | 2 where U schematically depicts a time evolution operator which is made more precise in this paper. Note that we write nφ → 2χ to denote the net energy flow from the φ field to χ field, but this can be different from underlying S-matrix amplitudes. For example, φ → 2χ has a contribution from the φφ → φχχ scattering process.
Although the quantum nature of the inflaton coherent-oscillation induced GPP has been known (e.g., see [65]), the present paper extends the previous ideas to graviton-mediated scattering, and to our knowledge is the first to articulate clearly and to compute analytically the quantum interference effects. It also clearly explains the previously unexplained "noise" in the particle production spectrum seen in the literature (see, e.g., [62,63]). The application of a novel perturbative technique to solve the background inflaton dynamics is a technical highlight of this paper.
The order of presentation is as follows. In section 2, we give a brief review of the GPP computation using the Bogoliubov transform technique. In section 3, we derive an analytic formula for the relevant Bogoliubov coefficient using a novel perturbation theory technique and a stationary phase approximation. The result is a sum of amplitudes analogous to nφ → 2χ, with section 4 presenting explicit results for n ≤ 4, and section 5 discussing the quantum interference between amplitudes. In section 6, we compare the analytic results with numerical computations. In section 7, we interpret the interference as a novel contribution to the Boltzmann collision equation arising from the initial inflaton field being a macroscopic state described as a coherent superposition of nφ states. We then conclude in section 8 with a summary and outlook.
The appendices contain some of the supporting technical details of this work. In appendix A we describe the background field evolution in polar coordinates. In appendix B, we summarize the novel perturbation technique used to solve the inflaton dynamics with asymptotic series involving functions λ(t) and θ(t) that describe slow and fast time scales, respectively. Appendix C explains the technical details of evaluating the terms formally set up in the stationary phase computation in section 3 using this technique. In appendix D, we remind the reader how the statistical ensemble factor enters the usual collision integral of a Boltzmann equation in a manner that is in contrast with the picture of section 7.

Gravitational particle production
Here, we focus on a background spacetime described by standard Einstein gravity with a spatially-flat Friedmann-Lemaître-Robertson-Walker (FLRW) metric ds 2 = dt 2 −a 2 (t)|d x| 2 = a 2 (η) dη 2 − |d x| 2 where dη = a −1 dt is conformal time. The dominant energy-momentum tensor for the dynamics of the scale factor a(t) comes from a minimally-coupled real scalar inflaton field φ with mass m φ and a slow-roll inflationary potential V (φ). We will assume that m 2 φ ≡ d 2 V (φ)/dφ 2 | φ=v = 0, where v is the minimum of V (φ) during the inflaton's coherent oscillation phase after the quasi-dS phase, and we will also assume that the nonlinearities in V (φ) can be captured as a Taylor expansion about φ = v.
The inflaton potential will be parameterized as where M P = 1/ √ 8πG is the reduced Planck mass. When specific examples are needed, we will consider two inflaton models denoted by where we take v = 0 as the minimum for the Quadratic model and v = M P /2 for the Hilltop model. 2 Note that α 3 = α 4 = 0 for the Quadratic potential while α 3 = 5 √ 6 and α 4 = 155 for the Hilltop potential.
We augment the standard inflationary picture with a scalar "spectator" fieldχ whose action is given by where χ ≡ aχ is the rescaled field, m χ is the particle mass, and R = −6∂ 2 η a/a 3 is the Ricci scalar. Following the usual procedure (e.g., [65][66][67][68]), we promote the scalar field to an operatorχ that satisfies the canonical equal-time commutation relations. The field operator is decomposed into mode functions χ k labeled by wavevector k aŝ where the mode functions satisfy the normalization condition χ k ∂ η χ * k − χ * k ∂ η χ k = i, and the creation and annihilation operators satisfy the canonical commutation relations. Due to the action in eq. (2.5), the mode equation is is the angular frequency of the k th Fourier mode. 3 The vacuum state |0 is defined aŝ α k |0 = 0 for all k, and particle creation is generated by the time-dependence of ω k (η). While one can solve the mode equation directly given initial conditions, for our purposes we use the Bogoliubov parameterization. The mode functions are expressed as ) where α k and β k are the Bogoliubov coefficients, which decompose the mode function into positive and negative-frequency components, respectively. The nearly-adiabatic conditions in the far past motivates the Bunch-Davies initial condition such that α k = 1 and β k = 0 at initial time t = t i . In the evolution of χ k from the initial negative-frequency solution, a positive-frequency component may appear, signaling particle creation. In the far late-time, the number density of produced particles is given by where f χ (k, t) ≡ |β k (t)| 2 denotes the produced χ-particle phase-space density. We therefore seek a solution for β k to compute GPP. The time-evolution of the Bogoliubov coefficients α k and β k can be written as 4α as is done for example in [11,65,69]. For the case of a scalar χ field, we use the definition with Ω k defined in eq. (2.9). The background evolution (assumed driven by the dynamics of the inflaton) enters the determination of N k through a, H, R, andṘ, while the spectator field enters through k and m χ . Setting α k ≈ 1 in eq. (2.12), we write which is valid for |β k | 1. This important integral expression is the staring point for the main results of this paper.
In situations where there are coherent oscillations of the inflaton, some contributions to GPP can be interpreted as coming from scattering of the inflaton quanta into χs via graviton exchange [16,21,59,62]. Below, we explain a novel computation of this coefficient that gives not only the amplitudes of nφ → 2χ during the inflaton coherent oscillations, but also the interference between amplitudes with different n ≥ 1.
In the scenario with coherent oscillations of the inflaton field after the quasi-dS phase, the Hubble expansion rate has two broad classes of components: intuitively, H = H slow + H fast , where the oscillatory fast component is smaller in amplitude but varies on a larger frequency scale compared to the monotonically decreasing slow component. The respective scales are d dt ln H slow ∼ H and d dt ln H fast ∼ m φ , and therefore H < m φ is required for a meaningful distinction. This condition holds in the oscillatory era for most single field inflationary models. This decomposition is made precise using the formalism summarized in appendix B. This approach differs from that of ref. [64], where time was partitioned into bins of size H −1/2 slow m −1/2 φ to find theφ behavior applicable to each of those bins. That approach is suitable for the particle production computation without the interference, but to find the interference, we need to keep track of the time phase of the Bogoliubov coefficients across different time bins ∆t. The approach we will take below is to use a novel method of expanding the background field evolution in H slow /m φ .
This formalism introduces functions λ(t) and θ(t) that partition the slow-time and fast-time dependence of a general quasi-periodic function of time such as N k (t). The monotonically decreasing slow-time variable λ can be thought of as where all quantities with a subscript-e index will refer to its value at the time t e when the quasi-dS era ends, also referred to as the end of inflation. For the two inflation models we will consider, Quadratic and Hilltop, H e /m φ 0.5 and 0.03, respectively (see table 1), and therefore λ(t) is less than unity for t ≥ t e . The fast-time variable θ can be thought of as a diffeomorphism of time to a monotonically increasing phase function such that In short, λ and θ describe time scales of H −1 and m −1 φ , respectively. This partitioning is the basis for the novel perturbation technique appearing in ref. [70], which allows us to resum secular effects and track β k (t) accurately for a long time (a time much longer than m −1 χ ). For example, the Hubble expansion rate is expanded systematically as where h (θ) contains the fast-time behavior as a sum of sinusoids that depend on integer multiples of θ. The higher-integer frequency components become increasingly negligible as they generally come with higher powers of λ(t). This accurate tracking for a long time is useful for capturing our sought-after interference effects, which develop on a time scale ∆t H −1 m −1 χ for the k/a e > m φ modes. The time evolution of these functions is defined byλ where the initial conditions and details of the O(λ 2 ) corrections are irrelevant for the arguments of this section. The λ, θ decomposition is explained in further detail in appendix B.
To describe the interference effects, we write β k as a sum over β (n→2) k contributions by partitioning the time dependence of the integrand from eq. (2.14) into λ(t) and θ(t), where an exact definition of β (n→2) k is presented below. In summary, the n → 2 label alludes to the nφ → 2χ scattering process, with the integer n denoting the frequency that stems from an oscillatory dependence on nθ. As noted below eq. (3.3) and shown in section 4, larger n are increasingly suppressed, and therefore β k is well approximated by the first few terms of this sum. This property motivates the use of this formalism, and is analogous to the suppression of higher particle number processes in perturbative QFT.
We now present definitions used in the specification and computation of β . Given a general quantity X(λ, θ), we define the convention where X (n) (λ) is the n th -frequency component of X. The slow-time component X s is defined as the n = 0 term, which does not contain any fast-time information, i.e., no θdependence, and the fast component X f is defined as the remainder (i.e., the sum of the n = 0 terms). As an example, it will be useful to separate the phase-factor Ω k into a slow component Ω s,k (λ) ≡ Ω (0) k (λ) and a fast component Ω f,k (λ, θ) ≡ Ω k (λ, θ) − Ω s,k (λ). We define β (n→2) k for the resonant scattering situations of current interest such that 8) and N k ≡ N k e −2iΩ f,k , with the transformation of t to λ, θ dependence as well as the expression of N k (λ, θ) described in appendices B and C. 5 The exponential in eq. (3.8) hints at the relation to the amplitude for nφ → 2χ as its phase is stationary when nm φ ≈ 2 k 2 /a 2 + m 2 χ , which corresponds to the energy condition of n inflatons at rest annihilating to produce two χ particles with momentum k/a. This correspondence to scattering motivates referring to β (n→2) k as the n → 2 resonance component, and suggests evaluation using the stationary phase approximation.
We compute the n → 2 resonance component of the Bogoliubov coefficient using the stationary phase approximation, which will ultimately lead to an expansion in powers of k −3/2 . For the purposes of explaining the computation, we write as the total (complex-valued) phase for the n → 2 resonance. The phase is stationary wheṅ Ψ (n) k as the resonance time, 6 which will usually have a small imaginary component due to the complex nature of the phase. The phase is expanded as which is the starting point of the stationary phase approximation. We now define a new variable z such that the quadratic term in the exponential becomes − 1 2 z 2 , and evaluate using an expansion of Gaussian integrals as where the coefficients of the higher powers of z are treated as increasingly negligible, an assumption that will be justified shortly. In going from eq. (3.10) to eq. (3.11), we moved the contour into the complex plane in addition to changing the integration limits. In particular, the contour was rotated by approximately 45 degrees asΨ (n) k at resonance is dominated by its imaginary component due to the first two terms of eq. (3.9).
We need an expansion parameter to truncate the expansion of eq. (3.12), which we will see is proportional to k −3/2 by the following argument. Sinceλ ∼ λ 2 and the derivatives of the phase depend on λ, we know that ∂ t Ψ (n) k ∼ λ −1 and therefore the z coefficient in eq. (3.11) scales as λ /2−1 , which is suppressed for ≥ 3. Hence, the next step is to evaluate φ was determined by definition in eq. (A.12). Therefore, we estimate φ a e for n = O(1). Since our present computation is focusing on scattering of particle-like modes at the end of the quasi-dS era, λ(t is small for all the modes of our present interest and thus is naïvely a useful expansion parameter. This statement will be made more sharp in section 6. We would like to solve for λ(t (n) k ) itself as a function of k. Hence, we define a separate expansion parameter for the stationary phase approximation as where ε is a bookkeeping parameter inspired by the smallness of eq. (3.13). We can parameterize a perturbation series solution of λ(t where the constant coefficients r (n) k,j are determined by solving the stationary phase conditioṅ Ψ (n) When using the replacement of eq. (3.15), it is important to write k in terms of ε −2/3 κ n using eq. (3.14) to cancel out the fractional powers of ε that appear due to a ∼ λ −2/3 . This is equivalent to assuming k and am φ have the same magnitude at resonance. Afterwards, we can use our solution to evaluate the phase-derivative coefficients appearing in eq. (3.12). Some of the technical details of this computation are given in appendix C.

Analytic formulas for the Bogoliubov coefficient
In this section, we explicitly list the analytic amplitudes β (n→2) k for n ∈ {1, 2, 3, 4} solved by the procedure described above, with the k-dependence expressed as an expansion in κ −3/2 n as defined by eq. (3.14). We choose a conformally-coupled (ξ = 1/6) scalar χ field because of the relative simplicity of the source of nonadiabaticity. To make the interference phase more manifest, we express our results as (4.1) Up to a global phase that is independent of n, which therefore affects neither the interference nor the magnitude of β k , the leading terms for the phase can be written as where r χ ≡ m χ /m φ , H and Ξ are boundary conditions defined in eqs. (A.10) and (A.11), respectively, and 2 F 1 is the hypergeometric function. We will give a physical interpretation of this leading order phase in section 7. If we define ∆Φ k,leading and r χ = m χ /m φ , then we can write our results as i , z (n) are merely notational variables to allow a visually manageable display of the results, with their explicit values given in appendix E. Given the generalized nature of α n as defined in eq. (2.2), these results are applicable to any inflationary potential which can expanded as a polynomial with a positive quadratic term at its minimum. Since A (2→2) k will be the dominant term, we have shown it to higher order in the expansion. Note that the higher order terms in these formulas are organized by time reversal symmetry, and therefore have relative corrections that go as ε 2 κ −3 n instead of εκ respectively. This method of writing the boundary conditions allows a cleaner set of analytic expressions. Numerical solutions were used to obtain the values found in table 1 for two inflationary models of interest, and this was done mostly for accuracy when comparing the analytic β k with numerically computed β k . Analytic expressions can be obtained for these integration constants as an expansion in H e /m φ for standard slow-roll inflationary scenarios entering the coherent oscillations period. The calculations were done without choosing a particular scheme, a type of gauge choice concept that is particular to our computational formalism briefly described in appendix B. Instead, the scheme dependence was kept general throughout and completely cancelled out in the final result. Checking scheme independence of observables was a robust tool to verify different steps of the calculation. Another feature to note is that all amplitudes vanish as they should when m χ /m φ → 0 since we are considering the conformally-coupled case. A related feature is that the leading phase Φ (n→2) k,leading diverges in the limit that m χ /m φ → 0. Let's now compare the squared amplitudes between the current computation and an earlier work by some of the present authors [64]. From the latter, we have the estimate where the definition ofH end and the above equation for f χ are given by eqs. (8.13) and (8.17) of ref. [64], respectively. This can be compared to our eq. (4.3c). For the Hilltop model of eq. (2.4), the leading expressions differ by a factor of where the value of a is found in table 1, and the value ofH end for this Hilltop model is given by eq. (8.27) of ref. [64]. The difference between a e and a end (also H e and H end ) is a result of different definitions for the end of inflation. This ratio can be used as an estimate of corrections that this paper represents to the computations of ref. [64] as far as the non-interference piece is concerned.

Discussion of the interference
Now, let's consider the interferences arising from the results of section 4. To focus the discussion to the physically most significant case, consider the interference between 2 → 2 and 3 → 2 amplitudes: is real, the interference phase between these two processes comes from where r χ = m χ /m φ , and ∆ is defined as m as numerical coefficients that depend only on the inflaton potential interaction strengths α 3 and α 4 , as can be seen in appendix E. The term proportional to 2κ k ), and each of these terms with the respective coefficients are effectively a rewriting of the resonance times. The hypergeometric functions correspond to the 2Ω k phases appearing in eq. (2.12) evaluated at the respective resonance times. Equation (5.2) is one of the main analytic results of this present work.
The ∆ term contains the leading higher-λ power correction to the leading stationaryphase result. This contains the nontrivial corrections to the phases coming from the cubic and quartic interaction terms of the inflaton potential: i.e., it depends on α 3,4 . It vanishes in the large k/m φ limit because this is just the property of an asymptotic expansion through the stationary phase method.
In section 7, we will discuss how the phases can be interpreted in terms of phases accumulating through the Hamiltonian energy driven time evolution. In this intuitive picture, for a given time interval, the inflaton background field self-interaction and self-gravitational interaction change the accumulated phase of the inflaton interpreted as a collection of oneparticle states because of the change in the effective free propagator Hamiltonian energy. For example, in the parameter region of {m φ m χ , α 3 = 0, m φ H}, one can easily check that ∆ increases as expected from the intuition that the steepening of the potential by the quartic potential contribution increases the effective oscillation mass. The κ −3 2 correction term in |A (2→2) k | 2 of eq. (4.3c) increases with m χ increasing, although the physical interpretation of this increase is not as obvious.
In the more generic region in the parameter space, ∆ is not monotonic with increasing m χ /m φ . For example, ∆ goes through a zero as m χ /m φ is increased if α 3 /α 4 O(1). Since ∆ generically diverges as m χ → m φ and decreases with increasing m χ for small m χ /m φ , there can be two zeroes if ∆ > 0 when m χ /m φ = 0 and α 3 /α 4 O(1).

Numerical examples
In this section, we employ the analytic results of section 4 to study GPP and quantum interference for two specific models of inflation: the Quadratic Potential model from eq. (2.3) and the Hilltop Potential model from eq. (2.4). We evaluate the absolute value of the Bogoliubov coefficients |β k | using the analytic expressions for β k in eqs.  4). We consider a range of dimensionless comoving wavenumbers k ∈ (10 −2 , 10 2 ), where we've set a e H e = 1 such that the modes with k = 1 leave the horizon at the end of inflation.
Our results for the Quadratic Potential model are presented in figure 1. The blue-dotted curve corresponds to our leading-order analytic calculation |β (2→2) k | on both the upper and lower panels, while the red-dashed curve on the lower panel includes the first sub-leading correction |β Since the quadratic potential has a Z 2 symmetry, φ → −φ, the nφ processes with odd n have vanishing amplitudes: e.g. β The spectrum at large k is approximately a power law |β k | ≈ |β (2→2) k | ∝ k −9/4 , but closer inspection reveals a sub-leading component that oscillates as k is varied. These oscillations are explained in this work as an interference effect. Using eq. (4.2), the oscillation period ∆k is controlled by the variation in the phase with respect to k, and can be explicitly written as for any n 1 and n 2 . For the Quadratic model used in figure 1, this evaluates to ∆k 2a e H e for n 1 = 2 and n 2 = 4. For comparison, the gray curve shows the result of calculating |β k | by direct numerical integration of the mode equations. The analytic results derived here agree very well with the numerical integration at large k in both the average power-law behavior and the oscillatory features. This agreement can be viewed as a validation of our analytic approximations. The exponentially dropping black curve 7 by contrast highlights the power-law behavior |β k | ∝ k −9/4 coming from the oscillating inflaton field that drives corresponding oscillations in the scale factor.
For the Hilltop Potential model, our results appear in figure 2. Once again, the leading power-law behavior at large k is |β k | ≈ |β (2→2) k | ∝ k −9/4 as seen from both the direct numerical integration (gray-solid) and our analytic approximation (blue-dotted). The subleading oscillatory components (green-dashed and red-dot-dashed) have a richer structure in this model, which is evident by comparing the lower panels of figures 1 and 2. This behavior can be understood as follows: for the Hilltop Potential model the components β (n→2) k have similar amplitudes with increasing n, leading to a pronounced interference pattern, whereas the amplitudes decrease more rapidly in the Quadratic Potential model, and the interference is dominated by just the first two terms. Moreover, since the Hilltop Potential model does not have a Z 2 symmetry at the minimum of the inflaton's potential, the processes with an odd number of inflatonsφ → 2χ, 3φ → 2χ, and so on -are not forbidden. It turns out that β (1→2) k amplitude is numerically less important than that of β (3→2) k for the interference partly owing to the suppression of (κ 3 /κ 1 ) 15/4 1 (see eq. (3.14)). By including up to the sub-sub-leading order in our analytic calculations, |β we obtain the red-dot-dashed curve that matches the result of direct numerical integration (gray-solid) very well at large k.  Figure 1. The Bogoliubov coefficient |β k | as a function of comoving wavenumber k (with a e H e = 1) in the Quadratic Potential model. We assume a conformally-coupled (ξ = 1/6) scalar spectator with mass m χ that experiences GPP due to an expanding spacetime background driven by an inflaton field φ on a quadratic potential with mass m φ . Top: We calculate |β k | using the analytic results derived in this work (blue-dotted) and using direct numerical integration of the mode equations (gray). Note that |β k | scales as k −9/4 at large k. As a comparative contrast to this power law behavior in k, the black curve shows an approximate expression for |β k | for GPP in a matter dominated (MD) universe, extrapolated to lower k values (beyond the range of strict validity) for visual completeness of the exponential behavior. Bottom: The Bogoliubov coefficient exhibits an oscillatory feature in the large k power-law tail of the spectrum, which is explained in this work as a result of quantum interference between β  Here we note the limits of applicability of our analytic results, using β −800, and the magnitude of the correction is approximately 800κ −3 2 . An upper limit on the magnitude of the correction results in a lower limit on κ 2 , which, in turn results, in an m φ /H e -dependent lower limit on k (see eq. (3.14)). For the figures we have assumed that the next-order corrections to β (2→2) k are no more than 30% (since the lower limit on k only depends of the third-root of the correction limit, the result is relatively insensitive to the choice of 30%). From the figures it is clear that the k −9/4 behavior extends to k somewhat lower than the cutoff in the convergence of our expansion.

A heuristic derivation
Our aim in this section is to describe semi-quantitatively the Bogoliubov computation of the resonance-induced GPP in terms of an approximate S-matrix perspective by showing how the Boltzmann equation would need to be modified to capture the interference effects.
Here, we will focus on the interference of 2 → 2 and 3 → 2 scattering as this is often the most interesting case, with other generalizations being straightforward.
Consider an incoherent gas of N ∼ ρ e V 3 /m φ number of φ particles, where ρ e ∼ M 2 P H 2 e is the energy density and V 3 ∼ H −3 e is the 3-volume of a large box approximating the causal Hubble patch. Usually, one first decoheres this large N state system into an ensemble of 2φ → 2χ and 3φ → 2χ, and then considers each process statistically independent. In this case, the macroscopic particle production of χ is described by a semiclassical 1-particle χ distribution obtained from integrating the collision term aŝ where U (t 2 , t 1 ) is the time-evolution operator from time t 1 to t 2 , and S n factors are φ initialstate dependent weighting factors (generalization of Bose-Einstein distribution), eventually leading to the cross section picture of the usual Boltzmann equations as shown explicitly in appendix D. This treats "typical" 2-body scatterings and 3-body scatterings to be additive incoherently. However, this type of computation neglects the nontrivial interference that can occur from Schrödinger time evolution phases between different scatterings.
Hence, we arrive at the main idea. The scattering perspective that we will construct below will simply replace the nonadiabatic period during which the 2χ particle frequencies are in resonance with an approximate S-matrix scattering description. The different scattering events are diagrams (e.g., see figure 3) that interfere because of the coherence of the waves entering the interaction region approximated by an S-matrix. This will allow us to compute the interference phase using the wave free-propagation phase. Thus, before we Figure 3. The Schrödinger propagator phase difference between 3φ → 2χ and 2φ → 2χ scatterings leads to interference. The disk region of diameter δt represents the usual collision region of Boltzmann equation, which is typically treated with an S-matrix taking the formal limit δt → ∞.
The observable interference phase of χ is the Schrödinger-picture free-particle propagator between t 3 and t 2 . The φ-interference phase can contain t e information as 3φ propagation phase from t e to t 3 does not cancel the 2φ propagation conjugate phase from t e to t 3 .
describe the scattering, let's divide the time period [t e , t f ] into 3 regions: where t n is the time at which nφ → 2χ resonance occurs, i.e., 2 k 2 a 2 (tn) + m 2 χ ≈ nm φ , which is the analog of the timet (n) k that satisfies the stationary-phase condition from section 3. From a scattering perspective, we work in the Schrödinger picture with metric inhomogeneities in time, with φ treated as a quantum field, and the interaction Hamiltonian coming from the metric fluctuation coupling to the φ energy-momentum tensor. To describe the spatially homogeneous classical inflaton field, imagine setting up a normalized coherent state |φ) ∼ (a state containing a macroscopic number of particles) 8 at time t e such that (φ|U (t e , t) φ(x, t) U (t, t e )|φ) = φ EOM (t) (7.2) for t > t e , where φ EOM (t) is the solution to eq. (A.1). Note that the quantum phase of U (t, t e ) has turned into the classical phases embedded in φ EOM (t), approximated as nm φ t for integers n. This is one source of the interference phase as we will see below.
We will assume for the semi-quantitative discussion that this coherent state can be generalized straightforwardly to the effective FLRW background of an expanding box: ds 2 = dt 2 − a 2 slow (t)|d x| 2 . Instead of treating H slow as a Minkowski graviton effect, the background is treated as an expanding box even in the scattering picture because a purely Minkowski treatment is inefficient in explaining k/a slow dilution.
The normalized coherent state |φ) of eq. (7.2) can be decomposed as a superposition of normalized wave packet states over r numbers of φ particles, written as [72] |r, controls the φ particle wave packets with central momenta {q i } ≡ {q 1 , ..., q r }. For illustration, suppose the initial state in the notation of eq. (7.3) is decohered into clusters of 2 and 3-particle states described by a density matrix ρ = ψ S(ψ)|ψ ψ|, where |ψ ≈ |0, t e χ ⊗ |Ψ φ ⊗ |0, t e δgµν , (7.4) and S partitions the macroscopic N -inflaton state into an ensemble of coherent superpositions of 2-particle and 3-particle states. The wave packet function F r appearing in eq. (7.3) is assumed to be peaked at close to zero spatial momentum since the inflatons are assumed to be cold. The amplitudes ζ 2 and ζ 3 control the mixing of 2-and 3-particle inflaton states. The φ state of eq. (7.5) can be intuitively considered a "classical" coherence because it represents a macroscopic state, 9 and the Bogoliubov vacuum does not contain the quantum data for φ in ζ n . However, this "classical" coherence itself is really part of the quantum coherence associated with the Schrödinger time evolution operator just as in photon time-phase coherence in lasers.
With the illustrative partition of eq. (7.5), the analog of eq. (7.1) becomeŝ is a multiplicity factor associated with the partition achieved through the density matrix probability factor S. Because of the resonant behavior, we know The sum over S is a macroscopic number, similar to V 3 where t ± n ≡ t n ±δt/2, with δt as the interaction time, i.e., the time scale of a Boltzmann collision term, which by construction is supposed to be much smaller than the free-streaming time scale. However, δt is viewed in the S-matrix picture as an asymptotically long time scale, as one formally takes δt → ∞ to take advantage of the properties associated with meromorphic matrix elements. 10 This is the usual requirement of the validity of the Boltzmann treatment. Hence, the squared amplitude in the modified Boltzmann collision analog of eq. (7.6) becomes where one notes in eq. (7.9) that the cross-term induced coefficient as part of the interference phase. More generically, the interference phase between n 1 φ → 2χ and in the limit that δt |t n 2 − t n 1 |. The phase of the Schrödinger-propagator independent quantity ζ * n 1 ζ n 2 is apparently independent of n in the case of our particular |β k | 2 computation.
Comparing with eq. (5.2), we see that where the hypergeometric function arises from integrals of the form which used the relationship a ≈ aλ −2/3 and the definition ofλ in eq. (B.3). The hypergeometric function term by itself has a divergent piece as m χ → 0, which is obviously spurious since the left hand side of eq. (7.11) is convergent for finite t n . Similarly, the remaining terms of eq. (5.2) can be identified with the inflaton phase: which also matches the interpretation of Ξ being the phase offset that depends on the properties of the inflaton at the end of the quasi-dS era at time t e .

Conclusions
In this article we report on our study of quantum interference in the phenomenon of gravitational particle production. Our main results appear in section 4. We have derived analytic expressions for the Bogoliubov coefficients β k describing the gravitational 10 Note that A (n→2) k is typically proportional to the scattering amplitude and An at tree level order.
production of conformally-coupled, massive scalar particles during the inflaton's coherent oscillations after inflation. By employing a novel perturbation technique (relying on a nonlinear field redefinition) and a stationary phase calculation, we have expressed β k as a sum over resonant contributions β (n→2) k . Oscillatory features in the spectrum |β k | 2 are understood to result from an interference among the resonant contributions, e.g. |β ; see also eq. (5.2) for details. These analytic results are in excellent agreement with a direct numerical integration of the mode equations; as shown in section 6, the agreement is within a few percent in certain kinematic regions. Our work explains much of the previously unexplained "noise" in numerically-computed spectra, seen for example in refs. [61][62][63]. As we discuss in section 7, the resonant contributions β (n→2) k are related to gravity-mediated inflaton scattering amplitudes nφ → 2χ corresponding to n inflaton particles with mass m φ at rest annihilating to 2 scalar particles with mass m χ < nm φ /2. This work also elucidates the quantum nature of gravitational particle production induced by classical inflaton coherent dynamics.
As noted in section 7, the interference phase can be understood as arising from the free propagator phases of the external legs of the scattering process. This means that the phases are dependent on the kinematics of the inflaton and the χ particles, as well as the scattering times of say n 1 φ → 2χ and n 2 φ → 2χ processes. Unlike the usual scattering situations where n 1 φ → 2χ and n 2 φ → 2χ are incoherent, the coherent oscillation nature of the initial inflaton state allows for the scattering amplitudes to interfere. This interference is efficiently captured using the Bogoliubov transformation formalism.
The modulations of the χ-particle momentum spectrum shown in figures 1 and 2 in principle can be probed by kinematic-dependent subsequent scattering dynamics of χ particles. For example, if interesting motivated scenarios exist for χ particle scattering resonances with judicious energy spacing, the interference pattern of χ energies may lead to enhanced production of final states compared to situations without this interference pattern in the χ particle spectrum. Investigations into possible applications will be left to future work.

A Background field evolution in novel polar coordinates
The evolution of the inflaton field φ is usually described by the second-order equation which is often referred to as the inflaton equation of motion. For our purposes we wish to exchange the second-order differentiation to a set of two first-order equations at the cost of introducing another dependent variable. There is freedom in choosing the two variables; our choice is the Hubble rate H and phase Ξ (both reals), defined in terms of φ andφ as and v is the field value where the potential is minimized. As explained in ref. [70], the change of variables from {φ,φ} → {H, Ξ} is analogous to switching from Cartesian to polar coordinates in phase space, with H representing the radial coordinate and Ξ the angular coordinate.
which is simply the real part of eq. (A.2). As done in ref. [70], we simplify the presentation of our problem by using the change of variables which is equivalent to setting √ 6M P = m φ = 1 and v = t e = 0. Using the expansion of the potential in eq. which allows us to express the EOMs entirely in terms of H and Ξ: ν n H n cos n Ξ , (A.8) for some constants ν n , with ν 1 = 2α 3 , ν 2 = − 7 2 α 2 3 +3α 4 , and so on. We write the derivatives for H −1 and Ξ to highlight that both grow linearly with time if H is small. This will be useful when deriving the constants H and Ξ, which are associated with the boundary conditions at t = +∞ for H and Ξ, respectively.

A.1 Defining boundary conditions in asymptotic far future
This subsection introduces the constants H, Ξ, and a that our formulas for β (n→2) k ultimately depend on. These constants quantify the boundary conditions in the far future as t → +∞, and contain the same information as the initial conditions. This limit is necessary to integrate the equations of motion starting at λ = 0 for our asymptotic expansions. We

A.2 A special property at the end of inflation
Here we note a special property about Ξ e . The end of inflation is defined bÿ a e = 0 ↔ ρ e + 3p e = 0 ↔φ 2 e = V (φ e ) , (A. 15) where the arrows represent equivalence between all three statements. Using the third statement and eq. (A.2), we can show that where s e ≡ sign(φ e − v) = ±1. This has the advantage of specifying the end of inflation in a closed-form and geometric manner using a single variable with no derivatives involved. Examples of both physically distinct solutions for Ξ e can be found in table 1.
Using the geometric expression of Ξ e in eq. (A.16), we can express the approximations for the boundary condition constants from the previous subsection as which only uses the initial conditions of a e , H e , and s e , along with information about the potential up to the quartic interaction, i.e., m φ , α 3 , and α 4 . We expect s e to only appear with quantities such as α 3 that break the reflection symmetry of the potential about its minimum, i.e., V (v + ∆φ) = V (v − ∆φ), i.e., α n = 0 for all odd n. If the potential is symmetric, then H, Ξ − Ξ e , and a must be independent of the s e initial condition. In this case, the action s e → −s e causes both Ξ e and Ξ shift by the same factor of π such that the difference is unaffected.

B Summary of the perturbative asymptotic series formalism
This section summarizes the formalism of ref. [70], in which the equations of motion are solved using asymptotic expansions of the variables as a function of "slow-" and "fast-" time variables λ and θ. As shown in section 3, this allows computation of the Bogoliubov coefficients as an analytic expansion in powers of k −3/2 . We begin by using eq. (A.5) to scale and shift away various constants. It is convenient to solve the dynamics using the polar coordinates H and Ξ defined by eq. (A.2), as the former is used in GPP calculations. We use perturbative series in powers of λ to write them as where h (θ) and ξ (θ) are oscillatory functions that must remain bounded in magnitude to maintain the stability of the expansion. This requirement will determine most of the constants of integration associated with solving for these functions at each order in λ. As one can see from eq. (B.1), the expansion in powers of λ is justified if H/m φ is small. To solve the EOMs from eq. (A.8), we must specify the time evolution of λ and θ. We define the derivatives of both variables aṡ where c j and ω j are constant coefficients. While the EOMs determine c 0 = ω 0 = 1 and c 1 = ω 1 = 0, the coefficients for all j ≥ 2 remain unfixed parameters. It will be shown that the parameter also remains undetermined due to the time translation invariance of the EOMs, i.e., the freedom to choose the origin of the time coordinate. We call every such choice of h 1 and {c j , ω j |j ≥ 2} a renormalization scheme (RS) because of the analogy with the coupling flow equations. Of course, given that this degree of freedom is a diffeomorphism choice, we could have also called it a gauge choice.
We are free to choose h 1 by the following argument. The derivatives in eqs. (B.3) and (B.4) are invariant under shifts in t or θ. Thus, if λ(t) and θ(t) are solutions for a given set of coefficients {c j , ω j }, then λ(t − t s ) and θ(t − t s ) − θ s must also be solutions for any shifts t s and θ s . Using the typical Taylor series, this can be expressed as making the replacements with derivatives at t determined by eqs. (B.3) and (B.4). Applying these shifts to our asymptotic expansions yield equivalent solutions due to the time translation invariance of the EOMs. However, the form of Ξ = θ + O(λ) is violated unless θ s + t s = 0. Therefore, we have one remaining symmetry parameter, which we will denote as δh 1 = 3 2 t s = − 3 2 θ s because of the shift it induces in h 1 . In summary, if we apply then our asymptotic solutions transform as which is equivalent to a change of h j and ξ j . Fixing the value of h 1 breaks this shift symmetry and therefore acts as an additional RS parameter. Our results in section 3 will be computed with a consistent truncation to render the results explicitly independent of the renormalization scheme. (Detailed proof of this appears in ref. [70].) This gives us the freedom to choose the RS such that for all ≥ 1, which is convenient as the expressions tend to be relatively compact in this scheme. The relevant results to O(λ 3 ) are given by h 1 = 0 and where m χ and E k are understood as m χ /m φ and E k /m φ , respectively. We now explain how to obtain Ω k (λ, θ). We start by subtracting terms that only contribute a time-independent global phase to β k , writing eq. (2.9) as where the slashed out terms are the neglected global phase, and the limit t f → ∞ was taken in the second line. Crucially, this neglected phase is scheme independent, which ensures the same about the remainder. Using the expression for t(λ) in eq. (B.12), we reduce this to with the task being to evaluate the integral on the right. Using the decomposition of eq. (3.6) on E k (t), we write the slow and fast components of Ω k (t) as k (λ(t ))e inθ(t ) , (C. 6) respectively. The slow component integral converges due to the m χ term, and the fast component integral can be solved perturbatively as follows. We write Ω f,k as a sum over Ω which can be solved recursively to obtain Ω (n) k as an expansion in powers of λ. This is because β λ ∂ λ Ω (n) k is always suppressed by an extra power of λ relative to Ω (n) k . To obtain the E (n) k components, it is convenient to write in which k was written in terms of E k . When expanding, it is important to not expand the implicit λ dependence of E k , and instead treat it as O(λ 0 ). Using the results of the RS below eq. (B.11), we write which we can now apply to eqs. (C.5) and (C.7) to solve for Ω k (λ, θ). The results are    up to O(λ 2 ) and O(λ 3 ), respectively. We did not include the λ 3 term of Ω s,k as its derivatives are suppressed by extra powers of λ relative to the λ 3 term of Ω f,k . This because only the latter depends on θ, which has an O(λ 0 ) derivative. In addition, the neglected term includes dependence on α 5 and α 6 , which this appendix does not cover.

E Coefficients in the Bogoliubov formulas
The relevant coefficients for the results of section 4 are listed here. The x j coefficients that appear in A