Gravitational Production of Superheavy Dark Matter and Associated Cosmological Signatures

We study the gravitational production of super-Hubble-mass dark matter in the very early universe. We first review the simplest scenario where dark matter is produced mainly during slow roll inflation. Then we move on to consider the cases where dark matter is produced during the transition period between inflation and the subsequent cosmological evolution. The limits of smooth and sudden transitions are studied, respectively. The relic abundances and the cosmological collider signals are calculated.


I. INTRODUCTION
ΛCDM cosmology has succeeded in explaining the existence and structure of the cosmic microwave background (CMB), the large-scale structure (LSS), the abundances of the light elements and the accelerating expansion of the universe. However, the nature of cold dark matter (CDM) [1], including its production mechanism and interactions, is still not known.
The amount of gravitational dark matter production during inflation is well known and it is roughly proportional to H 3 µ 3 e −2πµ , where µ ≡ m 2 X /H 2 − 9/4. In reality, after inflation, there will eventually be the radiation-dominated universe. There are a variety of mechanisms of the so-called preheating/reheating period sandwiched between the inflation and radiation-dominated universe [23,24] (see [25][26][27][28][29] for reviews). When the mass of the gravitationally produced dark matter particle is large, the production rate is in general exponentially suppressed in terms of the mass of the dark matter particle. In this case, the produced superheavy dark matter particle may not be sufficient to explain the existing dark matter relic abundance. However, as noted recently in [30], during an inflaton oscillation regime [31,32], where the scale factor oscillates very rapidly, dark matter particles with number density proportional to H 3 can be produced. An analytical approach for this case is developed recently in [33].
We discuss that there is another scenario where the number density of order H 3 superheavy dark matter can be produced, where there is a sharp transition between the inflation-radiation period. This may be realized in some inflationary models such as quintessential inflation [34], as discussed later.
We will also revisit smooth transition cases by introducing the Stokes line method, which was used in the analysis of the Schwinger effect [35]. In the community of the gravitational dark matter production, the common method is to calculate the Bogoliubov coefficients in the WKB approximation [36][37][38] which assumes that the mode functions of the fields are close to adiabatic mode functions. Quantum transitions such as excitations with slowly changing Hamiltonian [39][40][41] and particle productions in slowly changing background [42] are exponentially small, implying that the desired Bogoliubov coefficients can not be calculated by introducing a normal series solution with integer orders. On the other hand, the characteristics of mode functions change considerably during the universe's expansion, and we have to identify carefully the dominant and subdominant parts in the WKB ansatz [45] in order to interpret the Bogoliubov coefficients as particle productions. Some researchers [35,[43][44][45][46] considered these questions and argued that the evaluations of particle productions require more information besides WKB ansatz, i.e. the Stokes phenomenon which involves the emergence of subdominant component with negative frequency [47]. These studies imply that the Stokes phenomenon can indicate the production events and provides results with a reasonable precision, and therefore we may adopt this interpretation. On the other hand, the ref. [48] showed that the superadiabatic approximation can describe universal and smooth particle productions when the time evolution hits the Stokes line. This method utilizes Dingle's theory of asymptotic series [49]. We adopt this method to calculate the gravitational dark matter production in a smoothly changing background.
There are many models for gravitationally produced dark matter with different mass ranges and different types of interactions [50][51][52][53][54][55][56][57]. One may ask if it is possible to determine some properties of such dark matter independently of the details of these models. One such possibility is to make use of the cosmological-collider signals [58][59][60][61][62]. By measuring the squeezed-limit non-Gaussianity, from the frequency of its oscillation, we can in principle read off the mass of the dark matter particle. The spin information may also be read off from the angular dependence in the squeezed limit. Note though that the frequency depends on the inflationary Hubble scale and possible non-minimal gravitational couplings as well. Also, the coupling between dark matter and the inflaton may affect the relic abundance of dark matter, or correct the mass of the dark matter particle in a significant way, as we will discuss later. This paper is organized as follows: In Section II, we setup the model and study three scenarios of the dark matter gravitational production in the early universe. The dark matter relic abundance is also calculated. In Section III, we discuss associated cosmological collider signals which may help to determine the mass of the gravitationally produced super-heavy dark matter particle. We conclude in Section IV.

