Production of Purely Gravitational Dark Matter: The Case of Fermion and Vector Boson

TL;DR

This paper analyzes how purely gravitational interactions can produce fermionic and vector dark matter particles in the early universe, potentially explaining their relic abundance without other interactions.

Contribution

It provides a comprehensive analysis of gravitational particle production for fermionic and vector dark matter, including heavy particles relative to the Hubble parameter.

Findings

Particles lighter than the inflaton can be produced abundantly.

Heavy particles above the Hubble scale can also be generated efficiently.

Relic abundance depends on inflation scale and particle mass.

Abstract

We consider the simplest possibility for a model of particle dark matter in which dark matter has only gravitational interaction with the standard model sector. Even in such a case, it is known that the gravitational particle production in an expanding universe may lead to a correct relic abundance depending on the inflation scale and the mass of dark matter particle. We provide a comprehensive and systematic analysis of the gravitational particle production of fermionic and vectorial dark matter, and emphasize that particles which are much heavier than the Hubble parameter but lighter than inflaton can also be produced abundantly.

Figures

Equations

$$S=\int d^{4}xe\left[-\frac{1}{2}\overline{\psi }\left(e_a^{\mu }{\gamma }^a{D}_{\mu }-m\right)\psi \right],$$

$$D_{\mu }={\partial }_{\mu }+\frac{1}{4}{\omega }_{\mu }^{ab}{\gamma }_{[a}{\gamma }_{b]},$$

$${\omega }_{\mu }^{ab}=2e^{\nu [a}{\partial }_{[\mu }{e_{\nu ]}}^{b]}-e^{\nu [a}{e}^{b]\sigma }{e}_{\mu c}{\partial }_{\nu }{e}_{\sigma }^c.$$

$$d{s}^2=-d{t}^2+a^2(t){\delta }_{ij}d{x}^{i}d{x}^j=a^2(\tau )\left(-d{\tau }^2+{\delta }_{ij}d{x}^{i}d{x}^j\right),$$

$${{\omega }_i}^{j0}={\delta }_i^j\mathcal{H},$$

$$S=\int d\tau d^{3}x\left[-\frac{1}{2}\overline{\tilde{\psi }}\left({\delta }_a^{\mu }{\gamma }^a{\partial }_{\mu }-am\right)\tilde{\psi }\right].$$

$$\tilde{\psi }(\vec{x},\tau )=\int \frac{d^{3}k}{(2\pi {)}^3}{\psi }_{\vec{k}}(\tau )e^{i\vec{k}\cdot \vec{x}}.$$

$$\left({\partial }_{\tau }{\gamma }^0-i\vec{k}\cdot \vec{\gamma }-am\right){\psi }_{\vec{k}}(\tau )=0.$$

$${\psi }_{\vec{k}}(\tau )=\sum _{h=\pm }\left[u_{\vec{k},h}(\tau )b_{\vec{k},h}+v_{\vec{k},h}(\tau )b_{-\vec{k},h}^{†}\right],$$

$$u_{\vec{k},h}^{†}(\tau )u_{\vec{k},h^{\prime }}(\tau )=v_{\vec{k},h}^{†}(\tau )v_{\vec{k},h^{\prime }}(\tau )={\delta }_{h{h}^{\prime }},u_{\vec{k},h}^{†}(\tau )v_{\vec{k},h^{\prime }}(\tau )=0.$$

$$\left\{b_{\vec{k},h},b_{{\vec{k}}^{\prime },h^{\prime }}^{†}\right\}={\left(2\pi \right)}^3\delta (\vec{k}-{\vec{k}}^{\prime }){\delta }_{h{h}^{\prime }},\left\{b_{\vec{k},h},b_{{\vec{k}}^{\prime },h^{\prime }}\right\}=\left\{b_{\vec{k},h}^{†},b_{{\vec{k}}^{\prime },h^{\prime }}^{†}\right\}=0,$$

$$u_{\vec{k},h}(\tau )=\left(\begin{array}{c}{u}_{\vec{k},h}^{+}(\tau )\\ u_{\vec{k},h}^{-}(\tau )\end{array}\right)\otimes {\xi }_{\vec{k},h},$$

$$i{\partial }_{\tau }{u}_{\vec{k},h}^{\pm }(\tau )+hk{u}_{\vec{k},h}^{\mp }(\tau )\mp am{u}_{\vec{k},h}^{\pm }(\tau )=0,$$

$${\partial }_{\tau }^2{u}_{\vec{k},h}^{\pm }(\tau )+\left[{\omega }_k^2(\tau )\pm i(am{)}^{\prime }\right]u_{\vec{k},h}^{\pm }(\tau )=0,{\omega }_k^2(\tau )\equiv k^2+a^2{m}^2.$$

$${\left|u_{\vec{k},h}^{+}(\tau )\right|}^2+{\left|u_{\vec{k},h}^{-}(\tau )\right|}^2=1.$$

$$u_{\vec{k},h}^{\pm }(\tau )=iϵ\sqrt{\frac{-\pi k\tau }{4}}{e}^{\pm \frac{\pi m}{2H_{\inf }}}{H}_{{\nu }_{\pm }}^{(1)}(-k\tau ),{\nu }_{\pm }\equiv \frac{1}{2}\mp \frac{im}{H_{\inf }},$$

$$u_{\vec{k},h}^{+}\to \sqrt{\frac{{\omega }_k+am}{2{\omega }_k}}{e}^{-i\int {\omega }_{k}d\tau },u_{\vec{k},h}^{-}\to -h\sqrt{\frac{{\omega }_k-am}{2{\omega }_k}}{e}^{-i\int {\omega }_{k}d\tau }.$$

$$H_{\nu }^{(1)}(x)\to \sqrt{\frac{2}{\pi x}}{e}^{i(x-(2\nu +1)\pi /4)}\mathrm{for}x\to +\infty .$$

$$H_{\nu }^{(1)}(x)\to -\frac{i}{\pi }\Gamma (\nu ){\left(\frac{2}{x}\right)}^{\nu }\mathrm{for}x\to 0\mathrm{and}\mathrm{Re}(\nu )>0.$$

$$a^4(\tau ){\rho }_{\psi }(\tau )$$

$$f_{\psi }(\vec{k},\tau )$$

