On Thermal Gravitational Contribution to Particle Production and Dark Matter

We investigate the particle production from thermal gravitational annihilation in the very early universe, which is an important contribution for particles that might not be in thermal equilibrium or/and only have gravitational interaction, such as dark matter (DM). For particles with spin 0, 1/2 and 1 we calculate the relevant cross sections through gravitational annihilation and give the analytic formulas with full mass-dependent terms. We find that DM with mass between TeV and $10^{16}$GeV could have the relic abundance that fits the observation, with small dependence on its spin. We also discuss the effects of gravitational annihilation from inflatons. Interestingly, contributions from inflatons could be dominant and have the same power dependence on Hubble parameter of inflation as that from vacuum fluctuation. Also, fermion production from inflatons, in comparison to boson, is suppressed by its mass due to helicity selection.


I. INTRODUCTION
The accumulated firm evidence for dark matter (DM) has challenged modern particle physics for decades. From the galactic rotation curves to galaxy cluster, large scale structure (LSS) and cosmic microwave background (CMB), the existence of DM has been well-established, based only on the gravitational interaction. Numerous models for DM has also been proposed, see Refs. [1,2] for reviews. Broadly speaking, for DM as elementary particles, it either can be in thermal equilibrium with other particles and then freeze out, or was never in equilibrium but still produced gradually through various processes. The first class is usually referred as weakly-interacting massive particle (WIMP), while the second includes axion, sterile neutrino, gravitino and so on.
All the mentioned DM candidates above inevitably have interactions other than gravitation, therefore in principle could give rise to possible signatures in direct, indirect and collider searches. However, so far there is no confirmed evidence in all those searches for DM's new interaction, it is fair to ask what if DM only has gravitational interaction. Recent studies [3,4] have shown that it is viable to generate scalar DM abundantly with only gravitational annihilation, namely particles in the thermal both annihilate into DM through a graviton, the quantum of Einstein's gravity in the weak-field limit. Scenarios and phenomenologies in extended theories are also discussed, for example in Refs. [5][6][7][8][9][10].
In the view of effective field theories, microscopically gravity can be treated effectively as quantum field theory as long as the energy scale is much lower than Planck scale (M P = 1.12 × 10 19 GeV) [11,12]. Recently, it has also been shown that general relativity can be derived as an effective field theory of gravitational quantum field theory with spin and scaling gauge symmetries [13,14]. Since the energy scale during/after inflation has already been constrained to be 10 16 GeV which is much lower than M P , we would expect that the local scattering and/or annihilation through graviton can be described in effective field theory. Then these processes should in principle contribute to the cosmological evolution of all particle species, including DM. Particles with interactions much stronger than gravity would be in thermal equilibrium with other particles and short-range gravitational processes are essentially irrelevant for them. However, if DM is very weakly interacting and was never in equilibrium in the early universe, we should include the contributions from gravitational processes.
In this paper, we investigate the viable mass range for DM with spin 0, 1/2 and 1, produced by the gravitational annihilation of particles in the thermal bath with various spins. We compute all the possible, general annihilation cross sections analytically, including all the finite mass term. We find that for the production from particles in the thermal bath the abundance of DM is tightly related with the highest temperature T max after inflation, proportional to T 3 max /M 3 P if its mass m X < T max and m 3 We also discuss the effects from inflation dynamics and show that, gravitational annihilation from inflatons might be the dominant channel for scalar/vector DM production (there is a suppression factor for fermionic DM due to helicity selection) and interestingly has the same power dependence on Hubble parameter as production from vacuum fluctuation. This paper is organized as follows. In Sec. I we start with the standard Boltzmann equation to follow the cosmological evolution of particles and establish the convention and terminology for later discussions. Then in Sec. III we calculate the gravitational annihilation cross section for different initial and final states with spin 0, 1/2 and 1. Later in Sec. IV we apply our calculated cross section to DM and investigate the viable mass range. In Sec. V we discuss the effects from chaotic inflation and show that inflaton's contribution can be very important. Finally, we give the summary.