II. DARK MATTER AND ITS RELIC ABUNDANCE
Consider a FRW universe with the metric ds 2 = −dt 2 + a 2 (t)dx 2 . We denote the dark matter field as X. The action is where R is the Ricci scalar, and the last term represents the non-minimal coupling.
The equation of motion for X is We can write the field X in terms of modes as where a k and a † −k are the annihilation and creation operators that satisfy the commutation relations [a k , a k ] = 0 and [a k , a † k ] = (2π) 3 δ (3) (k − k ). We find where H is the Hubble parameter here. Note thatḢ is very small in the slow-roll inflation cases and vanishes in the case of exact de Sitter space.
We will consider three scenarios of the gravitational production of superheavy dark matter.
The production in an exact de Sitter phase is reviewed in Section II A. We use the Stokes line method to estimate the particle production in a toy universe in which inflation is connected to Minkowski spacetime in Section II B. Section II C is devoted to an analysis of a sudden transition between the inflation and radiation period.

A. Slow Roll Inflation
In this section, we review the gravitational particle production during slow roll inflation, approximated by a period of de Sitter expansion. This part is well known and recently reviewed in [63]. The particle production in the presence of a background electric field is discussed in [64].
During slow roll inflation, the scale factor is approximately a(t) = e Ht = −1/(Hτ ), and R = 12H 2 . When inflation ends, the scale factor starts to evolve in a non-accelerated way.
There are two classes of solutions to the massive field equation of motion, corresponding to "in" state and "out" state, respectively [58], where H (1) ν (x) is the Hankel function of the first kind, and J ν (x) is the Bessel function. The mode functions are related via a Bogoliubov transformation as Inserting the explicit expressions for the "in" state mode function and "out" state mode function, we obtain the following expressions for the Bogoliubov coefficients then |β k | 2 can be evaluated as Then gives the total number of particles produced per comoving three-volume, from the past infinity to future infinity, as shown below (see also [64]). We can focus on the time duration when |ω k /ω 2 k | is largest. By studying the maximum value of |ω k /ω 2 k |, we know that particles are produced mostly around the time −kτ ∼ µ .
Using this relation to rewrite the k integral into τ integral, we have Then the integral in (10) can be evaluated as The production rate per unit physical four-volume can be obtained as The physical number density of particles at each moment can be evaluated as: which is constant in time. Hence, as a particular case, n X can be regarded as the number density at the end of inflation. Then the current DM abundance is [16]: in the unit of GeV.
Here Ω R and T 0 are the radiation energy fraction and temperature today. This shows that the relic density is exponentially sensitive to µ. If we further assume that reheating happens soon enough after inflation ends, T RH HM pl /g R by energy conservation and µ = 1, we will need H 10 9 GeV in order to make X a major component of DM. Here we assume that production during the exact de Sitter phase is dominant and determines the eventual relic density, but it may be subdominant depending on post-inflationary phases, as we will discuss below.

B. Inflation Connected to Minkowski Spacetime
In this section, we discuss the case where inflation is smoothly connected to Minkowski spacetime. The particle production in this scenario would approximate the particle production in more realistic scenarios where inflation is connected to a stage of the universe with much lower Hubble scale (and thus approximately Minkowski), for example, radiation-dominated universe with transition time of order H −1 .
We introduce the Stokes line method to compute the dark matter relic abundance. The method utilizes the fact that the adiabatic expansion from the WKB approximation is a factorially divergent asymptotic series which can be resummed by the Borel's summation.
Under a proper truncation of the asymptotic series, smooth particle production is obtained while the physical time is crossing the Stokes line of mode functions [48].
We parametrize our toy model of the scale factor evolution as Note that in this subsection, H is the initial Hubble parameter, which is a constant and not the Hubble parameterȧ/a for all times.
For t 0, the scale factor a(t) approaches de Sitter spacetime. For t 0, the scale factor a(t) approaches Minkowski spacetime with a approaching unity. We would like to use the Stokes line method to calculate the particle production in this toy model. This study would shed some light on the particle production in the real universe. We present the main steps in the following and leave a detailed review of this method in Appendix A. Applications of this method to the inflationary universe is also given in this Appendix A.
We use t c to denote the complex time on lower half-plane satisfying ω k (t c ) = 0, and f k (t) can be expressed as where S k is the Stokes multiplier which indicates the moment when the Stokes line is hit, and φ is an unimportant constant phase. Thus, the Bogoliubov coefficients are We The equation describes n X estimated by numerical integration well for a wide range of m X /H, shown on the right panel of Fig. 2.
To be consistent with observations, we require In single-field slow-roll models of inflation, the initial Hubble parameter can be related to the tensor-to-scalar ratio r as [65] H = 8 × 10 13 r/0.1GeV . For r = 0.05, we have m X /H ∼ 4.93 to explain the full dark matter relic abundance.
The corresponding value for r = 0.001 is m X /H ∼ 4.125. We draw the parameter space compatible with Ω DM = Ω X in Figure. 3.