$$=\frac{1}{2{\omega }_k}\left[am+2\mathrm{Im}\left(u_{\vec{k},h}^{+\ast }(\tau ){\partial }_{\tau }{u}_{\vec{k},h}^{+}(\tau )\right)\right]+\frac{1}{2}$$

$$a^3(\tau )n_{\psi }(\tau )=2\int \frac{d^{3}k}{(2\pi {)}^3}{f}_{\psi }(k,\tau ).$$

$$u_{\vec{k},h}^{+}(\tau )=A_{k,h}(\tau )g_{+}{e}^{-i{\int }^{\tau }{\omega }_k({\tau }^{\prime })d{\tau }^{\prime }}+B_{k,h}(\tau )g_{-}{e}^{i{\int }^{\tau }{\omega }_k({\tau }^{\prime })d{\tau }^{\prime }},$$

$$A_{k,h}^{\prime }(\tau )=-\frac{g_{-}^{\prime }}{g_{+}}{e}^{2i{\int }^{\tau }{\omega }_k({\tau }^{\prime })d{\tau }^{\prime }}{B}_{k,h}(\tau ),B_{k,h}^{\prime }(\tau )=-\frac{g_{+}^{\prime }}{g_{-}}{e}^{-2i{\int }^{\tau }{\omega }_k({\tau }^{\prime })d{\tau }^{\prime }}{A}_{k,h}(\tau ),$$

$$u_{\vec{k},h}^{-}(\tau )=-h\left[A_{k,h}(\tau )g_{-}{e}^{-i{\int }^{\tau }{\omega }_k({\tau }^{\prime })d{\tau }^{\prime }}-B_{k,h}(\tau )g_{+}{e}^{i{\int }^{\tau }{\omega }_k({\tau }^{\prime })d{\tau }^{\prime }}\right],$$

$$f_{\psi }(\vec{k},\tau )=|B_k(\tau ){|}^2,B_k(\tau )\simeq {\int }^{\tau }d{\tau }^{\prime }\frac{a^{2}Hmk}{2{\omega }_k^2}{e}^{-2i{\int }^{{\tau }^{\prime }}{\omega }_k({\tau }^{\prime \prime })d{\tau }^{\prime \prime }}.$$

$$n_{\psi }^{(\mathrm{slow})}(\tau )\sim \mathcal{A}{a}^{-3}(\tau )k_c^3,\frac{k_c}{a_{\mathrm{end}}}\equiv m{\left(\frac{H_{\inf }}{m}\right)}^{\frac{2}{3(1+w)}},$$

$$n_{\psi }^{(\mathrm{fast})}(\tau )\sim \mathcal{C}{H}_{\inf }^3{\left(\frac{m}{m_{\inf }}\right)}^2{\left(\frac{a_{\mathrm{end}}}{a(\tau )}\right)}^3,$$

$$\frac{{\rho }_{\psi }}{s}\simeq \frac{m{H}_{\inf }{T}_{\mathrm{R}}}{4M_P^2}\left[\mathcal{C}{\left(\frac{m}{m_{\inf }}\right)}^2+\eta \frac{m}{H_{\inf }}\right],$$