II. BOLTZMANN EQUATION
To be self-contained, let us start with the standard Boltzmann equation in cosmology [15] for the evolution of number density n 3 through the 2 ↔ 2 process 1 , p 1 + p 2 ↔ p 3 + p 4 , where a is the scalar factor, Hubble parameter H =ȧ/a, p i denote the spatial momenta, p i for 4-vector, M is the matrix element, f i is the distribution for particle i without internal degree of freedom, +(−) sign in ± is for bosons (fermions) and pol means the sum of all polarizations. For particles that were in thermal equilibrium, such as WIMP, we need to keep both terms in the bracket of Eq. (2.1). This is due to the cross symmetry M 12→34 = M 34→12 and f 1 f 2 is compatible to f 3 f 4 for E i ∼ m 3 where m 3 is the mass for particle 3. In cases where f 3,4 is much smaller than 1 and/or f 1,2 , we can neglect the second term and the above Boltzmann equation becomes The term in the bracket can be replaced by 4Fg 1 g 2 σ 12→34 , where g i is the spin degree of freedom, σ ≡ σ 12→34 is the cross section and Changing to the integration variables E 1 , E 2 and s, we have where As will be shown in next section, throughout our discussion, we have m 1 = m 2 = m and m 3 = m 4 = M and the integration range then can be simplified to So far, the discussions have been quite general and apply for other very weakly interacting particles as well, see Ref. [16] for a recent review. It is evident that the key part is to calculate the annihilation cross section σ. After that we can perform numerical integration or analytic computation for some special cases. If f 1,2 have quantum statistical distributions, like Fermi-Dirac or Bose-Einstein distributions (e E/T ± 1) −1 , no compact analytic formulas can be derived. However, for E > T , we can use approximate Maxwell-Boltzmann distribution, e −E/T , and then integrate over E − and E + to get where K i is the modified Bessel function of the second kind with order i.

III. ANNIHILATION CROSS SECTION
In this section, we compute the annihilation cross section in the center-of-mass (CM) frame for various initial and final states in Fig. 1. Note that the initial particles are different from the final ones so that the process can change the number density and contribute to Boltzmann equation, although in a broader context for other physics problems they can be the same. Since the cross section is a Lorentz-invariant quantity, the results derived here will also be valid in other frames.
In effective field theory, the leading interactions between graviton and matter are described by where κ = √ 32πG(G is the Newton's constant), h µν is the graviton field and T µν is the energy-momentum tensor for matter fields. This term is linear on h µν but sufficient for our discussions in which less than 2 gravitons appear in the processes. We shall use the harmonic gauge fixing condition for gravity so that the graviton's propagator with momentum p has the following form where η µν is the metric for flat spacetime. Since we are considering the leading-order tree-level scattering processes, we do not need to include higher-dimensional operators with more graviton h µν , renormalization effects and ghost. The symmetric energy-momentum tensors T µν for complex scalar S, spin-1 2 Dirac fermion F , massive vector V and massless vector γ are listed in the following T µν for real scalar φ can easily be obtained by substituting Then we can get the Feynman rules to do the calculation of Fig. 1. To make the results as compact as possible, we extract the common factor for unpolarized collisions where g i is the degrees of freedom for initial state i, S is the symmetric factor (S = 2 for identical final states, for example, real scalars, neutral gauge bosons, otherwise S = 1), | p i | and | p f | are the lengths of three-momentum for initial and final states, respectively. As shown, M is the matrix element and we have defined A as the integration of polarization-summed |M| 2 over the scattering angle θ, with the κ 4 factor pulled out. Note that the kinematic variables in CM frame for m i = mī ≡ m and m f = mf ≡ M , and After some tedious calculations, we obtain A for different processes of initial states with mass m and final ones with M where both the initial and final states can be complex scalar S, fermion F (spin 1/2), massive vector V and massless vector γ. For processes involving final scalar S,

17)
A (γ → γ) = s 2 10 . (3.18) Note that we can use the cross symmetry, A (f → i) = A (i → f ) with interchanging m ↔ M , to get As for other processes, such as A (S → F ), A (S/F → V ) and A (S/F/V → γ). In the case s 4m 2 and s 4M 2 , we can neglect the mass-dependent terms and get very concise As which are just proportional to s 2 .
The above results have shown consistencies under several checks. For example, A is gauge invariant when involving massless vector γ where we have explicitly checked in R ξ gauge and the results are independent of gauge-fixing parameter ξ in the T µν γ (ξ), And the coefficient of s 2 term in A (V → S) is three times as that in A (S → S), which is due to three polarizations of V . Furthermore, A (i → f )s are symmetric over m and M when the initial and final states are the same, i = f . For later convenience, we also tabulate the case m 2 s in Table. I. One can easily check that A (V → f ) = A (S → f ) + A (γ → f ). Interestingly, we notice that A (γ → f ) = 4A (F → f ) which might be related with spin structures in gravitational interactions.