C. Sudden Transition between Inflation and Radiation Domination
In this subsection we aim at evaluating the particle production during a sudden inflationradiation transition, with the transition time scale ∆t H −1 . Such a situation may be realized by modifying the inflaton potential for quintessential inflation [34]. Radiation can also be produced gravitationally, or one can also introduce coupling between the inflaton and other fields to reheat the universe [7,66]. See also [67]. In this case, we can parametrize the evolution history of the universe a(τ ) in the following way, 1 such that a andȧ are continuous. Note that in this subsection, H is a constant parameter (initial Hubble parameter in the inflation stage), which is no longer interpreted as the Hubble parameterȧ/a after the sudden transition.
The dark matter X is quantized as Let us focus on subhorizon modes with k min ≡ am X < k, for which the mass term can be neglected. We also assume ξ = 0 for simplicity. The mode function satisfies the Mukhanov-Sasaki equation For the modes with k aδt −1 ≡ k max , the solutions can be approximated as follows. First, the above equation is solved by the coefficients c 1,2 can be determined by the connecting conditions u k (τ − end ) = u k (τ + end ) and u k (τ − end ) = u k (τ + end ) as follows: From these we find Inserting the expression of c 1 and c 2 into the expression for u k (τ ) when τ > τ end , On the other hand, the mode function u k (τ ) can be written as where u out is the positive-frequency component of the mode function with the normalization condition So the Bogoliubov coefficients are given by which satisfies the normalization condition |α k | 2 − |β k | 2 = 1.
We have to ensure that the produced particle does not have much backreaction on our background described by Eq. (23) and (24). So we can estimate the minimum of ∆t by requiring that the energy density of the produced particles is comparable to the inflationary vacuum energy as follows: That is, as long as ∆t , backreaction is negligible. Now we estimate the particle number produced in this comoving k ranging from k min to k max as Assuming k max No matter what value of µ we take, the de Sitter particle production is much less than the particle production from a sudden inflation-radiation transition era. The corresponding DM relic abundance can be estimated similarly: For instantaneous reheating case, we need m X ∼ H ∼ 10 8 GeV to produce enough DM, which is much smaller compared to the ones in smooth transition scenarios as expected.