$$\eta =\mathcal{A}\times \{\begin{array}{ll}1& \mathrm{for}m\lesssim H_{\inf }\\ (m/H_{\inf }{)}^2{e}^{-cm/H_{\inf }}& \mathrm{for}m\gtrsim H_{\inf }\end{array}.$$

Full text

1. 1 Introduction

2. 2 Fermion production

3. 2.1 Fermion action in the FRW Universe

4. 2.2 Fermion production

5. 3 Vector boson production

6. 3.1 Vector boson action in the FRW Universe

7. 3.2 Transverse mode production

8. 3.3 Longitudinal mode production

9. 3.3.1 Heavy vector boson case

10. 3.3.2 Light vector boson case

11. 3.4 Combined Results

12. 4 Conclusions and discussion

13. A Conventions

14. B Energy density of vector boson

15. B.1 Energy-momentum tensor and energy density

16. B.2 Energy scaling of the vector boson

17. B.2.1 Transverse mode

18. B.2.2 Longitudinal mode

19. C Thermal production

1 Introduction

Gravitational particle production is a universal phenomenon which happens in the expanding universe [1, 2]. If a field is not conformal, it feels the time-dependent background evolution of the metric and results in finite number density even if we start with the adiabatic vacuum with no particle excitation initially. Gravitational particle production in inflationary universe was discussed in Refs. [3, 4, 5], where it was shown that the particle production happens during the transition from inflationary (or de Sitter) era to the radiation-dominated or matter-dominated universe. Depending on the mass of the particle, these gravitationally produced particles can be dominant component of dark matter (DM). A typical number density of produced particles during the transition is given by $n\sim H_{\rm inf}^{3}$ where $H_{\rm inf}$ denotes the Hubble parameter during inflation for a scalar field with a minimal coupling to gravity if $m\lesssim H_{\rm inf}$ where $m$ denotes the scalar field mass.#1#1#1 Precisely, $H_{\rm inf}$ should be regarded as the Hubble scale at the end of inflation. In this paper we take the Hubble parameter $H$ to be constant during inflation, which is a good approximation for our purpose.

Note also that a light scalar field with $m\lesssim H_{\rm inf}$ gravitationally develops condensation during inflation, whose subsequent oscillation can account for the DM in the present universe (see e.g. Refs. [7, 6] for recent works on general spectator fields).

Recently it was argued that the gravitational particle production also happens during the reheating era in which the inflaton oscillates coherently, due to the fact that the cosmic scale factor $a(t)$ contains a small but rapidly oscillating part [8, 9]. Since the time scale of the coherent oscillation is given by $\sim m_{\rm inf}^{-1}$, where $m_{\rm inf}$ denotes the inflaton mass around its potential minimum, particles with mass smaller than $m_{\rm inf}$ are likely to be produced, along the line of preheating in models where the inflaton directly couples to other fields [10, 11, 12, 13]. The number density produced in this way is also given by $n\sim H_{\rm inf}^{3}$ for a minimally coupled scalar. It is comparable to the contribution discussed above, but now the particle mass can be as large as the inflaton mass $m_{\rm inf}$. Since many inflation models have $m_{\rm inf}\gg H_{\rm inf}$, it opens up a possibility that the super-Hubble mass particles are efficiently produced and become a dominant component of DM. This scenario was further studied in detail in Refs. [14, 15]. See also Refs. [16, 17] for some related works.#2#2#2 Gravitationally interacting particles are also produced by scatterings of the Standard Model (SM) particles in thermal bath through the graviton exchange [18, 19, 20, 21].

This class of scenario is interesting since we do not need any particular interaction of DM with SM particles to account for the present DM abundance. In this sense, it may be the simplest model of DM.

In this paper we consider production of a massive fermion and vector boson which only have gravitational interaction. A fermion or vector boson may be more suitable as a candidate of purely gravitational DM since their interactions with the SM particles can be naturally forbidden at the renormalizable level. For a singlet fermion $\psi$, the only possible renormalizable interaction in the Lagrangian is $\mathcal{L}\sim y_{i}\psi(L_{i}\Phi)+{\rm h.c.}$ where $L_{i}$ $(i=1,2,3)$ denotes the lepton doublet and $\Phi$ is the SM Higgs boson with $y_{i}$ being coupling constants. This coupling is forbidden by, e.g., assuming that the SM fermions are charged under $B-L$ gauge symmetry but $\psi$ is a singlet. Especially, $\psi$ is absolutely stable if the $Z_{2}$ subgroup of U(1) $B-L$ symmetry remains unbroken. For an Abelian hidden vector boson $A_{\mu}$, the only possible renormalizable interaction is the kinetic mixing with hypercharge photon, which may also be forbidden by imposing discrete symmetry. For a scalar field $\phi$, on the other hand, one can always introduce a Higgs-portal coupling $\mathcal{L}\sim-\lambda\phi^{2}|\Phi|^{2}$ which drives $\phi$ into thermal bath unless the coupling constant $\lambda$ is very small. Once thermalized, the DM abundance is determined by the standard freezeout scenario [22, 23]. Thus fermion or vector boson can be more likely to be purely gravitational DM.

The gravitational production of fermion or vector boson is qualitatively different from the case of scalar field with minimal coupling to gravity. This is because the fermion and (transverse) vector boson are conformal in the massless limit. One can choose a canonical basis such that the dependence on the cosmic scale factor completely disappears. Thus there is no particle production in the massless limit, which is similar to the case of a scalar field with conformal coupling. The case of massive vector boson is a bit complicated since the longitudinal component behaves rather like a minimally coupled scalar field and requires a careful investigation.

In Sec. 2 we study the gravitational fermion production. In Sec. 3 we investigate the gravitational production of massive vector boson in detail. In particular, we carefully distinguish the transverse and longitudinal mode and discuss how they behave under the time-dependent background. We conclude in Sec. 4.

2 Fermion production

2.1 Fermion action in the FRW Universe

Let us consider an action of free Majorana fermion $\psi$, which satisfies the Majorana condition $\psi=-C^{-1}\overline{\psi}^{T}$ with $C$ denoting the charge conjugation matrix,

[TABLE]

where $e^{\mu}_{a}$ denotes the vierbein with $a,b,\dots$ and $\mu,\nu,\dots$ represent local Lorentz and general coordinate indices respectively, and $e\equiv\det(e^{\mu}_{a})=\sqrt{-g}$. The covariant derivative is given by

[TABLE]

where the spin connection is defined as

[TABLE]

In the Friedmann-Robertson-Walker (FRW) background,

[TABLE]

the only non-zero components are

[TABLE]

where $\mathcal{H}=a^{\prime}/a$ is the conformal Hubble parameter, which is related to the Hubble parameter $H=\dot{a}/a$ through $\mathcal{H}=aH$. Here and in what follows, the prime (dot) denotes the derivative with respect to conformal time $\tau$ (physical time $t$). It is convenient to perform the rescaling as $\widetilde{\psi}\equiv a^{3/2}\psi$ so that the action becomes

[TABLE]

It is seen that the rescaled field has a canonical kinetic term and the action is independent of the scale factor $a$ in the massless limit $m\to 0$. In other words, a fermion is conformal in the massless limit. Therefore, the rate of gravitational particle production is suppressed by the fermion mass $m_{\psi}$.

It is convenient to work with the Fourier mode since we are interested in the free fermion:

[TABLE]

The equation of motion for the mode function is given by

[TABLE]

Now let us expand the mode function as

[TABLE]

where $v_{\vec{k},h}=-C^{-1}\overline{u}^{T}_{-\vec{k},h}$ and $h$ denotes the helicity degree of freedom. The normalization condition is taken as follows:

[TABLE]

The creation and annihilation operators satisfy the following anti-commutation relation:

[TABLE]

so that the original field satisfies the anti-commutation relation $\left\{\widetilde{\psi}(\tau,\vec{x}),\widetilde{\psi}^{\dagger}(\tau,\vec{y})\right\}=\delta(\vec{x}-\vec{y})$. Let us write the mode function as

[TABLE]

where $\xi_{\vec{k},h}$ denotes the eigenvector of the helicity, which satisfies $(\vec{\sigma}\cdot\hat{\vec{k}})\xi_{\vec{k},h}=h\xi_{\vec{k},h}$ with $\hat{k}\equiv\vec{k}/|\vec{k}|$ and $h=\pm 1$. Taking $\vec{k}$ to be the $z$-direction, we have $\xi_{\vec{k},+}=(1,0)^{T}$ and $\xi_{\vec{k},-}=(0,1)^{T}$. Adopting the Dirac representation for the gamma matrices, the equation of motion becomes

[TABLE]

which may be cast into the second order form,

[TABLE]

Note that $u_{\vec{k},h}^{+}$ and $u_{\vec{k},h}^{-}$ are not independent of each other and the normalization condition implies

[TABLE]

During the de Sitter phase in which the Hubble parameter is given by $H=H_{\rm inf}={\rm const}$, noting $(am)^{\prime}=a^{2}H_{\rm inf}m$ and $\tau=-1/(aH_{\rm inf})$, we find an exact solution to Eq. (14) as

[TABLE]

where $\epsilon=1$ for $u_{\vec{k},h}^{+}$ and $\epsilon=-h$ for $u_{\vec{k},h}^{-}$, and $H_{\nu}^{(1)}(x)$ is the Hankel function of the first kind with order $\nu$. In the far past $(k\tau\to-\infty)$, or the short wavelength limit, it approaches to the same mode function as the Minkowski vacuum:

[TABLE]

Here $a(\tau)\to 0$ limit should be understood. It is also evident that there is no significant growth in the superhorizon limit $(k\tau\to 0)$ during inflation.#3#3#3 Asymptotic form of the Hankel function in the short wavelength limit is

$\displaystyle H_{\nu}^{(1)}(x)\to\sqrt{\frac{2}{\pi x}}e^{i(x-(2\nu+1)\pi/4)}~{}~{}~{}{\rm for}~{}~{}~{}x\to+\infty.$

(18)

In the long wavelength limit, the Hankel function becomes

$\displaystyle H_{\nu}^{(1)}(x)\to-\frac{i}{\pi}\Gamma(\nu)\left(\frac{2}{x}\right)^{\nu}~{}~{}{\rm for}~{}~{}x\to 0~{}~{}{\rm and}~{}~{}{\rm Re(\nu)}>0.$

(19)

2.2 Fermion production

Now let us estimate the gravitational fermion production. The fermion production in the rapidly oscillating background or the fermionic preheating was studied in Refs. [24, 25, 26, 27]. The gravitational fermion production in the expanding universe was studied in Refs. [28, 29, 30] and also the gravitino preheating was extensively studied in Refs. [31, 32, 33, 34, 35, 36, 37]. Here we combine these knowledges to estimate the rate of gravitational fermion production, especially pointing out the contribution from the inflaton coherent oscillation.

The energy density of the fermion is given by

[TABLE]

where the prefactor $2$ counts the two helicity modes and

[TABLE]

denotes the occupation number or the phase space density. The last factor $+1/2$ cancels the negatively divergent energy density due to the fermionic zero-point fluctuations. One sees that $f_{\psi}=0$ for the Minkowski mode function (17). In the time-dependent background $(a^{\prime}>0)$ the mode function may be modified from this asymptotic form and hence we will obtain $f_{\psi}>0$ that signals particle production. The number density is also a useful quantity, which is then given by

[TABLE]

In order to estimate the particle production, we conveniently rewrite the mode function as

[TABLE]

where coefficients are assumed to satisfy

[TABLE]

where $g_{\pm}\equiv\sqrt{(\omega_{k}\pm am)/(2\omega_{k})}$. One can check that $u_{k,h}^{+}(\tau)$ satisfies the equation of motion (14). The other mode function is given by

[TABLE]

and hence the normalization condition (15) implies $|A_{k,h}(\tau)|^{2}+|B_{k,h}(\tau)|^{2}=1$, which ensures that the phase space density cannot exceed unity as expected from the Pauli exclusion principle. The initial condition (17) is equivalent to $A_{k,h}(\tau\to-\infty)=1$ and $B_{k,h}(\tau\to-\infty)=0$. Since the initial condition and the time evolution (25) do not explicitly depend on $h$, $A_{k,h}$ and $B_{k,h}$ are the same for both $h=\pm$. In what follows we omit the helicity subscript $h$ for this reason. Deviation from these initial values $A_{k}=1$ and $B_{k}=0$ indicate particle production. With these definitions, we have

[TABLE]

Thus what we have to do is to calculate $B_{k}(\tau)$ under the background evolution of the cosmic scale factor $a(\tau)$. The calculation is almost parallel to the case of scalar field with conformal coupling. Thus we do not write the detail of the calculation below. Readers are referred to Ref. [14] for more detailed discussion to evaluate this integral. There it was shown that there are two contributions to the particle production: one comes from the “slow” change of the background whose time scale is just parametrized by the Hubble expansion rate $H$ and the other comes from the “fast” change of the background related to the inflaton coherent oscillation whose time scale is characterized by the inflaton mass $m_{\phi}$ [8, 9].

For the “slow” contribution, the production rate is suppressed as $e^{-cm/H_{\rm inf}}$ with $c$ being order one numerical constant for $m\gtrsim H_{\rm inf}$. Thus we focus on the case $m\lesssim H_{\rm inf}$. In this case, the dominant production comes from the epoch around $k\sim am$, i.e., the fermion becomes non-relativistic. The number density produced in this way is estimated as

[TABLE]

where $\mathcal{A}\sim\mathcal{O}(10^{-3})$ is a numerical coefficient and $k_{c}$ is defined as a momentum such that $k/a(\tau)=m$ at $H=m$, $a_{\rm end}$ denotes the scale factor at the end of inflation and $w$ denotes the equation of state parameter after inflation until the epoch $H=m$. Hereafter we take $w=0$ for simplicity, having in mind a scenario that the inflaton coherent oscillation behaves as non-relativistic matter and it finally decays into radiation at $H=H_{\rm R}$ with the reheating temperature $T_{\rm R}=(90/\pi^{2}g_{*})^{1/4}\sqrt{H_{\rm R}M_{P}}$.#4#4#4 We implicitly assumed that reheating is competed after $H=m$, i.e., $H_{\rm R}<m$. Otherwise, $k_{c}$ should be replaced by $k_{c}/a_{\rm end}=m(H_{\rm inf}/m)^{2/3}(m/H_{\rm R})^{1/6}$. This result is consistent with Ref. [30]. The “fast” contribution is given by

[TABLE]

where $m_{\rm inf}$ denotes the inflaton mass and $\mathcal{C}\sim 10^{-3}-10^{-2}$ is a numerical constant. As emphasized in Refs. [8, 9, 14, 15], this fast contribution is not suppressed even for $m\gg H_{\rm inf}$ but gets a suppression for $m\gtrsim m_{\rm inf}$ and the majority of inflation models predicts $m_{\rm inf}\gg H_{\rm inf}$. The final fermion abundance is roughly given by the sum of these two contributions. In any case, the fermion abundance goes to zero in the massless limit $m\to 0$ as expected, since a fermion is conformal in this limit. Combining them, the present fermion abundance in terms of its energy density divided by the entropy density $s$ is given by

[TABLE]

where

[TABLE]

For $H_{\rm inf}\ll m\lesssim m_{\rm inf}$, the fast contribution likely dominates and we have

[TABLE]

For $m\lesssim H_{\rm inf}$, on the other hand, the slow contribution likely dominates and we have

[TABLE]

In this case, the abundance is independent of the inflationary scale $H_{\rm inf}$.#5#5#5 If $m<H_{\rm R}$, the expression (33) should be multiplied by the factor $(m/H_{\rm R})^{1/2}$ (hence becomes independent of the reheating temperature) due to the change of expression of $k_{c}$, as noticed in the previous footnote.

Comparing it with the present DM abundance $\rho_{\rm DM}/s\sim 4\times 10^{-10}\,{\rm GeV}$, it is possible that a pure singlet fermion having only the gravitational interaction becomes a dominant component of DM if its mass is relatively large.

In Fig. 1 we illustrate with several contours of the fermion abundance on the plane of $(m,H_{\rm inf})$ with two choices of inflaton mass and reheating temperature, the left panel with $m_{\textrm{inf}}=10^{15}\,$GeV and $T_{\rm R}=10^{11}\,$GeV,　 and right panel with $m_{\textrm{inf}}=10^{13}\,$GeV and $T_{\rm R}=10^{11}\,$GeV. Three different contours (gray solid, blue dashed, and purple dotted) correspond to $\Omega_{\psi}=(1,0.1,0.01)\times\Omega_{\textrm{DM}}$, where $\Omega_{i}$ $(i=\psi,{\rm DM})$ is the density parameter defined by $\rho_{i}/\rho_{\rm crit}$ with $\rho_{\rm crit}$ being the critical energy density of the present universe. Evidently, wide parameter space exists for the correct DM relic density. Note that we should also include the contribution from thermal production by gravitational annihilation of SM particles in the thermal bath [18, 19, 20, 21]. See Appendix C for details. In the parameter space we have shown, however, contributions from thermal production is negligible.

3 Vector boson production

3.1 Vector boson action in the FRW Universe

Let us consider an action of massive vector boson,

[TABLE]

where $F_{\mu\nu}=\partial_{\mu}A_{\nu}-\partial_{\nu}A_{\mu}$. In the FRW background, this action is rewritten as

[TABLE]

One can impose a $Z_{2}$ symmetry under which only $A_{\mu}$ changes its sign to forbid the kinetic mixing with the standard model hypercharge photon. Then $A_{\mu}$ is stable and a candidate of DM. See Refs. [38, 39] for concrete model buildings.

The vector boson mass can be regarded as a result of the Higgs mechanism. In this case, the radial component of the Higgs boson is a physical field but it can be decoupled from the dynamics if the radial component is heavy enough. This is achieved by assuming that the gauge coupling constant is much smaller than the Higgs self coupling constant, for example. Or we can rely on the Stuckelberg mechanism: let the gauge boson mass term be $m^{2}g^{\mu\nu}(A_{\mu}+c\partial_{\mu}\phi)(A_{\nu}+c\partial_{\nu}\phi)$ by introducing additional real scalar field $\phi$. This mass term respects the gauge symmetry $A_{\mu}\to A_{\mu}+\partial_{\mu}\chi$ if $\phi$ transforms as $\phi\to\phi-\chi/c$ with arbitrary function $\chi$. By setting $\phi=0$ using this gauge degree of freedom, we end up with the massive vector boson action (34). In this case, there is no physical degree of freedom other than the vector boson.

It is again convenient to work with the Fourier mode since we are interested in the free vector boson:

[TABLE]

Since $A_{\mu}(\vec{x},t)$ is a real field, $A_{\mu}(\vec{k},t)=A_{\mu}^{*}(-\vec{k},t)$ must be satisfied. The three physical components are divided into the transverse and longitudinal mode: $\vec{A}=\vec{A}_{T}+\hat{k}A_{L}$, where the transverse mode satisfies $\vec{k}\cdot\vec{A}_{T}=0$ with the momentum vector $\vec{k}$ and $\hat{k}\equiv\vec{k}/|\vec{k}|$. Note that $A_{0}$ is not dynamical but it is determined by the constraint equation as

[TABLE]

Substituting this into the action, we find that the transverse mode and longitudinal mode are decoupled from each other: [40]

[TABLE]

The longitudinal mode is further redefined using the canonical field $\widetilde{A_{L}}\equiv f(\tau)A_{L}$ with $f(\tau)\equiv am/\sqrt{k^{2}+a^{2}m^{2}}$,

[TABLE]

where the effective mass is given by#6#6#6 The same expression was also derived in Ref. [41] assuming that the vector boson mass is generated by the Higgs mechanism.

[TABLE]

3.2 Transverse mode production

From the action of the transverse mode (39) it is evident that it is conformal in the limit $m\to 0$, i.e., the scale factor dependence disappears in the limit $m\to 0$. Hence there is no particle production in this limit. It is similar to the situation of a scalar field with conformal coupling.

To see this in detail, let us expand $\vec{A}_{T}$ as

[TABLE]

where $\vec{\epsilon}_{h}$ denotes the polarization vector for two polarization modes $h=+$ and $-$ which satisfies $\vec{\epsilon_{h}^{*}}\cdot\vec{\epsilon}_{h^{\prime}}=\delta_{hh^{\prime}}$. A concrete expression is $\vec{\epsilon}_{\pm}=(1,\pm i,0)/\sqrt{2}$ if $\vec{k}$ points to the $z$-direction. The ladder operators satisfy

[TABLE]

The mode function obeys the same equation of motion as a scalar field with conformal coupling:

[TABLE]

The solution during inflation (in the de Sitter phase) is given by

[TABLE]

where we have chosen the boundary condition at $k\tau\to-\infty$ (deep inside the horizon) so that it approaches to the mode function in the Minkowski space:

[TABLE]

Since $\nu^{2}<1/4$ for $m^{2}>0$, there is no growth of the superhorizon modes during inflation. Therefore, the dominant contribution to the gravitational production of the transverse vector boson happens around/after the end of inflation. We can solve the equation of motion (45) with the initial condition (47) under the background evolution of the cosmic scale factor $a(\tau)$ due to the (spatially homogeneous) inflaton dynamics to determine the energy density of the transverse vector boson through (84). Instead of solving (45) directly, we can make useful parametrization as follows:

[TABLE]

where $\omega_{k}\equiv\sqrt{k^{2}+a^{2}m^{2}}$ and $\alpha_{k}(\tau)$ and $\beta_{k}(\tau)$ are assumed to satisfy

[TABLE]

It is checked that these set of equations satisfy the equation of motion (45). The initial condition is taken to be $\alpha_{k}\to 1$ and $\beta_{k}\to 0$ at $k\tau\to-\infty$. Using this parametrization, the energy density (84) is expressed as

[TABLE]

Here $f_{T}(k,\tau)$ denotes the occupation number of the transverse vector boson. Note that the normalization condition implies $|\alpha_{k}(\tau)|^{2}-|\beta_{k}(\tau)|^{2}=1$. In contrast to the case of fermion, this normalization condition does not limit the possible produced number density of vector boson.

Then what we have to do is to estimate $\beta_{T}(k,\tau)$. The calculation is the same as the gravitational production of a scalar field with conformal coupling as performed in detail in Ref. [14]. Here we present only the results. The number density of the transverse vector boson is given by

[TABLE]

where $\eta$ is given by the same expression as (31) after reinterpreting $m$ in (31) as the vector boson mass. It is assumed that $m\ll m_{\rm inf}$ since otherwise the vector boson production is suppressed. Taking account of the two polarization degrees of freedom, the numerical coefficient $\mathcal{C}_{T}$ is found to be $3/(256\pi)$ if the inflaton potential is well approximated by the quadratic one [15]. Assuming that the universe is matter-dominated before the completion of reheating, we obtain the energy to entropy density ratio as

[TABLE]

3.3 Longitudinal mode production

The longitudinal mode is more similar to a scalar field, as seen from the action (41). It is quantized as

[TABLE]

where the ladder operators satisfy

[TABLE]

The equation of motion of the mode function is

[TABLE]

For convenience, we also present the equation of motion in the original basis:

[TABLE]

During the de Sitter phase, the effective mass of the longitudinal mode is given by

[TABLE]

In the high momentum limit $k\gg am$, it is approximated as

[TABLE]

which is the same form as a massive scalar coupled to gravity minimally. Thus we can take the mode function as in the Minkowski form in the high momentum limit $(k/a\gg{\rm max}[m,H_{\rm inf}])$ as

[TABLE]

In the low momentum limit $k\ll am$, we obtain

[TABLE]

This is always positive definite even in the massless limit $m\to 0$, which is a unique feature of massive vector boson, different from a massive scalar. Now we consider two cases separately: heavy vector boson $m\gtrsim H_{\rm inf}$ and light vector boson $m\lesssim H_{\rm inf}$.

3.3.1 Heavy vector boson case

For the heavy vector boson case $m\gtrsim H_{\rm inf}$, it is evident that the effective mass squared $m_{L}^{2}$ is always positive independently of the wavenumber $k$, and hence there is no significant growth of the vacuum fluctuation. In particular, no superhorizon modes are enhanced during the de Sitter phase. Therefore, in this case, we should only take account of the production of high momentum modes after inflation. For $k\gg am$, we obtain

[TABLE]

This is the same form as the minimally coupled scalar field. This may be regarded as a consequence of the Goldstone boson equivalence theorem, which says that the longitudinal vector boson may be identified with the Goldstone boson in the high energy limit. Similar to the case of scalar field, we can again make the following parameterization:

[TABLE]

where $\alpha_{k}(\tau)$ and $\beta_{k}(\tau)$ are assumed to satisfy

[TABLE]

It is again checked that these set of equations satisfy the equation of motion (55). The initial condition is taken to be $\alpha_{k}\to 1$ and $\beta_{k}\to 0$ at $k\tau\to-\infty$. The energy density (85) is expressed as

[TABLE]

where $f_{L}(k,\tau)$ denotes the occupation number of the longitudinal vector boson. The normalization condition implies $|\alpha_{k}(\tau)|^{2}-|\beta_{k}(\tau)|^{2}=1$. In this case, therefore, the number density of produced longitudinal vector boson is estimated in the same as a minimal scalar field [14],

[TABLE]

where $\mathcal{C}_{L}\sim 3/(512\pi)$ [15] and we assumed $m<m_{\rm inf}$. Compared with the transverse mode, there is no suppression factor of $(m/m_{\rm inf})^{4}$. Thus the energy to entropy density ratio is evaluated as

[TABLE]

3.3.2 Light vector boson case

For the light vector boson case $m\lesssim H_{\rm inf}$, the situation is different. Superhorizon modes with $m<k/a<H_{\rm inf}$ experience tachyonic growth during inflation, similar to the case of light scalar and these inflation-generated long wavelength modes may give dominant contributions to the final vector boson abundance. Note that the growth is eventually terminated when the physical wavenumber $k/a$ becomes equal to the vector boson mass $m$. Thus the power spectrum of such a massive vector boson at large scale is much suppressed compared with the case of light scalar field with the same mass [40]. For modes with $m<k/a<H_{\rm inf}$, it is easily found from (55) that the mode function grows as#7#7#7 Note that the original field $A_{L}(k,\tau)$ before the canonical rescaling remains constant.

[TABLE]

Therefore, at the end of inflation, the mode function becomes

[TABLE]

In terms of the energy density per log frequency at the end of inflation, after the renormalization of the UV divergence as usual, we have

[TABLE]

At this stage, therefore, the shortest wavelength mode $(k\sim H_{\rm inf}a_{\rm end})$ gives the dominant contribution to the total energy density of the longitudinal vector boson. However, these modes are (highly) relativistic and receive larger suppression factor due to the redshift. As shown in Appendix, the energy density scales as

[TABLE]

Assuming again that the universe is matter-dominated before the completion of the reheating, the final energy to entropy density ratio is expressed in terms of the number density at $H=m$, defined by $n_{L}(k;H=m)=\rho_{L}(k;H=m)/\sqrt{m^{2}+k^{2}/a^{2}(H=m)}$:

[TABLE]

where $H_{\rm R}\equiv\sqrt{\pi^{2}g_{*}/90}\,T_{\rm R}^{2}/M_{P}$ is the Hubble parameter at the completion of reheating.

For $H_{\rm R}<m$, we find

[TABLE]

where $k_{*}\equiv a(H=m)m$ denotes the comoving wavenumber for which the corresponding mode becomes non-relativistic at $H=m$. Thus it is seen that the total energy density is dominated by modes around $k=k_{*}$. In this case we have $k_{*}=a_{\rm end}m(H_{\rm inf}/m)^{2/3}$. The total vector boson abundance is then given by

[TABLE]

It is independent of the vector boson mass $m$.

For $H_{\rm R}>m$, after straightforward but tedious calculations, we find

[TABLE]

Again we find that it is peaked around $k=k_{*}$ where $k_{*}\equiv a(H=m)m$, which is now evaluated as $k_{*}=a_{\rm end}m(H_{\rm inf}/H_{\rm R})^{2/3}(H_{\rm R}/m)^{1/2}$. Thus the total vector boson abundance is given by

[TABLE]

It is consistent with Ref. [40].

3.4 Combined Results

Let us summarize the results so far. The abundance of longitudinal vector boson which is purely gravitationally produced is given by

[TABLE]

Note that the origin of the vector boson is different between the case of $m>H_{\rm inf}$ and $m<H_{\rm inf}$. In the former case, the vector boson is dominantly produced at the end of inflation or during the early stage of reheating and the main produced mode is about the inflaton mass: $k\sim a_{\rm end}m_{\rm inf}$. In the latter case, the dominant contribution comes from the superhorizon mode generated during inflation, which eventually re-enters the horizon at $H\sim m$. The transverse modes are also produced at the end of inflation and during the reheating stage, but they are always subdominant compared with the longitudinal mode.

In Fig. 2, we show several contours of the vector boson abundance on the parameter space of $(m,H_{\rm inf})$ for two sets of inflaton mass and reheating temperature, $m_{\textrm{inf}}=10^{13}\,$GeV and $T_{\rm R}=10^{11}\,$GeV (Left panel), $m_{\textrm{inf}}=10^{12}\,$GeV and $T_{\rm R}=10^{10}\,$GeV (Right panel). Similar to the fermion case, thermal production is included, see Appendix C. Three contours (gray solid, blue dashed, and purple dotted) correspond to $\Omega_{A}=(1,0.1,0.01)\times\Omega_{\textrm{DM}}$ where $\Omega_{A}=\rho_{A}/\rho_{\rm crit}$ is the density parameter of the vector boson. We can see there are wide and viable parameter regions that can satisfy the current DM relic abundance. Contours without including thermal productions are shown in thin curves that however almost coincide with thick ones, which means thermal contributions are negligible in the showed parameter space.

4 Conclusions and discussion

We have studied the DM production mechanism in the case where the DM particle is a massive fermion or vector boson and has only the gravitational interaction. The production takes place through the so-called gravitational particle production under the standard inflationary cosmology.

For the case of a massive fermion, the presence of mass term violates the conformal invariance and it somehow feels the background time evolution, resulting in particle production. The dominant production process depends on the fermion’s mass $m$. For $m\lesssim H_{\rm inf}$, the non-adiabaticity of the fermion wave function is prominent when the fermion becomes non-relativistic $k\sim am$ for each wavenumber $k$. Those with momentum $k$ such that $k\sim am$ and $H\sim m$ gives the dominant contribution to the final fermion abundance as already pointed out in Ref. [30]. For $m\gg H_{\rm inf}$, such an effect of the universe expansion is negligible while the inflaton coherent oscillation produces excites the high momentum fermion modes, since the cosmic scale factor $a(\tau)$ includes a small but nonzero oscillating part. In both cases, we have the viable parameter regions that can reproduce the present DM abundance. All these features are similar to the case of a scalar field with conformally coupled to gravity [14].

For the case of a massive vector boson, the story is a bit complicated. The transverse mode is conformal in the massless limit, and hence the gravitational production proceeds only through the presence of mass term. Again it is similar to the case of conformally coupled scalar field. On the other hand, the longitudinal mode shows more non-trivial behavior. For $m\lesssim H_{\rm inf}$, during the de Sitter phase the vector obtains superhorizon quantum fluctuations and eventually behaves as non-relativistic matter. In contrast to the scalar field with minimal coupling, there is a limit for the growth of the superhorizon modes at $k\sim am$, and hence such a model is not constrained by the presence of DM isocurvature perturbation on cosmological scales [40]. For $m\gg H_{\rm inf}$, it is rather close to the minimally coupled scalar field, and the inflaton coherent oscillation excites the high-momentum longitudinal mode. In both cases we have correct parameter regions that reproduce the present DM abundance.

Several comments are in order. For analyses of both the fermion and vector boson, we have assumed that the fermion/vector boson mass parameter $m$ is just a constant. It is possible that the mass is given by the expectation of value of some other Higgs-like scalar field. In such a case the mass parameter $m$ can be dynamical. For example, it can have a different value in the inflationary era and the present universe. Our results crucially depend on this assumption. If there is a Higgs field, we should take account of the dynamics of the Higgs field if the mass of the Higgs field is smaller than the Hubble parameter $H$, which may significantly affect our estimate of the fermion and vector boson abundance.

Finally we comment on the detectability of our model. Since DM is only gravitationally interacting, it is hopeless to find it through a kind of direct detection experiments. If DM is long-lived but has finite lifetime, its decay would produce high-energy cosmic rays for indirect detections [42]. As seen from Figs. 1 and 2, the wide DM mass range is consistent with the DM abundance and correspondingly the energy of cosmic-rays induced by the DM decay also can take wide range of values. Another interesting possibility may be to search for the effect of heavy field (say, $m\sim H_{\rm inf}$) through the non-Gaussianity of the primordial curvature fluctuations [43, 44, 45, 46, 47, 17].

Acknowledgments

The work of Y.E. was supported in part by JSPS Research Fellowships for Young Scientists. This work was supported by the Grant-in-Aid for Scientific Research C (No.18K03609 [KN]) and Innovative Areas (No.26104009 [KN], No.15H05888 [KN], No.17H06359 [KN], 16H06490 [YT]).

Appendix A Conventions

Our convention follows Ref. [48]. The flat space metric is taken to be $\eta_{\mu\nu}={\rm diag}(-1,1,1,1)$. The gamma matrices in the flat space is denoted by $\gamma_{a}$ while those in the curved background are expressed by $\hat{\gamma}_{\mu}$. They are related by $\hat{\gamma}_{\mu}=e^{a}_{\mu}\gamma_{a}$ by using the vierbein $e^{a}_{\mu}$. In the FRW background, the vierbein is

[TABLE]

so that the metric is given by $g_{\mu\nu}=e^{a}_{\mu}e^{b}_{\nu}\eta_{ab}=a^{2}\eta_{\mu\nu}$. The Clifford algebra is then defined as

[TABLE]

Thus we have $(\hat{\gamma}^{0})^{\dagger}=-\hat{\gamma}^{0}$ and $(\hat{\gamma}^{i})^{\dagger}=\hat{\gamma}^{i}$. We also define $\gamma_{5}\equiv i\gamma_{0}\gamma_{1}\gamma_{2}\gamma_{3}$. The explicit Weyl representation of the gamma matrices is

[TABLE]

where $\sigma_{\mu}=(-1,\sigma_{i})$ and $\bar{\sigma}_{\mu}=(1,\sigma_{i})$ with $\sigma_{i}$ being the usual Pauli matrices that satisfy $\{\sigma_{i},\sigma_{j}\}=2\delta_{ij}$. In the Dirac representation the explicit form of the gamma matrices is

[TABLE]

Appendix B Energy density of vector boson

B.1 Energy-momentum tensor and energy density

The energy momentum tensor is defined as

[TABLE]

For a massive vector boson, it is given by

[TABLE]

After some calculations, we find that the energy density $\rho=\left<0|T_{00}|0\right>$ is given by

[TABLE]

where the prefactor $2$ in the expression of $\rho_{T}$ accounts for the two polarization states. Again we can rewrite the longitudinal mode by using the canonical field $\widetilde{A}_{L}(\vec{k})=fA_{L}(\vec{k})$ with $f(\tau)\equiv am/\sqrt{k^{2}+a^{2}m^{2}}$,

[TABLE]

Note that these expressions show UV divergence when the vacuum configuration for the mode function is substituted, which must be subtracted to obtain nearly vanishing cosmological constant observed now.

B.2 Energy scaling of the vector boson

In this appendix we present approximate solutions to the equation of motion of the vector boson mode function after inflation and show how the corresponding energy density scales with respect to the cosmic scale factor $a(\tau)$. Hereafter we assume the reheating phase of the universe with equation of state $w$ such that $H\propto a^{-3(1+w)/2}$.

B.2.1 Transverse mode

The equation of motion for the transverse mode is

[TABLE]

In the light vector boson regime $(m\ll H)$, the solution looks like#8#8#8 Actually there is also a solution $\mathcal{A}_{T}\propto a^{(1+3w)/2}$ for $k\ll am$ and this is a growing solution for $w>-1/3$. It is not trivial whether this growing solution is of physical importance or not, but detailed investigations show that the boundary condition (47) does not lead to this growing solution [40].

[TABLE]

Thus $\rho_{T}\propto a^{-4}$ for $k\gg am$ and $\rho_{T}\propto a^{-2}$ for $k\ll am$.

In the heavy vector boson regime $(m\gg H)$, the solution looks like

[TABLE]

Thus $\rho_{T}\propto a^{-4}$ for $k\gg am$ and $\rho_{T}\propto a^{-3}$ for $k\ll am$, as just expected from the scaling of relativistic and non-relativistic matter, respectively.

B.2.2 Longitudinal mode

The equation of motion for the longitudinal mode is

[TABLE]

where the effective mass is given by#9#9#9 Some relations in the universe with equation of state $w$ are: $a^{\prime\prime}/a=\mathcal{H}^{2}(1-3w)/2$, $\mathcal{H}\tau=2/(1+3w)$, $\mathcal{H}^{\prime}=-\mathcal{H}^{2}(1+3w)/2$, $a\propto\tau^{2/(1+3w)}$

[TABLE]

In the case of light vector boson regime $(m\ll H)$, it is a bit complicated:

[TABLE]

Thus the solution looks like

[TABLE]

with $c_{1}$–$c_{4}$ denoting numerical constants determined by the initial conditions at the end of inflation.#10#10#10 These solutions can be found by substituting the ansatz $\widetilde{\mathcal{A}_{L}}\propto a^{\ell}$ and requiring that leading terms vanish to determine the exponent $\ell$. Here we neglect the last terms in (92) since they do not much affect the value of $\ell$.

Note again that the growing solution for $k\ll am$ actually is not of physical importance if one properly connects the solution during inflation to that during the reheating phase [40]. Thus we only need to see terms with $c_{1}$ and $c_{3}$ for $k\ll aH$. As a result, we obtain the scaling of the longitudinal vector boson energy density as $\rho_{L}\propto a^{-4}$ for $k\gg aH$ (subhorizon modes) and $\rho_{L}\propto a^{-2}$ for $k\ll aH$ (superhorizon modes).

In the case of heavy vector boson regime $(m\gg H)$, on the other hand, we simply obtain

[TABLE]

It is the same as the transverse one and hence we obtain $\rho_{L}\propto a^{-4}$ for $k\gg am$ and $\rho_{L}\propto a^{-3}$ for $k\ll am$, just as expected.

Appendix C Thermal production

Here we briefly summarize the thermal production (TP) of gravitationally interacting DM. The production cross section of purely gravitational DM $X$, either $X$ is fermion or vector boson, through the scattering of SM particles in the thermal bath with the temperature $T$ is [18, 19, 20, 21]

[TABLE]

with some numerical constant $y,y\simeq 0.2$ for fermions and $y\simeq 0.8$ for vectors. See Ref. [20] for details. The DM number density created per Hubble time is given by

[TABLE]

assuming $T>m_{X}$. Otherwise, the production rate is suppressed by $e^{-m_{X}/T}$. Therefore, if $T_{\rm R}>m_{X}$, the production is dominated at $T\sim T_{\rm R}$. On the other hand, if $T_{\rm R}<m_{X}$, the production is dominated at $T\sim m_{X}$. Thus the DM abundance through the gravitational annihilation of SM particles is given by

[TABLE]

Here the subscript “prod” refers to the dominant production epoch. It is estimated as

[TABLE]

where $\tilde{y}$ includes various $\mathcal{O}(1)$ numerical factor neglected in the rough estimation. This should be compared with the gravitational production studied in the main text. For the case of vector boson or minimal scalar, it can dominate over the gravitational production if the reheating is close to instantaneous, i.e., if $H_{\rm R}$ is not very far from $H_{\rm inf}$. For the fermion case, it depends on the model parameters $m_{X},H_{\rm inf},T_{\rm R}$ and $m_{\rm inf}$ in a more complicated way.