IV. APPLICATION TO DARK MATTER
With the cross section in hand, we now proceed to compute the abundance for stable particles like DM X with mass M = m X . In the absence of entropy production, we have d(a 3 s)/dt = 0, where s is the entropy density. Therefore, we can rewrite the equation Eq. 2.3 in terms of a more convenient quantity, th yield Y X ≡ n X /s, In the radiation dominant era, we have the following relations, where g * is the total number of effectively massless degrees of freedom. Integrate over temperature from the minimal value to maximum one, we finally get The above result has negligible dependence on T min , so we can freely take T min as zero or the present temperature of CMB. The yield Y X is related with the observed energy fraction for DM Ω X at present time, where Ω b is the energy density fractions of baryon, m p 1GeV is proton mass, n γ is the number density of photon today, s 0 is the total entropy density of photon and neutrino, and η 6×10 −10 is baryon-to-photon ratio. Assuming a minimal particle content in thermal bath, namely only SM, and the temperature is higher than electroweak symmetry breaking, all SM particles are therefore massless and we can use the results in Table. I. Since our formalism is for Dirac fermions, there is a factor of 1/2 for Weyl particles (neutrino in SM). Taking all these into account, we have 2 complex scalars, 45/2(1/2 × 3 + 1 × 3 + 2 × 3 × 3) Dirac fermions and 12(8 + 3 + 1) massless gauge bosons. In Fig. 2, we illustrate the correlation between temperature T max and DM mass m X for fixed Ω X 0.258 [17] in several cases, by integrating Eq. (4.2) numerically. DM with spin zero, 1/2 and 1 are shown with solid, dashed and dot-dashed curves, respectively. The gray dotted line indicates m X = T max , while its left (right) side marks m X < T max (m X > T max ). The turnover of these curves at m X = T max is due to the following reasons. When m X < T max , Ω X is proportional to κ 3 T 3 max and increasing m X would require smaller T max for fixed Ω X , which is A for the case s = 4m 2 (initial states with mass m and final states with mass M ). As usual, m is the mass of initial particle and it is equal to zero for γ. For scalar as final state in the second row, we also include the non-mimimal coupling, ζRS † S.
exactly the reason why we see straight lines in the logarithmic plot. This temperature dependence is different from gravitino production [18] that depends on T max linearly because gravitino is mostly produced by the scattering of supersymmetric particles, not by the diagram in Fig. 1. For m X > T max , the production is exponentially suppressed due to the Boltzmann distribution at high energy, roughly scaling as κ 3 m 3 X exp [−2m X /T max ]. In such a circumstance, Ω X = 0.258 would require T max ∼ m X /10. In case T max 10 12 GeV the produced X is negligible, we shall see in next section that inflation could then play an important role.