III. DARK MATTER ON THE COSMOLOGICAL COLLIDER
Let us consider an effective field theory where the dark matter sector couples to the primordial curvature perturbation ζ in the following way [68]: These two operators can be regarded as coming from the effective field theory of inflation, or two terms coming from the full Lagrangian of the original quasi-single field inflation with a constant turning trajectory [59]. Since if we consider other types of interactions, the procedure of obtaining the characteristic signals will be similar. Thus, we focus on these two terms for simplicity.
Due to the above interactions, the DM leaves its imprints in the primordial non-Gaussianity.
Here let us assume that m X H. The dominant contribution to the non-Gaussianity is from integrating out the dark matter field [69][70][71][72][73]. The resulting non-Gaussianity is of the shape of general single field inflation [74] and of equilateral shape. The magnitude of non-Gaussianity is power-law suppressed in the large-mass limit [75].
There is another type of signal, which, in magnitude, is smaller than the equilateral shape non-Gaussianity from integrating out heavy fields. Usually they are subject to Boltzmann suppression e −πµ in the large mass limit. This signal, however, is more informative and can tell us directly the mass of the dark matter field in terms of Hubble rate during inflation.
This is the cosmological collider signal in the squeezed limit non-Gaussianity from which we can also measure the spin of the dark matter particle. This mechanism is known as the The Feynman diagram that contributes to the cosmological collider signal in the squeezed cosmological collider [62] and is closely related to quasi-single field inflation [58][59][60] and the primordial quantum standard clocks [76][77][78].
Given the form of interactions, the ζ k 1 ζ k 2 ζ k 3 (the prime indicates that we ignore the momentum-conservation factor (2π) 3 δ 3 (k 1 + k 2 + k 3 )) is contributed by the loop diagrams.
The dark matter loop can be attached to any of the k 1 , k 2 and k 3 legs. In the squeezed diagram is initially proposed in [62]. The application to the Standard model mass spectrum can be found in [79][80][81][82], and it had also been applied to complex scalar fields in [83]. The crucial point is that although in general, loop diagrams may suffer from UV or IR divergence [84], and usually we need to introduce local counter terms to cancel the UV divergence, the clock-signal part of the contribution does not suffer from UV or IR divergence. The reason is that the contribution is coming from the non-local process. In this case, we only have to use a double Fourier transformation to deal with the momentum integral. We present the details of the computation in Appendix B. It gives the standard cosmological collider signal in the squeezed limit as where the factor g(µ) is We can express the bispectrum in terms of the dimensionless shape function S(k 1 , k 2 , k 3 ) [85] as where P where k 1 ∼ k 2 k 3 . In this limit, the shape function is and we plot k 1 /k 3 × S(k 1 , k 3 ) in FIG. 6.
The estimator of non-Gaussianity can be defined as [81] f NL 8c 1 c 2 g(µ) The modest lower bound of f NL can be estimated for an effective field theory with Planck scale cutoff. For such a theory, c 1 ∼ H and c 2 ∼ 1 (note that a 1/( √ M p ) factor emerges when we canonically normalize ζ). This will lead to vanishingly small non-Gaussianities. For the effective field theories which are not Planck mass suppressed, c 1 can be larger than H, c 2 can be larger than 1. But too large couplings with the inflaton may change the relic abundance. To estimate the upper bound, and thus cause inconsistencies for the dark matter production.
Note that the two operators in (36) actually originates from the following two actions where c 1 andc 1 , c 2 andc 2 are related viã Now we want to considerc 1 andc 2 as order one constants. The bound on the scale of the EFT should be given by the constraint that there are no overproduction of the dark matter relic abundance. Note that the decay chanel from a single inflaton to dark matter is kinematically forbidden since the inflaton is lighter than the dark matter. Thus, to check the relic abundance, the leading contribution comes from the dark matter produced by collisions of the thermalized inflaton particles.
The explicit form of the constraint depends on the details of reheating. If the inflaton is thermalized rather efficiently during reheating, to get the correct relic abundance for dark matter, we need For the typical parameter space, we have Λ ≥ 10 5 T RH . Depending on the reheating temperature, Λ can take different values.
On the other hand, if the inflaton has no chance to be thermalized, the above bound does not apply. For example, one may consider brane inflation for the sudden transition [86,87], where the inflaton just disappears at the moment of reheating. In addition, for smooth transition, one may consider quintessential inflation, where reheating can also be realized by gravitational particle production, without assuming the coupling between the inflaton and the standard model sector. In principle, dark matter can be created from the inflaton during kination as well, but during kination the energy density of the inflaton gets redshifted rapidly ∝ a −6 , and hence this production would not cause a problem for our estimations. One can also consider a scenario where the inflaton decays to a heavy field, while dark matter is created gravitationally at the end of inflation, and the former heavy field may later decay to radiation to reheat the Universe. In this case the reheating temperature is determined by the decay rate of the heavy field, and then the reheating temperature can be chosen to be sufficiently low so that dark matter cannot get produced at the reheating.
In general, either way, we need to discuss gravitational particle production by a smooth transition from inflation to kination or an early matter domination, and applying the Stokes line method to such situations would also be possible, though we used a smooth transition from inflation to Minkowski spacetime for simplicity. In gravitational reheating one also needs to worry about an overproduction of gravitons, but this issue can be circumvented by spinodal instability of the Higgs or by introducing a sufficient number of degrees of freedom for radiation gravitationally produced [67]. One can also use a scenario of [7]. We estimated the relic abundance assuming smooth transition from inflation to radiation domination, but an extension of that analysis to cases with a kination phase inserted would be straightforward.
Another lower bound on Λ (and thus upper limit on f NL ) can be obtained from the consideration that (42) does not change the dark matter mass too much (otherwise we lose the prediction power for the mass of the dark matter particle on the cosmological collider).
For instance, if ξ ∼ 1 and R ∼ H 2 this means if Λ > 10 4 H, then this additional term is negligible. Otherwise, the mass of dark matter on the cosmological collider will significantly differ from the probes that we expect to use to observe dark matter now. This limit will generically give f NL < 10 −6 . Thus, we expect that if the cosmological collider signal is observed for such dark matter production scenarios, it is very likely that the observed mass of the dark matter is strongly corrected by the inflaton coupling. However, this conclusion depends on the coupling details between the inflaton and the dark matter. Since we haven't exhausted all possible couplings between the inflaton and dark matter, rather only studied the simplest ones, it remains interesting to see if there can be dark matter couplings with the inflaton which can naturally take large values.
The dark matter field may also leave some imprints on the power spectrum through the relation which relates ζζX 2 and ζζ when the comoving momentum of ζ is soft [88].
However, these relations may not show characteristic features of dark matter, but rather are universal for all matter components [89]. For example, one can expand the correlation function in the series of the ratio between the soft leg and the hard leg. The leading order gives the Maldacena consistency relation which is fixed by dilatation symmetry [90]. The next-to-leading order is fixed by special conformal symmetry [91]. It is thus unclear how to extract dark mater properties from these correlators.
Finally, we briefly comment possible isocurvature fluctuations from dark matter. A field like X that is not directly coupled to the inflaton will introduce isocurvature perturbation.
From ref. [50], we can estimate the size of P δX as Since X is a heavy field (µ > 0), X(k) ∝ ( k aH ) iµ , which has no significant k dependence. At large scale, k → 0, P δX will be suppressed by k 3 . In other words, the power of isocurvature mode is diluted by inflation with 40 e-folds. Therefore we do not expect any large scale isocurvature signal to be observed.

IV. CONCLUSION
We considered three cosmological scenarios and calculated the dark matter relic abundance produced gravitationally. Our conclusions are as follows: For exact de Sitter space, the number density is proportional to H 3 µ 3 2π 3e 2πµ −1 , where µ ≡ m 2 X /H 2 − 9/4. For a universe where inflation is connected to a Minkowski universe, we used the Stokes line method to calculate the dark matter relic abundance. The resulting number density is exponentially suppressed. We fit the results numerically by Eq. (20). For a sudden transition where inflation is immediately connected to a radiation-dominated universe, the produced particle number density is proportional to H 3 .
The cosmological collider signal of the dark matter particle is discussed. We note that with the simplest couplings between the inflaton and dark matter, for the inflaton coupling to satisfy two conditions: (1) does not overproduce dark matter and (2) does not correct the dark matter mass significantly, the non-Gaussianity produced is too small to observe in future observations. It is interesting to see if there are mechanisms to boost the non-Gaussianity of the scenario to observational range to test the scenarios of gravitationally produced superheavy dark matter. In this appendix, we review the divergence series method and use it to study the particle production for general FRW backgrounds. We will first summarize the result in A 1. The summarized result will then be derived in A 2. In A 3, as a simple example, we use this method to recover the known results for particle production in de Sitter space.

Summary of the Result
To study the particle production problem, we start from the equation of motion of the massive field, which can be rewritten as This equation, together with the instant Minkowski initial condition, defines the problem of particle production on FRW backgrounds.
As inspired by scattering problems in quantum mechanics, the sign of ω 2 k (t) determines whether the mode function is oscillating or decaying in a "potential barrier", and the moments when ω k (t) = 0 are turning points. For large enough mass m, the turning points can be complex. The phase integral along a particular complex line crossing the real axis is connected to an exponentially small component in the mode function [48].
We denote t c as the complex time located at lower half-plane which makes ω k (t c ) = 0, together with a phase accumulated from t c The contour to carry out the integral will be specified in Eq. (A6).
We define the Dingle's singulant variable With these definitions, we write down the form of the approximate solution of Eq. (A1) [48]: where S k is called the Stokes multiplier function The constant factor exp[−Θ k (t i )] in the first line of Eq. (A4) is for matching the adiabatic ω k dt when S k → 0, and the amplitude of the negative-frequency part is then suppressed by the exponential The moment when ImF k (t) = 0 corresponds to the emergence of the negative-frequency part of the mode function, and the set of complex t satisfying this condition forms the Stokes line [48]. In the above equation, t m is where the Stokes line intersects the real time axis. In the integration, we first integrate from t c to t m along the Stoke line, then integrate from t m to t along the real axis.
From Eq. (A4), the approximated Bogoliubov coefficients can be extracted as which agrees with the results shown in [35,45,46].