V. EFFECTS OF INFLATION DYNAMICS
It is widely believed that in the very early universe there was an exponential expansion called inflation. After inflation, there was a short matter-dominant time as the inflation field φ oscillates around the potential minimum. Then inflatons decay into radiation with decay width Γ φ and reheat the universe with a temperature T R ∼ Γ φ M P . In the simplest approximation, we may just take T max = T R in our above discussion. However, realistically the effects from inflation are model-dependent since different inflationary scenarios could give various cosmological evolutions. More importantly, inflatons can also annihilate gravitationally into other particles and contribute the production. Here, we discuss some possible effects and only focus on the simplest chaotic inflation 2 , for example, with quadratic potential, although our formalism might also apply for other cases.
For inflation field φ with canonical kinetic term and general potential V (φ), its energy-momentum tensor is given by During inflation, we can use homogeneous field configuration φ(t) for the background evolution and get the energy and pressure densities Using the equation of motion for φ [15]φ where Γ φ is the decay width that closely connects to the reheating that inflatons decay into other particles and reheat the Universe, we can obtain the evolution equation for ρ φ , Usually, averaging over several oscillations is performed so that one can use Virial theorem to replaceφ 2 with averaged ρ φ , andρ φ just follows the evolution equation for non-relativistic matter. However, we do not perform such an average as we shall see immediately that the dominant production from inflaton happens at the transition time, not at the oscillation time. To compute the particle production from inflaton annihilation, we treat inflatons as particles with zero spatial momentum, namely the distribution of φ particle is where m φ is the mass of inflaton. We are aware that inflation field can not be always treated as collection of inflaton particles, for example, if some particle couples to inflaton non-gravitationally, its production should be calculated by solving the equation of motion and regarding inflation field as classical background [20]. Since in our case there is no direct coupling between other particles and inflaton, we simply use Eq. (5.5) for our estimations. This reduces Eq. (2.3) to a very compact formula Note that generally Fσ at s = 4m 2 φ does not vanish because the factors (s − 4m 2 φ ) in both F and |p i | of σ in Eq. 3.7 cancel with each other, unless there might be helicity/parity selection rules for different initial and final states so that the integrated squared matrix elements A is identically zero. We can see several examples in the second row of Table. II with scalar as the initial state. For instance, for conformal coupled massless scalar (ζ = 1/6), massless fermions and vectors, A vanishes, which means that they are not produced by inflaton's gravitational annihilation during inflation.
For massive scalars and vectors, we have Fσ κ 4 m 4 φ /256π and can roughly estimate how much particles are produced at inflation during a Hubble time interval, ∆t ∼ 1/H * , where H * is the Hubble parameter around the transition era between inflation and oscillation time. Interestingly, the above formula has the same power dependence on Hubble parameter as particle creation from vacuum fluctuations during inflation [20][21][22][23][24][25] and oscillation [26,27]. Moreover, the feature that conformal coupled massless scalars with ζ = 1/6, massless fermions and vectors are not produced by inflaton's gravitational annihilation during inflation also agrees with the results for vacuum fluctuations. This might be just a coincidence, or imply some deep underlying connection between two mechanisms, which however is beyond our scope here. In some sense the calculation with non-minimal coupling ζ serves as additional check of our computation.
To get the yield, take H * = m φ and assume instantaneous reheating, we have T max T R = Γ φ M P and From the above estimations, we can also learn that gravitational annihilation from inflatons might be dominant over the contributions from thermalized particles after reheating since the later one goes like T 3 R /M 3 P . From Fig. 2, it is obvious that for T max 10 12 GeV the thermally produced X is negligible. On the other hand, annihilation from inflations could still produce X too abundantly for large m φ , unless m X and/or Γ φ is small enough, for instance, m X 1TeV for Γ φ ∼ 10 −9 m φ 10 −14 M P . However, for massive fermions there would be a suppression factor m 2 X /m 2 φ from the annihilation cross section. This feature is in sharp contrast with the contributions from thermalized particles discussed in previous section where production for particles with different spins are at the same order. In Fig. 3, we show how the contour Ω X = 0.258 goes in the m X -T R plane for spin-0, 1/2 and 1 with m φ = 10 15 GeV and T R 10 12 GeV as an example. One can easily see the dramatic difference discussed just above between spin-1/2 and spin-0/spin-1, which shows that generally spin-0/1 DM would require much lower reheating temperature or lighter mass. We also notice that spin-0 coincides with spin-1 case except in the high mass regime when m X is close to m φ , which can be easily understood from A in Table. II because the longitude mode dominates in the high energy limit.

VI. SUMMARY
We have investigated the particle production from gravitational annihilation of thermal particles in the very early universe. In the case that dark matter (DM) particle might only have gravitational interaction, we have calculated the relic abundance and the possible viable mass range for DM with spin 0, 1/2 and 1. We have computed the 0.258 on the reheating temperature TR and DM mass mX , in spin-0 (solid black line), spin-1 2 (blue dashed line) and spin-1 (purple dot-dashed) cases. We have used m φ 10 15 GeV as example. Spin-0 coincides with spin-1 case except in the high mass regime when mX is close to m φ , which can be easily understood from A in Table. II. analytical cross section for general gravitational annihilation processes through a graviton. DM could be produced by gravitational annihilation of all other particles in the thermal path after inflation or inflatons during inflation. The first contribution crucially depends on the highest temperature T max after inflation, proportional to T 3 max /M 3 P if DM mass m X < T max and m 3 X /M 3 P exp [−2m X /T max ] if m X > T max . Particles with different spins produced by thermal bath are at similar order and can give the correct abundance for DM with mass between TeV and 10 16 GeV. While the second contribution from inflaton depends on the inflation scale, reheating temperature and also the spin of DM (spin 1/2 case is suppressed due to helicity selection, compared with scalar and vector particles). We have shown that in simplest chaotic inflation model the contribution from inflaton's gravitational annihilation could be the dominant production mechanism for stable particles like DM.