The Divergent Asymptotic Series Method
As suggested by Dingle [49] and Berry [48], to get the particle production in (A4), we assume that the dominant part of Eq. (A1) has an asymptotic series solution y in terms of a large parameter m (or using µ = m 2 − 9H 2 /4 in the cases with constant H).
where b (j) k ∝ m −j is the j-th order correction. Substitute the series in Eq. (A1), yielding It is convenient to introduce a variable W k (t) = ω 2 k (t) to investigate the analytic property of the differential equation. With this substitution, Eq. (A9) becomes By rearranging the equation and integration by part, Eq. (A10) generates Using the fact that W k corresponds to O(m 2 ) and matching the orders of m in the series (A8), we obtain a recurrence relation for b Since the series is expanded around t c , it is reasonable to investigate the recurrence relation around this point. Assuming that W k is analytical around t c , we can expand in terms of (t − t c ) n : By defining a variable Eq. (A12) reduces to a polynomial differential equation where we keep only the dominant term in the integrand. This recurrence relation can be solved by applying an ansatz b This solution is quite complicated, but it is enough to know the divergent behavior when j → ∞. We rewrite Eq. (A17) in terms of the singulant variable in the large-j limit as b (j) One can easily check that b which agrees with the ansatz of asymptotic series (A8), and this equation indicates that the series B k is divergent in factorial form. We will show that the factorial divergence is crucial since it is related to the Borel summation which generates a subdominant term in the solution.
Now we apply Berry's theory [48] to derive the Stokes multiplier function S k (t) shown in Eq. (A5) which describes the moment when massive particles emerge. To handle a divergent asymptotic series such as Eq. (A8), the standard method [47] is to truncate the series into a finite sum and a divergent tail as assuming n is large enough so that Eq. (A18) is applicable. By comparing with Eq. (A4), the terms inside the square bracket is the Stokes multiplier, and we then apply the Borel summation to make the series sum meaningful. The first step is to convert the factorial into the integral of Gamma function, and we denote the multiplier as where the variable changes to z = s F k (t) − 1 in the last line, and the contour is deformed that the upper limit is +∞. The denominator z indicates that the magnitude of the integrand dominates at z ≈ 0, and we can choose n ≈ |F k (t)| + 1 such that the phase is stationary at z ≈ 0, implying that the integral can be evaluated with saddle-point approximation.
Before proceeding the calculation of S k (t), we first justify two conditions required by the validity of the saddle-point approximation: |F k (t)| has to be sufficiently close to |ReF k (t)| since adjusting n can only cancel the real part of the exponent, and |ReF k (t)| has to be large enough such that Eq. (A18) is applicable. The first condition is valid when the evolution is close to the Stokes line when ImF k (t) = 0, implying that the saddle point approximation is valid around the particle production, and it is sufficient to explain the emergence of the sub-dominant term in the mode function. To justify the second condition, we use the fact that ReF k (t) is a constant on real axis and proportional to the large parameter m. For the cases we considered in this paper, ReF k (t) is sufficiently large as long as m H.
With the notations F R k = ReF k and F I k = ImF k (t), we expand the integrand in Eq. (A20) around z = 0 and approximate the integral from −∞ to +∞: where the constant and we apply the fact in the last line that |n − 1 − F R k | O(1) as long as n is sufficiently large. Physically, the initial condition requires S k (t) = 0 when t is small, which indicates that the correct S − should be 1/2. In the original method [48] this is handled by taking the principal value of the integral. Here we show that one can also evaluate S − by deforming the contour near the origin to be an infinitesimally small semicircle below the singularity z = 0.
Since the integrand is odd, only the integral over the semicircular contour contributes to where Θ denotes the Heaviside step function. We can analogously define four types of propagators for ζ, G ++ , G +− , G −+ , G −− in terms of u k (τ ).

(B18)
Since p and q are constrained by momentum conservation p + q − k 3 = 0, First we write the delta function δ (3) (p + q − k 3 ) into an integration of exponential function in x space and integrate out q and p by Fourier transform. To write this procedure more explicitly, we have and the Fourier transform is given by Then, we integrate out x. After this procedure, we obtain where we have used the fact that u k (0) = u * k (0). Finally, we obtain the non-Gaussianity as where the factor g(µ) is given by