# Gravity as a Portal to Reheating, Leptogenesis and Dark Matter

Basabendu Barman  Simon Cléry  Raymond T. Co  Yann Mambrini  Keith A. Olive

###### Abstract

We show that a minimal scenario, utilizing only the graviton as an
intermediate messenger between the inflaton, the dark sector and the Standard
Model (SM), is able to generate simultaneously the observed relic density of
dark matter (DM), the baryon asymmetry through leptogenesis, as well as a
sufficiently hot thermal bath after inflation. We assume an inflaton potential
of the form $V(\phi)\propto\phi^{k}$ about the minimum at the end of
inflation. The possibility of reheating via minimal gravitational interactions
has been excluded by constraints on dark radiation for excessive gravitational
waves produced from inflation. We thus extend the minimal model in several
ways: i) we consider non-minimal gravitational couplings–this points to the
parameter range of DM masses $M_{N_{1}}\simeq 2-10$ PeV, and right-handed
neutrino masses $M_{N_{2}}\simeq(5-20)\times 10^{11}$ GeV, and
$T_{\text{RH}}\lesssim 3\times 10^{5}$ GeV (for $k\leq 20$); ii) we propose an
explanation for the PeV excess observed by IceCube when the DM has a direct
but small Yukawa coupling to the SM; and iii) we also propose a novel
scenario, where the gravitational production of DM is a two-step process,
first through the production of two scalars, which then decay to fermionic DM
final states. In this case, the absence of a helicity suppression enhances the
production of DM and baryon asymmetry, and allows a great range for the
parameters including a dark matter mass below an MeV where dark matter
warmness can be observable by cosmic 21-cm lines, even when gravitational
interactions are responsible for reheating. We also show that detectable
primordial gravitational wave signals provide the opportunity to probe this
scenario for $T_{\text{RH}}\lesssim 5\times 10^{6}$ GeV in future experiments,
such as BBO, DECIGO, CE and ET.

PI/UAN-2022-720FT, UMN–TH–4202/22, FTPI–MINN–22/26

##  1 Introduction

Since the first calculation of the mass of the Milky Way by Henry Poincaré in
1906 [1], and the conclusion that the “dark matter is present in much greater
amount than luminous matter” by Fritz Zwicky in 1933 [2], the virial method
has been used frequently to compute the amount of dark matter (DM) in the
Universe. The presence of dark matter is generally deduced from its
gravitational effects. The precise abundance of DM is obtained from
observations of the anisotropy spectrum of the cosmological microwave
background (CMB) [3]. Initially, a neutrino, or more generally a heavy neutral
lepton was thought to be an ideal dark matter candidate [4]. This candidate
was assumed to interact weakly with the Standard Model (SM) and required a GeV
scale mass to satisfy relic abundance constraints [5, 6]. Generalizations of
the heavy neutrino DM candidate are referred to as WIMPs (weakly interacting
massive particles). Perhaps the best studied WIMP is the lightest
supersymmetric particle in supersymmetric extensions of the SM [7]. Other well
studied candidates include those generated by a Higgs-portal [8, 9, 10, 11,
12, 13, 14, 15, 16, 17, 18] or a $Z$-portal [19, 20]. However, these minimal
constructions are now heavily constrained (see, for example, Refs. [21, 22]
for reviews), and even extensions to $Z^{\prime}$-portal [23, 24, 25, 26, 27]
are in tension with the electroweak nature of dark matter.

The WIMP relic density is often determined by thermal freeze-out. WIMPs are
assumed to be in thermal equilibrium at temperatures higher than the mass, and
the relic density is determined by the equilibrium density when DM
annihilations can no longer keep up with the expansion rate of the Universe.
It is also possible that particles with interactions much weaker than
electroweak interactions and are never fully in equilibrium, but are
nonetheless produced in the thermal bath after inflation. An example of such a
candidate is the gravitino [28, 7, 29, 30]. More generally, these Feebly
Interacting Massive Particles (FIMPs) have been proposed [31, 32, 33] as an
alternative to WIMPs (see [34] for a recent review). In this framework, the DM
candidate is never thermalized due its extremely weak coupling to the SM, so
weak that they evade the current accelerator constraints.

In the early Universe, FIMPs can be produced from either the decay or
annihilation of states in the visible sector. When the SM temperature becomes
smaller than the typical mass scale of the interaction (i.e. the maximum of
the DM and the mediator mass), the generation process becomes suppressed,
leaving a constant comoving DM number density. Such a scenario is often
referred as the freeze-in mechanism [32]. In contrast to the “WIMP-miracle”
which produces the observed relic density with near weak-scale couplings and
masses, a “FIMP-miracle” occurs when one considers renormalizable couplings of
order $\sim\mathcal{O}(10^{-11})$ independent of the mass of the DM. If a
priori such couplings seem unatural, UV versions of the freeze-in mechanism
may invoke effective couplings, suppressed by a large mass scale above the
temperature of the thermal bath. This can be achieved via non-renormalizable
operators [35], suppressed by a high mass scale, e.g., in models where the
mediators between the visible sector and the dark sectors are very massive.
This is the case in unified theories like $SO(10)$ with a heavy $Z^{\prime}$
gauge boson [33, 36, 37], moduli fields [38], high scale SUSY [39, 40, 41, 42,
43] or heavy spin-2 constructions [44]. In other examples, freeze-in of DM may
proceed via loops [45, 46] or 4-body final states [47]. All of these scenarios
are particularly interesting, as the DM yield is sensitive to the highest
temperature $T_{\text{max}}$ reached by the SM plasma [48, 49, 50, 51, 52,
53], controlled by the dynamics of the inflaton decay.

Even feebler interactions are possible when the only effective coupling at the
UV scale is gravity. Indeed, the minimal irreducible interaction that should
exist between DM and the Standard Model (SM) is mediated by graviton exchange
[54, 55, 56, 57, 58, 59, 44, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71]
which can lead to the observed amount of DM through the scattering of the
particles in the thermal bath or directly through the gravitational transfer
of the energy stored in the inflaton condensate, as already been discussed in
detail Refs. [65, 66, 67, 68, 69].

DM requires an extension to the SM, but it is not the only reason why an
extension is necessary. As is well known, the visible or baryonic matter
content of the Universe is asymmetric. One interesting mechanism to produce
the baryon asymmetry of the Universe (BAU) via the lepton sector physics is
known as leptogenesis [72], where, instead of creating a baryon asymmetry
directly, a lepton asymmetry is generated first and subsequently gets
converted into baryon asymmetry by the $(B+L)$-violating electroweak sphaleron
transitions [73]. In thermal leptogenesis [74, 75, 76, 77], the decaying
particles, typically right-handed neutrinos (RHNs), are produced thermally
from the SM bath. However, the lower bound on the RHN mass in such scenarios
(known as the Davidson-Ibarra bound), leads to a lower bound on the reheating
temperature $T_{\text{RH}}\gtrsim 10^{10}$ GeV [78] so that the RHNs can be
produced from the thermal bath. One simpler alternative is the non-thermal
production of RHNs [79, 80, 81, 82, 83] originating from the decay of
inflaton. This interaction is necessarily model dependent as it depends on the
Yukawa interaction between the inflaton and the RHNs.

In addition to providing the DM abundance, gravitational interactions can also
be the source of baryogenesis. As shown in [84], it is possible to have a
model-independent theory of non-thermal production of RHNs from inflation,
once the inflaton potential is specified.111The simultaneous generation of
gravitational DM and the baryon asymmetry was also discussed in [85]. Our
results differ, as their choices of parameters are in conflict with the
tensor-to-scalar ratio bound from Planck. The abundance of RHNs is calculated
in the same manner as the dark matter abundance and can lead to observed BAU
from the out-of-equilibrium CP violating decay of the RHNs, produced during
the reheating epoch.

As noted above, the DM and RHNs may be produced through gravitational
interactions emanating from the thermal bath or directly from the inflaton
condensate. It has also been argued that the thermal bath itself may be
generated from gravitational interactions [86, 68, 84]. However, reheating the
Universe from graviton exchange processes alone requires a steep inflaton
potential during reheating, resulting in a low reheating temperature and a
massive enhancement of tensor modes after inflation. Hence, the minimal
scenario of gravitational reheating is excluded by an excessive generation of
dark radiation in the form of gravitational waves (GWs) during BBN, as already
noted in [87]. This limitation of minimal gravitational reheating is one
motivation to introduce, as a natural generalization, non-minimal couplings of
fields with gravity.

Motivated by these arguments, we derive a simultaneous solution for the DM
abundance, the baryon asymmetry, and the origin of the thermal bath from
purely gravitational interactions. In this sense, our scenario can be
considered as the most minimal possible, since we do not introduce any new
interactions for any process beyond the SM, except for gravity. The only new
fields required are the dark matter candidate and the RHNs (which are anyway
needed for the generation of neutrino masses). Our only model dependence comes
from the choice of the particular inflaton potential. However we are mostly
sensitive to the shape of the potential about the minimum after inflation. To
be definite, we adopt the class of inflationary models called T-models [88].
But, as will be shown, even this dependence proves to be weak when it comes to
combining the constraints of reheating, baryogenesis, and the dark matter
relic density. We further show that the present framework can give rise to a
detectable inflationary GW background, that in turn excludes the minimal
gravitational reheating scenario which leads to an excess of the present-day
GW energy density, in conflict with the BBN prediction. However, a large part
of the parameter space still remain within the reach of several futuristic GW
detection facilities.

The paper is organized as follows. After presenting our framework in Sec. 2,
we review the gravitational production from inflaton scattering and the
thermal bath in Sec. 3, where we also discuss the effect of non-minimal
gravitational interactions. We derive the set of parameters (dark matter mass,
RHN mass, and reheating temperature, $T_{\text{RH}}$, which simultaneously
provide the correct relic density and BAU in Sec. 4. If the dark matter is not
absolutely stable, we are able to propose an explanation for the PeV events
observed at IceCube in the case of a long-lived candidate. Finally, we propose
a novel scenario where the gravitational production is a two-step process
passing through a scalar singlet which couples with the RHN sector in Sec. 5,
before concluding in Sec. 6.

##  2 The framework

If the metric is expanded around Minkowski space-time:
$g_{\mu\nu}\simeq\eta_{\mu\nu}+\frac{2h_{\mu\nu}}{M_{P}}$, then the
gravitational interactions are described by the Lagrangian [89, 90]

$\displaystyle\sqrt{-g}\,\mathcal{L}_{\rm
int}=-\frac{1}{M_{P}}\,h_{\mu\nu}\,\left(T^{\mu\nu}_{\rm
SM}+T^{\mu\nu}_{\phi}+T^{\mu\nu}_{X}\right)\,,$ (2.1)

where $\phi$ is the inflaton and $X$ is a particle which does not belong to
the SM.222$M_{P}=(8\pi G_{N})^{-1/2}\simeq 2.4\times 10^{18}$ GeV is the
reduced Planck mass. In the present context we consider $X$ to be a spin 1/2
Majorana fermion which can be associated with the dark matter or a right-
handed neutrino. The graviton propagator for momentum $p$ is

$\Pi^{\mu\nu\rho\sigma}(p)=\frac{\eta^{\rho\nu}\eta^{\sigma\mu}+\eta^{\rho\mu}\eta^{\sigma\nu}-\eta^{\rho\sigma}\eta^{\mu\nu}}{2p^{2}}\,.$
(2.2)

The form of the stress-energy tensor $T^{\mu\nu}_{i}$ depends on the spin of
the field and, for Majorana spin-1/2 fermions, takes the form

$T^{\mu\nu}_{1/2}=\frac{i}{8}\left[\bar{\chi}\gamma^{\mu}\overset{\leftrightarrow}{\partial^{\nu}}\chi+\bar{\chi}\gamma^{\nu}\overset{\leftrightarrow}{\partial^{\mu}}\chi\right]-g^{\mu\nu}\left[\frac{i}{4}\bar{\chi}\gamma^{\alpha}\overset{\leftrightarrow}{\partial_{\alpha}}\chi-\frac{m_{\chi}}{2}\,\overline{\chi^{c}}\chi\right]\,,$
(2.3)

whereas for a scalar $\varphi$,

$T^{\mu\nu}_{0}=\partial^{\mu}\varphi\partial^{\nu}\varphi-g^{\mu\nu}\left[\frac{1}{2}\partial^{\alpha}\varphi\,\partial_{\alpha}\varphi-V(\varphi)\right]\,.$
(2.4)

There are of course many possible scalar potentials $V(\phi)$ which can
account for inflation. However, the calculations relevant in this paper are
largely independent of the potential during inflation and depend only on the
shape of the potential about the minimum. Without loss of generality, we will
assume that $V(\phi)$ is among the class of $\alpha$-attractor T-models [88]

$V(\phi)\;=\;\lambda
M_{P}^{4}\left|\sqrt{6}\tanh\left(\frac{\phi}{\sqrt{6}M_{P}}\right)\right|^{k}\,,$
(2.5)

which can be expanded about the origin333Our discussion is general and not
limited to T-models of inflation as the way we express the minimum of the
potential is generic.

$V(\phi)=\lambda\frac{\phi^{k}}{M_{P}^{k-4}}\,;\quad\phi\ll M_{P}\,.$ (2.6)

In this class of models, inflation occurs at large field values
($\phi>M_{P}$), and after the period of exponential expansion, the inflaton
begins to oscillate about the minimum and the process of reheating begins.
]The end of inflation may be defined when $\ddot{a}=0$ where $a$ is the
cosmological scale factor. The inflaton field value at that time is given by
[91, 52] as

$\displaystyle\phi_{\rm end}$
$\displaystyle\simeq\sqrt{\frac{3}{8}}M_{P}\ln\left[\frac{1}{2}+\frac{k}{3}\left(k+\sqrt{k^{2}+3}\right)\right].$
(2.7)

It is easy to show that at the end of inflation, the condition $\ddot{a}=0$ is
equivalent to $\dot{\phi}_{\rm end}^{2}=V(\phi_{\rm end})$ and thus the
inflaton energy density at $\phi_{\rm end}$ is $\rho_{\rm
end}=\frac{3}{2}V(\phi_{\rm end})$. The overall scale of the potential
parameterized by the coupling $\lambda$, can be determined from the amplitude
of the CMB power spectrum $A_{S}$,

$\displaystyle\lambda$
$\displaystyle\simeq\frac{18\pi^{2}A_{S}}{6^{k/2}N^{2}},$ (2.8)

where $N$ is the number of e-folds measured from the end of inflation to the
time when the pivot scale $k_{*}=0.05~{}{\rm Mpc}^{-1}$ exits the horizon. In
our analysis, we use $\ln(10^{10}A_{S})=3.044$ [92] and set $N=55$. This leads
to an inflaton mass of $m_{\phi}\simeq 1.2\times 10^{13}$ GeV for $k=2$. More
generically, $m_{\phi}\simeq 1.2\times 10^{13}$ is also the inflaton mass at
the end of inflation for any larger $k$ when the full potential in Eq. (2.5)
is used. While $N=55$ is appropriate for reheating temperatures of order
$10^{12}$ GeV, for lower reheating temperatures (between $10-10^{7}$ GeV)
$N=45-50$ [93]. However, we have checked that our results are very insensitive
to the value of $N$.

In addition to the inflationary sector and the SM, neutrino masses and mixing
require at least two (heavy) right-handed neutrino states for the seesaw
mechanism [94, 95, 96, 97, 98, 99]. One of these, if produced and remaining
out-of-equilibrium until its decay, can produce a lepton asymmetry. In order
to account for the dark matter in a most economic way, we assume three RHNs,
$N_{i}$, where for now, the lightest of these, $N_{1}$ is decoupled from the
other two and has a vanishing Yukawa coupling. Aside from the Yukawa
couplings, the only couplings we consider between the SM, the RHNs, and the
inflaton are gravitational of the form in Eq. (2.1). Needless to say, such
interactions are unavoidable, and must be taken into account in any extensions
beyond the SM.

As a concrete example, we consider the renormalizable interaction Lagrangian
between the Majorana RHNs and the SM

$\displaystyle\mathcal{L}\supset-\frac{1}{2}\,M_{N_{i}}\,\overline{N_{i}^{c}}N_{i}-(y_{N})_{ij}\,\overline{N}_{i}\,\widetilde{H}^{\dagger}\,L_{j}+{\rm
h.c.}\,.$ (2.9)

Here $H$ and $L$ are the SM Higgs and lepton doublet respectively. Lepton
number is clearly violated in this Lagrangian.444We consider the RHNs to be
mass diagonal. For now, we assume that $(y_{N})_{1i}=0$ for all $i$ and that
$N_{1}$ is stable. As a result, $N_{1}$ is a viable DM candidate. Later, we
will relax this condition and consider a metastable DM candidate with
$(y_{N})_{1i}\neq 0$, allowing for $N_{1}$ to decay into neutrinos that could
be observed at IceCube. The preservation of the lepton asymmetry will provide
a limit on $(y_{N})_{1i}/M_{N_{1}}$. The other two RHNs, namely $N_{2,3}$ are
assumed to be heavier and they participate in leptogenesis.

We would like to remind the readers that there are three types of seesaw
models, which differ by the properties of the exchanged heavy particles, e.g.,

  * (i)

Type-I: SM gauge fermion singlets

  * (ii)

Type-II: SM $SU(2)_{L}$ scalar triplets

  * (iii)

Type-III: SM $SU(2)_{L}$ fermion triplets.

In the present case we are considering the Type-I scenario, which is evident
from the Lagrangian in Eq. (2.9). The Type-I seesaw mechanism can indeed be
realized with only two active right-handed neutrino [100, 101, 102, 103]. In
this context, only the normal ($m_{1}<m_{2}<m_{3}$) and inverted
($m_{3}<m_{2}<m_{1}$) hierarchies are relevant, where $m_{i}$ are the light
neutrino masses. With only two RHN playing a role in the seesaw mechanism, we
expect $m_{1}=0$ or $m_{3}=0$, depending on the hierarchy. Indeed, due to the
reduced rank of the mass matrices (a $3\times 2$ Dirac matrix and a $2\times
2$ Majorana matrix) one neutrino remains massless, while the others acquire
their light mass through the usual seesaw suppression of the order
$m_{i}\sim\frac{(y_{N})_{ii}^{2}\langle H\rangle^{2}}{M_{N_{i}}}$.

We further assume the absence of any direct coupling between the inflaton
$\phi$ and the RHNs, such that there is no perturbative decay of the inflaton
into the RHN final state; in other words we do not attribute a lepton number
to the inflaton. Thus, the only possible production of the RHNs is the 2-to-2
gravitational scattering of the inflatons and of the particles in the
radiation bath. As we will show, these production channels dominate in
different regions of the parameter space. In Fig. 1, we show the $s$-channel
exchange of a graviton obtained from the Lagrangian in Eq. (2.1) for the
production of RHNs from the inflaton condensate, to which we can add a similar
diagram for the production of SM fields during the reheating process. Despite
the a priori Planck reduced interactions, we will show that this framework is
perfectly capable of simultaneously explaining the dark matter relic abundance
and the observed baryon asymmetry, while also reheating the Universe. The
Planck suppression due to graviton exchange is indeed partially compensated by
the energy available in the inflaton condensate at the end of inflation.

Figure 1: Feynman diagram for the production of RHN through the gravitational
scattering of the inflaton condensate. A similar diagram also exists with
Standard Model particles in the initial state.

##  3 Gravitational production of RHNs

In this section, we compute the gravitational interactions and resulting
abundances of the dark matter candidate, $N_{1}$ as well as the abundance of
the RHN neutrino responsible for leptogenesis. We are particularly interested
in the interactions of the type in Fig. 1 between the inflaton and the
$N_{i}$. In addition, we are interested in the gravitational interactions
between the inflaton and SM particles which make up the thermal bath.
Furthermore, we will show that it is possible to produce the thermal bath
assuming the absence of any inflaton decay mode leading to reheating. We will
quantify how such interactions can give rise to reasonable relic density and
baryon asymmetry. The DM candidate $N_{1}$ can be produced during reheating
from inflaton scattering $\phi\phi\rightarrow N_{1}N_{1}$ as well as from the
thermal bath (mediated by a massless graviton in both cases). For
$(y_{N})_{1i}=0$, $N_{1}$ couples only gravitationally, and the SM processes
will not lead to its production. On the other hand, for the generation of the
baryon asymmetry, we will cater to non-thermal leptogenesis, where the RHNs
$N_{2,3}$ are too weakly coupled to reach thermal equilibrium. Hence they are
predominantly produced only during reheating from gravitational inflaton
scattering. To summarize, we consider the following production via graviton
exchange

  * •

$\phi\,\phi\to N_{1}\,N_{1}$, SM SM $\to N_{1}\,N_{1}$ for production of the
DM candidate $N_{1}$.

  * •

$\phi\,\phi\to N_{2,3}\,N_{2,3}$ for production of $N_{2,3}$ that will lead to
non-thermal leptogenesis. (Contributions from SM SM $\to N_{2,3}\,N_{2,3}$ are
negligibly small.)

  * •

$\phi\phi\rightarrow$ SM SM for the reheating process.

In the following subsections we determine the production rate of DM and RHNs
from the processes listed above.

###  3.1 Gravitational dark matter

We start by computing the DM number density via 2-to-2 scattering of the bath
particles, mediated by graviton exchange. In this case the interaction rate is
given by [44, 85, 66, 68]

$\displaystyle
R_{N_{i}}^{T}=\frac{1}{2}\times\beta_{1/2}\,\frac{T^{8}}{M_{P}^{4}}\,,$ (3.1)

where we have $\beta_{1/2}=11351\pi^{3}/10368000\simeq 3.4\times 10^{-2}$
[68], the explicit factor of $\frac{1}{2}$ accounting for the Majorana nature
of $N_{i}$. The evolution of RHN number density $n_{N_{i}}$ (with $i=1,2,3$)
is governed by the Boltzmann equation

$\displaystyle\frac{dn_{N_{i}}}{dt}+3\,H\,n_{N_{i}}=R_{N_{i}}^{T}\,,$ (3.2)

where $H=\dot{a}/a$ is the Hubble parameter. Defining the comoving number
density as $Y_{N_{i}}=n_{N_{i}}\,a^{3}$ we can re-cast the Boltzmann equation
as

$\displaystyle\frac{dY^{T}_{N_{i}}}{da}=\frac{a^{2}}{H}\,R_{N_{i}}^{T}\,,$
(3.3)

where $i=1$ for DM production, and the superscript $T$ refer to the thermal
source of production. In order to properly capture the evolution of the
inflaton and radiation energy density (and hence temperature) we solve the
following set of coupled equations

$\displaystyle\frac{d\rho_{\phi}}{dt}+3H\,(1+w_{\phi})\,\rho_{\phi}=-(1+w_{\phi})\,\Gamma_{\phi}\,\rho_{\phi}\,,$
$\displaystyle\frac{d\rho_{R}}{dt}+4H\,\rho_{R}=+(1+w_{\phi})\,\Gamma_{\phi}\,\rho_{\phi}\,,$
(3.4)

where $w_{\phi}\equiv\frac{p_{\phi}}{\rho_{\phi}}=\frac{k-2}{k+2}$ [52] is the
general equation of state parameter, $\rho_{\phi}$ the inflaton energy
density, $\rho_{R}=\frac{\pi^{2}g_{*}}{30}\,T^{4}\equiv c_{*}\,T^{4}$, and
$g_{*}$ is the effective number of degrees of freedom for the thermal plasma
at the temperature $T$. We recall that in our framework the potential
$V(\phi)$ is proportional to $\phi^{k}$ [see, Eq. (2.6)].

During reheating, the total energy density is dominated by the inflaton and we
can approximate the Hubble parameter by $H^{2}\simeq\rho_{\phi}/3\,M_{P}^{2}$.
In this case, it is possible to analytically solve Eq. (3.1) and obtain

$\rho_{\phi}(a)=\rho_{\rm end}\left(\frac{a_{\rm
end}}{a}\right)^{\frac{6k}{k+2}}\,.$ (3.5)

Recall that we are assuming that the radiation bath is produced
gravitationally through inflaton scattering; namely, we do not rely on a
specific decay channel $\phi\rightarrow$ SM particles for reheating. In this
case, due to helicity conservation, the production of SM fermions from
inflaton scattering is strongly suppressed by the mass of the fermions,
whereas massless vectors are not produced because of the antisymmetry of
$F_{\mu\nu}$. However, scattering into scalars, especially Higgs scalars, is
always allowed and dominates the production rate. In [52], the inflaton
dissipation rate was parameterized as $\Gamma_{\phi}\propto\rho_{\phi}^{l}$.
For a quartic interaction with constant coupling, $l=(3/k)-(1/2)$. However,
for the effective gravitational coupling between the inflaton and SM Higgs,
the coupling is proportional to $m_{\phi}^{2}\propto\rho_{\phi}^{(1-2/k)}$.
This leads to $l=(3/2)-(1/k)$. More accurately, expanding the potential energy
in terms of the Fourier modes [104, 105, 68, 84, 52, 70]

$V(\phi)=V(\phi_{0})\sum_{n=-\infty}^{\infty}{\cal P}_{n}^{k}e^{-in\omega
t}=\rho_{\phi}\sum_{n=-\infty}^{\infty}{\cal P}_{n}^{k}e^{-in\omega t}\,,$
(3.6)

the production rate of radiation is given by [68, 69, 84]

$(1+w_{\phi})\,\Gamma_{\phi}\,\rho_{\phi}=R^{\phi^{k}}_{H}\simeq\frac{N_{h}\rho_{\phi}^{2}}{16\pi
M_{P}^{4}}\sum_{n=1}^{\infty}2n\omega|{\mathcal{P}}^{k}_{2n}|^{2}=\alpha_{k}M_{P}^{5}\left(\frac{\rho_{\phi}}{M_{P}^{4}}\right)^{\frac{5k-2}{2k}}\,,$
(3.7)

where $N_{h}=4$ is the number of internal degrees of freedom for one complex
Higgs doublet and we have neglected the Higgs bosons mass. The frequency of
oscillations of $\phi$ is given by [52]

$\omega=m_{\phi}\sqrt{\frac{\pi
k}{2(k-1)}}\frac{\Gamma(\frac{1}{2}+\frac{1}{k})}{\Gamma(\frac{1}{k})}\,,$
(3.8)

with $m_{\phi}^{2}=\frac{\partial^{2}V(\phi)}{\partial\phi^{2}}$ being the
inflaton mass squared. The definition of $\alpha_{k}$ follows the analysis in
[84], with the values given in Table 1. For $l=(3/2)-(1/k)$ the results of
[52] yields for the evolution of the radiation density

$\displaystyle\rho_{R}(a)=$
$\displaystyle\rho_{\text{RH}}\,\left(\frac{a_{\text{RH}}}{a}\right)^{4}\,\left[\frac{1-(a_{\text{end}}/a)^{\frac{8k-14}{k+2}}}{1-(a_{\text{end}}/a_{\text{RH}})^{\frac{8k-14}{k+2}}}\right]\,,$
(3.9)

which can be obtained by combining Eqs.(3.1), (3.5) and (3.7). The evolution
in Eq. (3.9) is valid when $a_{\text{end}}\ll a\ll a_{\text{RH}}$ where
$a_{\text{end}}$ marks the end of inflation (or the onset of reheating), while
$a_{\text{RH}}$ indicates the end of reheating defined as
$\rho_{\phi}(a_{\text{RH}})=\rho_{R}(a_{\text{RH}})=\rho_{\text{RH}}$. To
obtain Eq. (3.5), we have supposed $H\gg\Gamma_{\phi}$, which is valid for all
$a$ because $\Gamma_{\phi}<H$ at the end of inflation and $\Gamma_{\phi}$
decreases faster than $H$ for all $k\geq 2$. As a result, we note that
gravitational reheating is only possible for $\frac{6k}{k+2}>4$ [cf. Eq.
(3.5)], i.e., when $\rho_{\phi}$ redshifts faster than $\rho_{R}$. This limits
our parameter space to $k>4$. It is also important to ensure that a
sufficiently large reheating temperature is attained to allow big bang
nucleosynthesis to occur at $T\sim 1$ MeV as we will discuss in more detail in
Sec. 3.3.

$k$ | $\alpha_{k}$ | $\alpha_{k}^{\xi}$
---|---|---
6 | 0.000193 | $\alpha_{k}+0.00766\,\xi_{h}^{2}$
8 | 0.000528 | $\alpha_{k}+0.0205\,\xi_{h}^{2}$
10 | 0.000966 | $\alpha_{k}+0.0367\,\xi_{h}^{2}$
12 | 0.00144 | $\alpha_{k}+0.0537\,\xi_{h}^{2}$
14 | 0.00192 | $\alpha_{k}+0.0702\,\xi_{h}^{2}$
16 | 0.00239 | $\alpha_{k}+0.0855\,\xi_{h}^{2}$
18 | 0.00282 | $\alpha_{k}+0.0995\,\xi_{h}^{2}$
20 | 0.00322 | $\alpha_{k}+0.112\,\xi_{h}^{2}$
Table 1: Relevant coefficients for the gravitational reheating [cf. Eq. (3.7)
and Eq. (3.50)].

Using Eqs. (3.1) and (3.9) and relating $T^{8}$ to $\rho_{R}^{2}$, we obtain
the thermal rate of DM production

$\displaystyle
R^{T}_{N_{i}}=\frac{1}{2}\times\beta_{1/2}\,\frac{\rho_{\text{RH}}^{2}}{c_{*}^{2}\,M_{P}^{4}}$
$\displaystyle\left(\frac{a_{\text{RH}}}{a}\right)^{8}\,\left[\frac{1-(a_{\text{end}}/a)^{\frac{8k-14}{k+2}}}{1-(a_{\text{end}}/a_{\text{RH}})^{\frac{8k-14}{k+2}}}\right]^{2}\,.$
(3.10)

The DM number density at the end of reheating can then be computed by
integrating Eq. (3.3), leading to

$\displaystyle n^{T}_{N_{i}}(a_{\rm
RH})=\frac{\beta_{1/2}\,(k+2)\,\rho_{\text{RH}}^{\frac{3}{2}}}{12\,\sqrt{3}M_{P}^{3}\,c_{*}^{2}}\left(\frac{1}{1-r^{\frac{14-8k}{k+2}}}\right)^{2}\,$
$\displaystyle\left[\frac{2(7-4k)^{2}}{(k+5)(k-1)(5k-2)}\,r^{\frac{10+2k}{k+2}}-\frac{9}{(k+5)}+\frac{18}{(5k-2)}\,r^{\frac{14-8k}{k+2}}-\frac{1}{(k-1)}\,r^{\frac{28-16k}{k+2}}\right]\,,$
(3.11)

where $r=a_{\text{RH}}/a_{\text{end}}$. Since the gravitational reheating
temperature is generally quite low as discussed in Sec. 3.3, we can consider
the limit $r\gg 1$ and the dominant term in the expression above is

$n^{T}_{N_{i}}(a_{\rm
RH})\simeq\frac{\beta_{1/2}\,(k+2)\,\rho_{\text{RH}}^{\frac{3}{2}}}{12\,\sqrt{3}M_{P}^{3}\,c_{*}^{2}}\frac{2(7-4k)^{2}}{(k+5)(k-1)(5k-2)}\,r^{\frac{10+2k}{k+2}}\,.$
(3.12)

The contribution of gravitational scattering of the particles in the
primordial plasma to the DM relic abundance can then be determined using [106]

$\displaystyle\Omega_{N_{1}}^{T}\,h^{2}=1.6\times 10^{8}\,\frac{g_{0}}{g_{\rm
RH}}\,\frac{n_{N_{1}}^{T}(a_{\text{RH}})}{T_{\text{RH}}^{3}}\,\frac{M_{N_{1}}}{\text{GeV}}\,,$
(3.13)

which gives

$\displaystyle\Omega_{N_{1}}^{T}\,h^{2}$ $\displaystyle\simeq 1.6\times
10^{8}\times\frac{g_{0}\,\beta_{1/2}}{g_{\rm RH}}\,\times\frac{M_{N_{1}}}{\rm
GeV}\,\frac{c_{*}^{-\frac{5}{6}-\frac{5}{3k}}\,(7-4k)^{2}\,(k+2)}{6\,\sqrt{3}\,(k+5)(k-1)(5k-2)}\left(\frac{T_{\text{RH}}}{M_{P}}\right)^{\frac{5k-20}{3k}}\left(\frac{\rho_{\rm
end}}{M_{P}^{4}}\right)^{\frac{k+5}{3k}}\,,$ (3.14)

where $g_{0}=g_{*}(T_{0})=43/11$ and $g_{\rm RH}=g_{*}(T_{\text{RH}})=427/4$
are the number of relativistic degrees of freedom at present and at the end of
reheating respectively. In (3.14), we have used Eq. (3.5) evaluated at
$a=a_{\rm RH}$ to relate $r$, $T_{\text{RH}}$, and $\rho_{\rm end}$.

The DM candidate $N_{1}$ can also be produced directly from inflaton
scattering. In fact, particle production directly from inflaton scattering can
be much more efficient than gravitational thermal scattering when
$T_{\text{RH}}$ is low [68]. As we shall see in the next section, the same
process is also involved in the production of the baryon asymmetry. For the
production of $N_{1}$ through the scattering of the inflaton condensate, we
consider the time dependent oscillations of a classical inflaton field
$\phi(t)$555As it has been pointed out in [107, 108, 109], for all potentials
steeper than quadratic near the origin, the oscillating inflaton condensate
may fragment due to self-resonance. The equation of state in that case
approaches $w\to 1/3$ at sufficiently late times. If this occurs after
$T=T_{\text{RH}}$, then this effect is not important since the inflaton energy
would already be subdominant. Furthermore, $N$ is predominantly produced at
$a_{\text{end}}$ and so possible fragmentation at later times would not affect
the calculation of the baryon asymmetry. Reheating may be affected; however,
the exact time when $w$ transitions $1/3$ depends on $k$ and requires
dedicated lattice simulations. Self-resonance becomes less efficient for
larger $k$ as shown by the lattice results for $k$ up to 6 [107, 108, 109],
while most of our viable results are for $k>6$. Performing such lattice
simulations for larger $k$ is beyond the scope of the present analysis.. The
oscillating inflaton field with a time-dependent amplitude can be parametrized
as

$\phi(t)=\phi_{0}(t)\cdot\mathcal{Q}(t)=\phi_{0}(t)\sum_{n=-\infty}^{\infty}\,{\cal
Q}_{n}e^{-in\omega t}\,,$ (3.15)

where $\phi_{0}(t)$ is the time-dependent amplitude that includes the effects
of redshift and $\mathcal{Q}(t)$ describes the periodicity of the oscillation.
Furthermore, we assume a mass hierarchy $M_{N_{1,2,3}}<m_{\phi}$ such that the
$s$-channel graviton mediated process (as shown in Fig. 1) is kinematically
allowed. Note that, since $N_{1}$ is effectively decoupled from $N_{2,3}$, it
does not necessarily need to be the lightest of the three. The production rate
for $N_{i}$ from inflaton scattering mediated by gravity is given666Note the
difference of factor 2 with [68], comes from the Majorana nature of the RHNs.
by [68]

$R^{\phi^{k}}_{N_{i}}=\frac{\rho_{\phi}^{2}}{4\pi
M_{P}^{4}}\,\frac{M_{N_{i}}^{2}}{m_{\phi}^{2}}\,\Sigma_{N_{i}}^{k}\,,$ (3.16)

where

$\Sigma_{N_{i}}^{k}=\sum_{n=1}^{+\infty}|{\cal
P}_{2n}^{k}|^{2}\,\frac{m_{\phi}^{2}}{E_{2n}^{2}}\,\left[1-\frac{4\,M_{N_{i}}^{2}}{E_{2n}^{2}}\right]^{3/2}\,,$
(3.17)

accounts for the sum over the Fourier modes of the inflaton potential, and
$m_{\phi}^{2}=\lambda\,k\,(k-1)\,\left(\rho_{\phi}/(\lambda\,M_{P}^{4})\right)^{\frac{k-2}{k}}$.
Here $E_{n}=n\omega$ is the energy of the $n$-th inflaton oscillation mode and
$M_{N_{i}}$ is the mass of the produced RHN. Then, the number density of RHN
is obtained by solving a Boltzmann equation analogous to that in Eq. (3.3) as

$\frac{dY^{\phi^{k}}_{N_{i}}}{da}=\frac{\sqrt{3}M_{P}}{\sqrt{\rho_{\rm
RH}}}a^{2}\left(\frac{a}{a_{\rm
RH}}\right)^{\frac{3k}{k+2}}R_{N_{i}}^{\phi^{k}}(a).$ (3.18)

Integration of Eq. (3.18), leads to the following expression for the RHN
density [68, 84]

$n_{N_{i}}^{\phi^{k}}(a_{\rm
RH})\simeq\frac{M_{N_{1}}^{2}\,\sqrt{3}\,(k+2)\,\rho_{\rm
RH}^{\frac{1}{2}+\frac{2}{k}}}{24\,\pi\,k(k-1)\lambda^{\frac{2}{k}}\,M_{P}^{1+\frac{8}{k}}}\left(\frac{\rho_{\rm
end}}{\rho_{\rm RH}}\right)^{\frac{1}{k}}\,\Sigma_{N_{1}}^{k}\,,$ (3.19)

evaluated at the end of reheating. In order to obtain the DM relic abundance,
one can again follow Eq. (3.13), but now replacing
$n_{N_{1}}^{T}(a_{\text{RH}})$ with $n_{N_{1}}^{\phi}(a_{\text{RH}})$, and
obtain [68]

$\displaystyle\frac{\Omega_{N_{1}}^{\phi^{k}}h^{2}}{0.12}$ $\displaystyle=$
$\displaystyle\frac{\Sigma_{N_{1}}^{k}}{2.4^{\frac{8}{k}}}\frac{k+2}{k(k-1)}\left(\frac{10^{-11}}{\lambda}\right)^{\frac{2}{k}}\left(\frac{10^{40}{\rm
GeV}^{4}}{\rho_{\rm RH}}\right)^{\frac{1}{4}-\frac{1}{k}}$ (3.20)
$\displaystyle\times$ $\displaystyle\left(\frac{\rm\rho_{end}}{10^{64}{\rm
GeV}^{4}}\right)^{\frac{1}{k}}\left(\frac{M_{N_{1}}}{1.1\times
10^{7^{+}\frac{6}{k}}{\rm GeV}}\right)^{3}\,.$

The total DM relic abundance is a sum of the gravitational contribution from
thermal bath $(\Omega_{N_{1}}^{T}\,h^{2})$ and from inflaton scattering
$(\Omega_{N_{1}}^{\phi^{k}}\,h^{2})$.

###  3.2 Gravitational leptogenesis

Since $N_{1}$ is the stable DM candidate, in the present scenario the lighter
of $N_{2,3}$ can undergo out-of-equilibrium decay to SM final states. We
denote $N_{2}$ to be the lighter of these, and we must require that the mixing
of $N_{1}$ and $N_{2}$ to be sufficiently small so as to prevent the decay of
$N_{2}$ to $N_{1}$. For now, we take this coupling to be absent. The resulting
CP asymmetry from the decay of $N_{2}$ is of the form [110, 111, 112, 113,
114]

$\displaystyle\epsilon_{\Delta L}=\frac{\sum_{\alpha}[\Gamma(N_{2}\rightarrow
l_{\alpha}+H)-\Gamma(N_{2}\rightarrow\overline{l}_{\alpha}+H^{*})]}{\sum_{\alpha}[\Gamma(N_{2}\rightarrow
l_{\alpha}+H)+\Gamma(N_{2}\rightarrow\overline{l}_{\alpha}+H^{*})]}\,.$ (3.21)

The resulting lepton asymmetry depends on the out-of-equilibrium abundance of
$N_{2}$ as computed in the previous subsection. So long as $M_{N_{2}}\ll
m_{\phi}$ and any kinematic suppression can be ignored, the number density of
$N_{2}$ (at $a_{\rm RH}$) will be given by Eq. (3.19) with the substitution
$N_{1}\to N_{2}$. The CP asymmetry can be expressed as [113, 43, 84]

$\displaystyle\epsilon_{\Delta L}\simeq\frac{3\delta_{\rm
eff}}{16\,\pi}\,\frac{M_{N_{2}}\,m_{\nu\,,\text{max}}}{v^{2}}\,,$ (3.22)

where $\langle H\rangle\equiv v\approx 174$ GeV is the SM Higgs doublet vacuum
expectation value, $\delta_{\text{eff}}$ is the effective CP violating phase
in the neutrino mass matrix with $0\leq\delta_{\rm eff}\leq 1$, and, we take
$m_{\nu,\text{max}}=0.05$ eV as the heaviest light neutrino mass.

Here we are interested in non-thermal leptogenesis [79, 81, 82, 115, 116, 83,
117]. The gravitationally produced $N_{2}$ should not be thermalized into the
bath for the consistency of the calculation. To check this, we note that the
thermalization rate $\Gamma_{\rm th}\simeq y_{N_{2}}^{2}T/8\pi$ decreases
slower than the Hubble rate $H=\frac{\sqrt{\rho_{\phi}}}{\sqrt{3}M_{P}}$ based
on Eqs. (3.5) and (3.9). Thermalization is potentially dangerous until
$T\simeq M_{N_{2}}$ when the $N_{2}$ out-of-equilibrium decay rate dominates
over the thermalization rate. Using Eq. (3.5), $a_{\text{end}}/a\simeq
T/T_{\text{max}}$ based on Eq. (3.9), and $y_{N_{2}}\simeq
m_{\nu}M_{N_{2}}/v^{2}$, we find that $\Gamma_{\rm th}$ is always less than
$H$ at $T=M_{N_{2}}$ in the parameter space of interest.777More precisely,
$T/T_{\text{max}}=a_{\rm max}/a=(a_{\rm
max}/a_{\text{end}})(a_{\text{end}}/a)=((6k-3)/(2k+4))^{((k+2)/(8k-14))}(a_{\text{end}}/a)\approx
3^{1/8}(a_{\text{end}}/a)$ for large $k$ [68]. Thus the resulting lepton
asymmetry will not be suppressed by inverse decays.

The produced lepton asymmetry is eventually converted to baryon asymmetry via
electroweak sphaleron processes leading to

$\displaystyle Y_{B}=\frac{n_{B}}{s}=\frac{28}{79}\epsilon_{\Delta
L}\,\frac{n^{\phi}_{N_{2}}(T_{\text{RH}})}{s}\,,$ (3.23)

where $n_{N_{2}}^{\phi}(T_{\text{RH}})$ is the number density from Eqs. (3.12)
and (3.19) at the end of reheating and $s=2\pi^{2}g_{\rm
RH}T_{\text{RH}}^{3}/45$ is the entropy density. The final asymmetry then
becomes

$\displaystyle Y_{B}\simeq 3.5\times 10^{-4}\,\delta_{\rm
eff}\left(\frac{m_{\nu,\text{max}}}{0.05\,\text{eV}}\right)\left(\frac{M_{N_{2}}}{10^{13}\,\text{GeV}}\right)\,\left.\frac{n_{N_{2}}^{\phi}}{s}\right|_{T_{\text{RH}}}\,,$
(3.24)

while the observed value, as reported by Planck [3], is
$Y_{B}^{\text{obs}}\simeq 8.7\times 10^{-11}$.

We note that the lepton asymmetry is not washed out because the lepton-number
violating process involving the Yukawa scattering and the electroweak
sphaleron processes are never in equilibrium at the same time.

###  3.3 Gravitational reheating temperature

In section 3.1, we computed the energy density in radiation from purely
gravitational process. However, to avoid conflict with the Big Bang
nucleosynthesis (BBN), that requires the reheating temperature
$T_{\text{RH}}\gtrsim 1$ MeV, one needs to consider888This requirement of
having large $w_{\phi}$ can be relaxed with non-minimal gravitational
couplings as discussed in [84, 69]. $w_{\phi}\gtrsim 0.65$ [86, 84], or
$k=\frac{2\,(1+w_{\phi})}{(1-w_{\phi})}\gtrsim 9$. This lower bound comes from
the fact that, for higher $k$, the inflaton energy density redshifts faster so
the transition to radiation domination is achieved sooner, at a higher
temperature.

The precise bound on $T_{\text{RH}}$ is in fact more involved especially for
$k\leq 8$ for the following reason. As noted below Eq. (3.9), the inflaton-
dominated era ends when $\rho_{\phi}$ redshifts below $\rho_{R}$ as opposed to
a complete transfer of the $\phi$ energy to radiation in conventional
reheating by perturbative decays. Still, we define the reheating temperature
$T_{\text{RH}}$ by $\rho_{\phi}(T_{\text{RH}})=\rho_{R}(T_{\text{RH}})$. The
difference in the scale factor dependence between $\rho_{\phi}$ in Eq. (3.5)
and $\rho_{R}$ in Eq. (3.9) increases with $k$. In other words, for smaller
$k$, the inflaton energy density does not redshift significantly more than
radiation. Thus, $T_{\text{RH}}$ for low $k$ needs to be substantially higher
than $T_{\rm BBN}\approx 1$ MeV so that the inflaton energy density does not
excessively modify the expansion rate of the universe at BBN. We can recast
the BBN bound on the extra energy density in the form of $\Delta
N_{\nu}<0.226$ [118] into a bound on $T_{\text{RH}}$ as a function of $k$
using the following expression

$\displaystyle
1=\left.\frac{\rho_{\phi}}{\rho_{R}}\right|_{T_{\text{RH}}}=\left.\frac{\rho_{\phi}}{\rho_{R}}\right|_{T_{\rm
BBN}}\frac{(a_{\rm BBN}/a_{\text{RH}})^{\frac{6k}{k+2}}}{(T_{\text{RH}}/T_{\rm
BBN})^{4}},$ (3.25)

with entropy conservation
$g_{*}(a_{\text{RH}})\,T_{\text{RH}}^{3}\,a_{\text{RH}}^{3}=g_{*}(a_{\rm
BBN})\,T_{\rm BBN}^{3}\,a_{\rm BBN}^{3}$ within the SM sector, as well as
$\rho_{\phi}(T_{\rm
BBN})=\frac{7}{8}\,\left(\frac{4}{11}\right)^{\frac{4}{3}}\,\Delta
N_{\nu}\,\frac{\pi^{2}}{30}\,T_{\rm BBN}^{4}$. The resulting bounds for each
$k$ are $T_{\text{RH}}\lesssim 40\,{\rm MeV}~{}(k=6),\,9\,{\rm
MeV}~{}(k=8),\,5\,{\rm MeV}~{}(k=10),\,4\,{\rm MeV}~{}(k=12)$ and $3\,{\rm
MeV}~{}(14\leq k\leq 20)$, which we will show as red-shaded regions in the
subsequent figures concerning $T_{\text{RH}}$. We note that this estimate is
in good agreement with a more rigorous treatment performed in Ref. [119] using
a BBN computing package for the case of kination (large $k$ limit).

The reheating temperature can be determined by solving the Friedmann equation
(3.1) for the radiation energy density. This yields [84]

$\displaystyle\rho_{R}(a)$
$\displaystyle\simeq\alpha_{k}\frac{k+2}{8k-14}\sqrt{3}M_{P}^{4}\left(\frac{\rho_{\rm
end}}{M_{P}^{4}}\right)^{\frac{2k-1}{k}}\left(\frac{a_{\text{end}}}{a}\right)^{4}\,,$
(3.26)

and evaluating this at $a_{\rm RH}$ we have

$\displaystyle T_{\text{RH}}^{4}=\frac{30}{\pi^{2}\,g_{\rm
RH}}\,M_{P}^{4}\,\left(\frac{\rho_{\text{end}}}{M_{P}^{4}}\right)^{\frac{4k-7}{k-4}}\,\left(\frac{\alpha_{k}\,\sqrt{3}\,(k+2)}{8k-14}\right)^{\frac{3k}{k-4}}\,.$
(3.27)

From Eq. (3.27) we find $T_{\text{RH}}\simeq 60$ MeV for $k=10$ and
$\rho_{\text{end}}\simeq(4.8\times 10^{15}\,\text{GeV})^{4}$. Note that, due
to the logarithmic dependence of $\phi_{\rm end}$ on $k$ in Eq. (2.7),
$\rho_{\text{end}}$ changes very slowly with $k$ and remains approximately
fixed to the above value.

In Fig. 2, we show in the left panel the reheating temperature for minimal
gravitational interactions by the curve labeled $\xi_{h}=0$ (other values of
$\xi_{h}$, non-minimal coupling of the Higgs to the Ricci scalar, are
discussed in the next subsection). As one can see, $T_{\text{RH}}$ rises to
$\simeq 1$ TeV, for $k=20$. This minimal case with $\xi_{h}=0$ is excluded by
excessive gravitational waves as dark radiation as will be discussed in Sec.
3.5, so a non-minimal coupling is ultimately required.

Figure 2: Variation of $T_{\text{RH}}$ (left) and $T_{\text{max}}$ (right) as
a function of $k$, for different choices of $\xi_{h}$. Triangles highlight the
physical choices of even $k$. The red-shaded region is excluded by BBN because
low reheating temperatures lead to an excessive inflaton energy density during
BBN. The blue-shaded region is similarly excluded by BBN for excessive
gravitational waves produced during inflation. The gray-shaded region is
excluded as the lowest reheating temperature from gravitational reheating is
that from minimal gravity (pure graviton exchange), i.e., $\xi_{h}=0$.

At the start of the reheating process, the Universe quickly heats to a maximum
temperature, $T_{\rm max}$. As discussed in [68], the maximum temperature
attained through purely gravitational processes is of order $10^{12}$ GeV,
decreasing slightly with $k$. The maximum radiation density which determines
$T_{\text{max}}$ was found to be [68]

$\displaystyle\rho_{\text{max}}\simeq\sqrt{3}\,\alpha_{k}\,M_{P}^{4}\,\left(\frac{\rho_{\rm
end}}{M_{P}^{4}}\right)^{\frac{2k-1}{k}}\,\frac{k+2}{12k-16}\,\left(\frac{2k+4}{6k-3}\right)^{\frac{2k+4}{4k-7}}\equiv
c_{*}\,T_{\rm max}^{4}\,.$ (3.28)

Asymptotically at large $k$, $T_{\rm max}\approx 8\times 10^{11}$ GeV and
$\left(\frac{T_{\text{max}}}{T_{\text{RH}}}\right)_{k\gg
4}\sim\left(\frac{1}{\alpha_{k\gg
4}}\right)^{1/2}\left(\frac{M_{P}^{4}}{\rho_{\rm end}}\right)^{1/2}\gg 1$.
This represents a minimal maximum temperature, as other process such as decays
(not considered here), may lead to a higher maximum temperature. The value of
$T_{\text{max}}$ is shown in the right panel of Fig. 2. For minimal
gravitational interactions, corresponding to the simple exchange of a
graviton, $T_{\text{max}}\simeq 10^{12}\,{\rm GeV}$. In fact, as we noted
before and will elaborate in Sec. 3.5, gravitational reheating with
$\xi_{h}=0$ (graviton exchange) is already ruled out by the BBN bound on dark
radiation in the form of GWs. Thus, in order to account for the reheating
mechanism in a gravitational framework, it is necessary to introduce non-
minimal couplings of fields to gravity. We compute the reheating and maximum
temperatures for non-minimal gravitational interactions in the next
subsection. The value of $T_{\text{max}}$ will be relevant when we discuss the
DM warmness constraint, because DM is produced relativistically and
predominantly at $T_{\text{max}}$.

###  3.4 Non-minimal gravitational production

As pure gravitational particle production can be insufficient, we also
consider the possibility that scalar fields have non-minimal couplings to
gravity which generate effective couplings between these scalar fields and the
RHNs. Thus, we allow both the inflaton $\phi$ and the Higgs field $H$ to be
non-minimally coupled. We denote the complex Higgs doublet as $h$ throughout
the following section. One can then write the action as

$\mathcal{S}_{J}=\int\,d^{4}x\sqrt{-\tilde{g}}\left[-\frac{M_{P}^{2}}{2}\,\Omega^{2}\,\widetilde{\mathcal{R}}+\widetilde{\mathcal{L}}_{\phi}+\widetilde{\mathcal{L}}_{h}+\widetilde{\mathcal{L}}_{N_{i}}\right]\,,$
(3.29)

where

$\displaystyle\widetilde{\mathcal{L}}_{\phi}=\frac{1}{2}\,\partial_{\mu}\phi\,\partial^{\mu}\phi-V(\phi)$
(3.30)
$\displaystyle\widetilde{\mathcal{L}}_{h}=\partial_{\mu}h\,\partial^{\mu}h^{\dagger}-V(hh^{\dagger})$
(3.31)
$\displaystyle\widetilde{\mathcal{L}}_{N_{i}}=\frac{i}{2}\,\overline{\mathcal{N}_{i}}\,\overleftrightarrow{\not{\nabla}}\,\mathcal{N}_{i}-\frac{1}{2}\,M_{N_{i}}\,\overline{(\mathcal{N})^{c}}_{i}\,\mathcal{N}_{i}+\widetilde{\mathcal{L}}_{\text{yuk}}$
(3.32)
$\displaystyle\widetilde{\mathcal{L}}_{\text{yuk}}=-y_{N_{i}}\,\overline{\mathcal{N}_{i}}\,\widetilde{h^{\dagger}}\,\mathbb{L}+\text{h.c.}\,,$
(3.33)

where $\mathcal{N}\,,\mathbb{L}$ are the RHN and SM lepton doublet fields in
Jordan frame. The conformal factor $\Omega^{2}$ is given by

$\Omega^{2}\equiv
1+\frac{\xi_{\phi}\,\phi^{2}}{M_{P}^{2}}+\frac{\xi_{h}\,|h|^{2}}{M_{P}^{2}}\,.$
(3.34)

It is convenient to remove the non-minimal couplings by performing the
redefinition of the metric field via the usual conformal transformation from
the Jordan frame to the Einstein frame,

$g_{\mu\nu}=\Omega^{2}\,\tilde{g}_{\mu\nu}\,,$ (3.35)
$\displaystyle\mathcal{S}_{E}=\int
d^{4}x\sqrt{-g}\Biggl{[}-\frac{M_{P}^{2}\,\mathcal{R}}{2}+\frac{K^{ab}}{2}\,g^{\mu\nu}\partial_{\mu}S_{a}\,\partial_{\nu}S_{b}+\frac{i}{2\,\Omega^{3}}\,\overline{N_{i}}\,\overleftrightarrow{\not{\nabla}}\,N_{i}-\frac{1}{\Omega^{4}}\,\left(\frac{M_{N_{i}}}{2}\,\overline{N_{i}^{c}}\,N_{i}+\mathcal{L}_{\text{yuk}}\right)$
$\displaystyle-\frac{3i}{4\Omega^{4}}\,\overline{N_{i}}\,\left(\overleftrightarrow{\not{\partial}}\,\Omega\right)\,N_{i}-\frac{1}{\Omega^{4}}\,\left(V_{\phi}+V_{h}\right)\Biggr{]}\,,$
(3.36)

where we have used

$\sqrt{-\tilde{g}}\rightarrow\frac{\sqrt{-g}}{\Omega^{4}}$ (3.37)
$\tilde{\not{\nabla}}\rightarrow\Omega\not{\nabla}-\frac{3}{2}\Omega^{2}(\not{\partial}\Omega)\,,$
(3.38)

and the indices $a,b$ enumerate the fields $\phi$, and the real components of
$h$. Then making spinor field redefinition $L\rightarrow\Omega^{3/2}L$,
$N_{i}\rightarrow\Omega^{3/2}N_{i}$ and
$\overline{N_{i}}\rightarrow\Omega^{3/2}\overline{N_{i}}$ we recover the
following action with canonical kinetic term for the RHN

$\displaystyle\mathcal{S}_{E}=\int
d^{4}x\sqrt{-g}\Biggl{[}-\frac{M_{P}^{2}\,\mathcal{R}}{2}+\frac{K^{ab}}{2}\,g^{\mu\nu}\partial_{\mu}S_{a}\,\partial_{\nu}S_{b}-\frac{1}{\Omega^{4}}\,\left(V_{\phi}+V_{h}\right)+\frac{i}{2}\,\overline{N_{i}}\,\overleftrightarrow{\not{\nabla}}\,N_{i}$
$\displaystyle-\frac{1}{2\,\Omega}\,M_{N_{i}}\,\overline{N_{i}^{c}}\,N_{i}+\frac{1}{\Omega}\,\mathcal{L}_{\text{yuk}}\Biggr{]}\,.$
(3.39)

The kinetic function is given by

$K^{ab}=6\,\frac{\partial\log\Omega}{\partial
S_{a}}\frac{\partial\log\Omega}{\partial
S_{b}}+\frac{\delta^{ab}}{\Omega^{2}}\,.$ (3.40)

In what follows, we will be interested in the small-field limit

$\frac{|\xi_{\phi}|\phi^{2}}{M_{P}^{2}},\,\frac{|\xi_{h}||h|^{2}}{M_{P}^{2}}\ll
1\,.$ (3.41)

Since we necessarily consider large field values for the inflaton coming out
of inflation, this constrains $\xi_{\phi}\ll 1$. However, there is no such
constraint on $\xi_{h}$ which can take relatively large values.

Expanding the kinetic and potential terms in the action in Eq. (3.4) in powers
of $1/M_{P}^{2}$, we obtain a canonical kinetic term for the scalar fields,
and deduce the leading-order interactions between scalars and the RHNs induced
by the non-minimal couplings. Note that, the non-minimally coupled Yukawa
interaction in Eq. (3.4) gives rise to a 3-to-2 (or 2-to-3) process with RHN
in the final state above electroweak symmetry breaking temperature. These
processes thus will be kinematically suppressed and have a subdominant
contribution to the RHN yield. The kinetic terms for the RHNs can be expressed
in the form

$\mathcal{L}_{\rm{non-
min.}}\;=\;-\sigma_{hN_{i}}^{\xi}\,|h|^{2}\,\overline{N_{i}^{c}}N_{i}-\sigma_{\phi
N_{i}}^{\xi}\,\phi^{2}\,\overline{N_{i}^{c}}N_{i}\,,$ (3.42)

with

$\sigma_{\phi N_{i}}^{\xi}\;=\;\frac{M_{N_{i}}}{2M_{P}^{2}}\xi_{\phi}$ (3.43)
$\sigma_{hN_{i}}^{\xi}\;=\;\frac{M_{N_{i}}}{2M_{P}^{2}}\xi_{h}\,.$ (3.44)

These non-minimal interactions open up additional channels [68]

  * •

RHN production from inflaton scattering: $\phi\phi\to N_{i}\,N_{i}$

  * •

RHN production from Higgs scattering: $hh^{\dagger}\to N_{i}\,N_{i}$

  * •

Higgs production from inflaton scattering: $\phi\phi\to hh^{\dagger}$.

Interestingly, as can be seen from the interaction terms, the production of
RHNs is systematically proportional to the mass of the fermion. Then, for the
thermal production of RHNs, the production rate is

$\displaystyle R^{T,\xi}_{N_{i}}\simeq
N_{h}\,\frac{\zeta(3)^{2}\,\xi_{h}^{2}}{32\,\pi^{5}}\,\frac{M_{N_{i}}^{2}\,T^{6}}{M_{P}^{4}}\,,$
(3.45)

where $\zeta(x)$ is the Riemann-zeta function. For both minimal and non-
minimal gravitational couplings, the leading term in the production rate for
scalar dark matter scales as $T^{8}$ [69]. Similarly, the production rate for
fermions in minimal gravity also scales as $T^{8}$ as seen in Eq. (3.1).
However, for non-minimal gravitational interactions, after the conformal
rescaling to obtain canonical kinetic terms, there is no non-minimal coupling
to the kinetic terms (in contrast to the scalars where this coupling is found
in Eq. (3.4)). Instead, we are left with only the mass term coupled to
$|h|^{2}$ and the thermal production rate is proportional only to $\propto
M_{N_{i}}^{2}T^{6}$.

Using the rate in Eq. (3.45) we obtain the number density at the end of
reheating due to the non-minimal interaction as

$\displaystyle n_{N_{i}}^{T(\xi_{h}\neq
0)}=\left(\frac{\sqrt{3}\,N_{h}\,\zeta(3)^{2}\,\xi_{h}^{2}}{32\,\pi^{5}\,c_{*}^{3/2}}\,\frac{M_{N_{i}}^{2}\,\rho_{\text{RH}}}{M_{P}^{3}}\right)\,\frac{(k+2)\,\left(1-\left(\frac{\rho_{\text{end}}}{\rho_{\text{RH}}}\right)^{\frac{7}{3\,k}-\frac{4}{3}}\right)^{-3/2}}{72\,(5-4k)\,\Gamma\left(\frac{29-20k}{14-8k}\right)}$
$\displaystyle\times\Biggl{[}9\sqrt{\pi}(5-4k)\left(\frac{\rho_{\text{end}}}{\rho_{\text{RH}}}\right)^{1/k}\Gamma\left(\frac{4k-4}{4k-7}\right)+4\,\left(\frac{\rho_{\text{end}}}{\rho_{\text{RH}}}\right)^{\frac{16k^{2}+4k+169}{21k-12k^{2}}}\,\Gamma\left(\frac{29-20k}{14-8k}\right)\,\mathcal{G}\Biggr{]}\,,$
(3.46)

with

$\displaystyle\mathcal{G}=\left(\frac{\rho_{\text{end}}}{\rho_{\text{RH}}}\right)^{\frac{4(k+30)}{3k(4k-7)}}\,\left[3(4k-5)\left(\frac{\rho_{\text{end}}}{\rho_{\text{RH}}}\right)^{\frac{16k^{2}+49}{3k(4k-7)}}-6\left(\frac{\rho_{\text{end}}}{\rho_{\text{RH}}}\right)^{\frac{56}{3(4k-7)}}\right]\,{}_{2}F_{1}(...)\,,$

where
${}_{2}F_{1}\left(-\frac{1}{2},\frac{3}{4k-7},\frac{4k-4}{4k-7},\left(\frac{\rho_{\text{end}}}{\rho_{\text{RH}}}\right)^{\frac{7}{3k}-\frac{4}{3}}\right)$
is the hypergeometric function.

For the inflaton scattering process $\phi\phi\rightarrow N_{i}\,N_{i}$, on the
other hand, we find

$\displaystyle R^{\phi,\xi}_{N_{i}}$ $\displaystyle=$
$\displaystyle\frac{M_{N_{i}}^{2}\xi_{\phi}^{2}\phi_{0}^{4}\omega^{2}}{32\pi
M_{P}^{4}}\sum\limits_{n=1}^{\infty}(2n)^{2}|\mathcal{Q}^{(2)}_{2n}|^{2}\times\sqrt{1-\frac{4M_{N_{i}}^{2}}{E_{2n}^{2}}}\,,$
(3.47)

where we define
$\phi_{0}=\left(\frac{\rho_{\phi}}{\lambda\,M_{P}^{4-k}}\right)^{\frac{1}{k}}$
and ${\cal Q}^{(2)}_{n}$ by

$\phi^{2}(t)=\phi_{0}^{2}(t)\cdot\mathcal{Q}^{2}(t)=\phi_{0}^{2}(t)\sum_{n=-\infty}^{\infty}\,{\cal
Q}^{(2)}_{n}e^{-in\omega t}\,.$ (3.48)

This rate is restricted by the small field limit that imposes a stringent
bound on $\xi_{\phi}$ from $\sqrt{|\xi_{\phi}|}\lesssim
M_{P}/\langle\phi\rangle$. In particular, at the beginning of inflaton
oscillations $\langle\phi\rangle\sim M_{P}$ and $|\xi_{\phi}|\lesssim 1$. When
we compare this rate with the one due to inflaton scattering mediated by
minimal gravitational interactions (Eq. (3.16)), we obtain

$\frac{R^{\phi,\xi}_{\phi
N_{i}}}{R^{\phi^{k}}_{N_{i}}}\approx\left[\frac{k\,(k-1)\,\xi_{\phi}\,(\omega/m_{\phi})^{2}}{\sqrt{8}}\right]^{2}\,\frac{\,\sum\limits_{n=1}^{\infty}(2n)^{2}|\mathcal{Q}^{(2)}_{2n}|^{2}}{\sum\limits_{n=1}^{\infty}\frac{1}{(2n)^{2}}\,\left|\mathcal{P}^{k}_{2n}\right|^{2}}\,.$
(3.49)

This ratio takes values between
$\\{32\,\xi_{\phi}^{2},\,242\,\xi_{\phi}^{2}\\}$ for $k\in[6,20]$. Hence, the
non-minimal contribution from inflaton scattering dominates over the graviton
exchange for $\xi_{\phi}>\frac{1}{2\sqrt{8}}$ when $k=6$ and for
$\xi_{\phi}>0.06$ when $k=20$. In what follows, we will neglect this
contribution as it dominates for values of $\xi_{\phi}$ close to the small
field limit, making the assumption of canonical kinetic terms of the fields
invalid.

The non-minimal coupling of Higgs bosons to gravity provides an additional
channel to reheat the Universe through gravitational processes, with the
following rate [84]

$(1+\omega_{\phi})\Gamma_{\phi}=R^{\phi,\xi}_{H}\simeq\frac{\xi_{h}^{2}N_{h}}{4\pi
M_{P}^{4}}\sum_{n=1}^{\infty}2n\omega\left|2\times{\mathcal{P}}^{k}_{2n}\rho_{\phi}+\frac{(n\omega)^{2}}{2}\phi_{0}^{2}|\mathcal{Q}_{n}|^{2}\right|^{2}=\alpha_{k}^{\xi}M_{P}^{5}\left(\frac{\rho_{\phi}}{M_{P}^{4}}\right)^{\frac{5k-2}{2k}}\,,$
(3.50)

where $\mathcal{Q}_{n}$ has been defined in Eq. (3.15) and $\alpha_{k}^{\xi}$
is given in Table 1. If we solve Eq. (3.1) for $\rho_{R}$, the reheating
temperature in the presence of the non-minimal coupling is then given by

$\displaystyle\left(T_{\text{RH}}^{\xi}\right)^{4}=\frac{30}{\pi^{2}\,g_{\rm
RH}}\,M_{P}^{4}\,\left(\frac{\rho_{\text{end}}}{M_{P}^{4}}\right)^{\frac{4k-7}{k-4}}\,\left(\frac{\alpha_{k}^{\xi}\,\sqrt{3}\,(k+2)}{8k-14}\right)^{\frac{3k}{k-4}}\,.$
(3.51)

The reheating temperature as a function of $k$ is shown in the left panel of
Fig. 2 for several values of $\xi_{h}$. The maximum temperature in this case
is determined from

$\displaystyle\rho_{\text{max}}^{\xi}\simeq\sqrt{3}\,\alpha_{k}^{\xi}\,M_{P}^{4}\,\left(\frac{\rho_{\text{end}}}{M_{P}^{4}}\right)^{\frac{2k-1}{k}}\,\frac{k+2}{12k-16}\,\left(\frac{2k+4}{6k-3}\right)^{\frac{2k+4}{4k-7}}\equiv
c_{*}\,(T_{\text{max}}^{\xi})^{4}\,.$ (3.52)

Contours of $T_{\text{max}}^{\xi}$ in the $(k,T_{\text{RH}})$ plane are shown
in the right panel of Fig. 2, where the appropriate values of $\xi_{h}$ taken
from the left panel are used to calculate $T_{\text{max}}^{\xi}$.

###  3.5 Gravitational waves generated during inflation

In this section, we review the calculation of gravitational waves generated by
quantum fluctuations during inflation, followed by a cosmological era where
the the inflaton energy dominates and redshifts faster than radiation. This
results in an enhancement of gravitational waves, which places a constraint
from excessive gravitational waves as dark radiation and offers a
gravitational wave signal with a distinctive spectrum.

The ratio of the gravitational wave (GW) energy density to that of the
radiation bath is given by [120]

$\displaystyle\frac{\rho_{\rm GW}}{\rho_{R}}=\frac{1}{32\pi
G\rho_{R}}\frac{k_{\rm GW}^{2}}{2}\mathcal{P}_{T}(k_{\rm GW})~{}~{}~{}{\rm
with}~{}~{}~{}\mathcal{P}_{T}(k_{\rm GW})\equiv\frac{2H_{I}^{2}(k_{\rm
GW})}{\pi^{2}M_{P}^{2}}\,,$ (3.53)

where $k_{\rm GW}$ is the momentum mode of the GW, $\mathcal{P}_{T}$ is the
primordial tensor power spectrum, $H_{I}(k_{\rm GW})$ is the Hubble scale
during inflation when the mode $k_{\rm GW}$ exits the horizon, $T_{\rm hc}$ is
the horizon-crossing temperature when the mode re-enters the horizon at
$k_{\rm GW}=H(T_{\rm hc})$, and the factor of $1/2$ accounts for the time
average of the rapidly oscillating metric perturbations. In our case,
$\rho_{\rm GW}/\rho_{R}$ is redshift invariant up to the change of $g_{*}$
after $T=T_{\text{max}}$ because entropy is only efficiently produced at
$T=T_{\text{max}}$ as discussed in Sec. 3.3. Therefore, the final
gravitational wave strength is given by $\Omega_{\rm
GW}h^{2}=\Omega_{\gamma}\,h^{2}\left(\rho_{\rm
GW}/\rho_{R}\right)\times\left[g_{*s}^{4}({\rm eV})/g_{*}(T_{\rm
hc})g_{*}^{3}({\rm eV})\right]^{\frac{1}{3}}$ where
$\Omega_{\gamma}=\rho_{\gamma,0}/\rho_{{\rm crit},0}$ is the fraction of the
photon energy density today. Here $g_{*s}({\rm eV})=\left[2+\frac{7}{8}\times
2\times 3\times\left(\frac{4}{11}\right)\right]\simeq 3.91$ and $g_{*}({\rm
eV})=\left[2+\frac{7}{8}\times 2\times
3\times\left(\frac{4}{11}\right)^{4/3}\right]\simeq 3.36$ denote the effective
number of relativistic degrees of freedom relevant for the entropy density and
the energy density, respectively.

As one can see, if horizon crossing occurs during radiation domination $k_{\rm
GW}^{2}=H^{2}(T_{\rm hc})=\rho_{R}/(3M_{P}^{2})$, then the GW spectrum becomes
scale invariant. On the other hand, if horizon crossing occurs during the
inflaton-dominated era, the GW strength is enhanced by a factor of
$\rho_{\phi}/\rho_{R}$ evaluated at $T_{\rm hc}$. As a result, the largest
enhancement is for the mode that re-enters the horizon right after inflation
at $T_{\text{max}}$. For minimal gravitational reheating $(\xi_{h}=0)$, the
enhancement in this mode is
$\rho_{\text{end}}/\rho_{R}(T_{\text{max}})\simeq(4-6)\times 10^{13}$ for
$k\in[6,20]$, based on Eqs. (2.5) and (2.7). This gives $\Omega_{\rm
GW}h^{2}\simeq(8-10)\times 10^{-6}$, which corresponds to the high frequency
points of the blue curves, which fix $\xi_{h}=0$, in the left panel of Fig. 3.
These values are excluded by the BBN bound of $\Omega_{\rm GW}h^{2}\simeq
1.3\times 10^{-6}$ [118], shown by the blue-shaded region at the top.
Therefore, the case with minimal gravitational interactions is excluded, which
was previously pointed out by Ref. [87]. The constraint is relaxed when
$T_{\text{max}}$ is increased, e.g., by non-minimal gravitational interactions
via $\xi_{h}$ [cf. Eq. (3.52)], because the GW energy density relative to that
of radiation is smaller in this case. The blue region in the right panel of
Fig. 3 (and subsequent figures) shows the constraint in this non-minimal
scenario, which excludes $\xi_{h}\lesssim 0.5$, as can be seen from the left
panel of Fig. 2. Alternatively, Ref. [121] offers the solution where the
inflaton energy is more efficiently transferred to radiation via a tachyonic
growth of a new field.

Figure 3: Gravitational wave constraints and future prospects. Left: The blue-
shaded region is excluded by BBN for excessive dark radiation. The blue and
red curves fix $\xi_{h}=0$ and $T_{\text{RH}}=300\,{\rm TeV}$, respectively.
Various curves of the same color use different values of $k$ as labeled and in
increments of 2. The sensitivity of several future experiments as a function
of frequency is also shown. Right: The blue region is excluded by BBN for the
excessive GW energy as dark radiation. The regions below the gray dashed
curves can be probed by the GW experiments as specified.

In addition to setting a constraint, such enhanced gravitational waves offer
an exciting signature to search for [122]. The amount of enhancement depends
on $\rho_{\phi}/\rho_{R}$ at the time of horizon crossing, implying that the
GW spectrum depends on $k$ via $\rho_{\phi}$. By analyzing modes that re-enter
the horizon after $T_{\text{max}}$ and using $\rho_{\phi}\propto
a^{-6k/(k+2)}$ from Eq. (3.5), we find the GW spectrum scales with the
frequency as $\Omega_{\rm GW}h^{2}\propto f^{\frac{k-4}{k-1}}$, which is
consistent with Ref. [87]. The enhanced GW spectra are demonstrated in the
left panel of Fig. 3 for the different values of $k$ in the blue and red
curves. The blue (red) curves correspond to the minimal scenario
($T_{\text{RH}}=300\,{\rm TeV}$), and allow for $k$ to vary from 10 (6) to 20
in increments of 2 (for $\xi_{h}=0$, $k<10$ is excluded by BBN for low
$T_{\text{RH}}$ according to Fig. 2.) Here, the frequency is obtained by
redshifting the initial momentum mode at $T_{\rm hc}$ to today’s photon
temperature $T_{\gamma,0}$ as

$\displaystyle f=\frac{k_{\rm GW}}{2\pi}\frac{T_{\gamma,0}}{T_{\rm
hc}}\left(\frac{g_{*}({\rm eV})}{g_{*}(T_{\rm hc})}\right)^{\frac{1}{3}}.$
(3.54)

Therefore, by measuring the slope of $\Omega_{\rm GW}h^{2}$ as a function $f$,
one can determine $k$ and thus reveal the shape of the inflaton potential
energy near the minimum. Note that, the end-point frequencies for the red
curves are different for different choices of $k$ (and different $\xi_{h}$),
as for a given $k$ and $T_{\text{RH}}$, the maximum possible frequency is
dictated by

$\displaystyle
f_{\text{max}}=\frac{H(T_{\text{max}})}{2\pi}\,\frac{a_{\text{end}}}{a_{0}}\,,$
(3.55)

which for $\xi_{h}=0$ turns out to be $\simeq 7\times 10^{10}$ Hz for all $k$,
while modes with frequencies $f>f_{\text{max}}$ are never produced. In the
right panel of Fig. 3, the regions below the gray dashed curves can be probed
by the future gravitational wave observatories—BBO [123, 124, 125], DECIGO
[126, 127, 128], CE [129, 130] and ET [131, 132, 133, 134]. Here we use the
sensitivity curves derived in Ref. [135]. Since these GW observatories probe
frequencies that correspond to modes that exit the horizon early in the
inflation period, we use the large-field asymptotic value of $V(\phi)$ in Eq.
(2.5) to obtain $H_{I}$. In the left panel of Fig. 3, we illustrate the GW
spectra in the red curves for a fixed $T_{\text{RH}}=300\,{\rm TeV}$, which
can be detected by DECIGO for $k\geq 8$ and by BBO for all $k\geq 6$.
Remarkably, in the right panel, a large region in the parameter space with
$T_{\text{RH}}<5\times 10^{6}\,{\rm GeV}$ can be probed by future GW
detectors. We emphasize that this potential GW signal is generic for our model
and applicable throughout this work, although we do not show these sensitivity
curves in subsequent figures for clarity of presentation.

##  4 Results and discussion

###  4.1 A stable DM candidate

As we have seen from the previous two subsections, for each value of $k$ and
$\xi_{h}$ (in the non-minimal case), there is a unique value for
$T_{\text{RH}}$. These are shown in the left panel of Fig. 2. The
gravitational thermal production of DM generally requires reheating
temperatures much larger than can be obtained with $\xi_{h}=0$. In this
section, we will consider $T_{\text{RH}}$ and $k$ as free parameters and it
should be understood that we are implicitly assuming that $\xi_{h}\neq 0$ and
takes the necessary value to achieve a particular reheating temperature for a
given value of $k$. For the production of DM, both minimal and non-minimal
thermal contributions are included, whereas for the generation of a lepton
asymmetry only minimal contributions from inflaton scattering are considered.

The results presented in this section depend on the underlying class of
inflationary models. As noted earlier, we consider T-models of inflation [88]
for which we have determined $\lambda$ and $\rho_{\text{end}}$. As discussed
above, there are two contributions to the DM relic density: from gravitational
scattering within the newly formed primordial plasma and directly from
inflaton scattering. These two contributions are presented separately in the
upper two panels of Fig. 4. In the upper left panel, we show two contours of
the yield, $n_{N_{1}}/s=10^{-22}$ and $10^{-24}$, for both minimal
gravitational interactions (dotted curves) using Eq. (3.1) and non-minimal
interactions (dot-dashed curves) using Eq. (3.46).999Note that minimal
gravitational interactions ($\xi_{h}=0$) are not actually possible at these
reheating temperatures which require $\xi_{h}\neq 0$.  Note the latter yield
is proportional to $M_{N_{1}}^{2}$ as shown by the contour labels, and we have
normalized these contours by choosing $M_{N_{1}}=10^{8}$ GeV. Also note,
$M_{N_{1}}n_{N_{1}}/s\simeq 0.44$ eV is needed to explain the observed dark
matter density, $\Omega_{N_{1}}h^{2}=0.12$. For $k>4$ and minimal
gravitational interactions, the relic density increases with reheating
temperature, $n_{N_{1}}/s\sim T_{\text{RH}}^{\frac{5k-20}{3k}}$. The scaling
of $n_{N_{1}}/s$ for non-minimal interactions is more complicated but also
increases with $T_{\text{RH}}$.

Figure 4: Top left: Contours of fixed comoving number density,
$n_{N_{1}}/s=10^{-22}$ and $10^{-24}$ in the $(k,T_{\text{RH}})$ plane.
$M_{N_{1}}n_{N_{1}}/s\simeq 0.44$ eV is needed to explain the observed dark
matter density, $\Omega_{N_{1}}h^{2}=0.12$. Dotted curves assume DM production
solely from minimal gravitational scattering in the thermal bath. Dot-dashed
curves correspond to non-minimal gravitational scatterings. The latter are
scaled with $M_{N_{1}}^{2}$. Top right: Contours of
$n_{N_{1}}/s=10^{-20},10^{-18},10^{-16}$ and $10^{-14}$ each scaled with
$M_{N_{1}}^{2}$ assuming DM production only from inflaton scattering. In both
upper panels, the gray-shaded region is excluded as minimal gravitational
interactions necessarily produce larger reheating temperatures. Low reheating
temperatures shaded in red (blue) are excluded by BBN for an excessive
inflaton (GW) energy. Bottom left: The total relic abundance
$(\Omega_{N_{1}}^{T}+\Omega_{N_{1}}^{\phi^{k}})h^{2}/0.12$ as a function of
reheating temperature for three choices of DM masses $\\{10^{6}$, $10^{7}$,
$10^{8}\\}$ GeV for fixed $k=14$. Individual contributions to the dark matter
density are distinguished by line types as indicated. Bottom right: Coloured
regions correspond to values of $(k,T_{\text{RH}})$ with
$(\Omega_{N_{1}}^{T}+\Omega_{N_{1}}^{\phi^{k}})h^{2}\leq 0.12$ for the three
choices of $M_{N_{1}}$ used in the bottom left panel, and the lines styles
indicate the dominant contribution.

In the upper right panel of Fig. 4, we provide four contours of the yield,
$n_{N_{1}}/s$, produced from inflaton scattering, which also scales as
$M_{N_{1}}^{2}$. The gravitational production process from inflaton scattering
is complementary to the thermal production process just discussed. Recall that
we are assuming $\xi_{\phi}$ is small enough that non-minimal scattering
processes can be ignored. In this case, from Eq. (3.19), we see that
$n_{N_{1}}/s\sim T_{\text{RH}}^{-1+\frac{4}{k}}$ and for $k>4$, the relic
density decreases with increasing $T_{\text{RH}}$. Recalling that
$M_{N_{1}}n_{N_{1}}/s\simeq 0.44$ eV is needed to explain the observed DM
density, $\Omega_{N_{1}}h^{2}=0.12$, we find indeed that higher reheating
temperatures require lighter DM candidates to fit with the relic abundance
constraint.

Combining the two constraints shown in the top panels of Fig. 4 we see that
for a given $k$ and $M_{N_{1}}$, there are both upper (from thermal
scattering) and lower (from inflaton scattering) limits to $T_{\text{RH}}$ so
as to avoid exceeding the observed cold DM abundance. The resulting relic
density as a function of $T_{\text{RH}}$ is shown in the bottom left panel of
Fig. 4, where we show the total relic abundance
$(\Omega_{N_{1}}^{T}+\Omega_{N_{1}}^{\phi^{k}})\,h^{2}$ relative to the
observed abundance for a fixed $k=14$ and three choices of the DM mass
$M_{N_{1}}=\\{10^{6}$, $10^{7}$, $10^{8}\\}$ GeV. We clearly see that the
desired relic density $(\Omega_{N_{1}}=0.12)$ is obtained twice: (i) at a
lower reheating temperature, where inflaton scattering dominates, and (ii) for
a higher reheating temperature, when we are in the thermal production regime.
The allowed region corresponds to the parameter space at or below the line
$\Omega h^{2}/0.12=1$ in the bottom left panel of Fig. 4. For
$M_{N_{1}}>3\times 10^{8}$ GeV, there are no values of ($T_{\text{RH}}$, $k$)
that result in an acceptable density of DM, and the allowed range in
$T_{\text{RH}}$ is larger with lighter DM. This is understandable as, the
thermal relic requiring a larger upper bound on $T_{\text{RH}}$ for lighter
DM, while the inflaton scattering requires a smaller lower bound on
$T_{\text{RH}}$ for lighter DM.

A two-dimensional version of the lower left panel of Fig. 4, over a range in
$k$, is shown in the lower right panel. Low values of $T_{\text{RH}}$ are
excluded by BBN. Once again, the gray-shaded region in the lower right corner
of this panel is also excluded since minimal gravitational interactions would
produce a reheating temperature larger than the values in that region. Within
each shaded band (the color corresponds to a specific choice of $M_{N_{1}}$),
the total relic density is below the observed DM density. The observed value
is reached on the border of the colored bands. For DM of masses very close to
$1$ PeV, there exists a viable parameter space for $k\geq 9$ (along the
boundary of the excessive GWs region), requiring $\xi_{h}\simeq 0.5$. For
larger masses, the range in $k$ extends to lower values, and higher reheating
temperatures are possible and require larger non-minimal coupling to gravity.

Having identified the regions of the $(k,T_{\text{RH}})$ parameter space with
a suitable DM density, we turn to the production of the baryon asymmetry
through gravitationally induced leptogenesis. This analysis was performed in
[84] and therefore we only briefly summarize the results found there. We note,
however, Ref. [84] neglected the kinematic suppression in Eq. (3.17) to
maintain the model independence of the analysis, though this effect is
included in the present work. In Fig. 5, we show contours of some benchmark
values of the mass of $N_{2}$ that reproduce the observed baryon asymmetry
$Y_{B}^{\text{obs}}$. We find that the gravitational contribution to the
baryon asymmetry is essentially entirely due to inflaton scattering rather
than the thermal particles in the SM bath. Since minimal gravitational
interactions are excluded by excessive GWs, non-minimal interactions are
required to produce a sufficiently large thermal bath so that GW fractional
energy is consistent with BBN. Leptogenesis via $N_{2}$ is therefore possible
above the border of the blue-shaded region in Fig. 5, indicating a mass
$M_{N_{2}}\gtrsim 3\times 10^{11}$ GeV is required. Larger values of
$M_{N_{2}}$ can produce the correct asymmetry so long as $\xi_{h}>0$.
Nonetheless, when $M_{N_{2}}\gtrsim 3\times 10^{12}\,{\rm GeV}$, the baryon
asymmetry starts to become suppressed for the following reason. The inflaton
mass obtained from Eqs. (2.5) and (2.7), $m_{\phi}\simeq 1.2\times
10^{13}\,{\rm GeV}$ across all $k$ values, is no longer much larger than
$M_{N_{2}}$ and the kinematic suppression in Eq. (3.17) becomes important.
This explains the existence of the green region as well as why the curve for
$M_{N_{2}}=10^{13}\,{\rm GeV}$ is at a lower $T_{\text{RH}}$ than that for
$M_{N_{2}}=3\times 10^{12}\,{\rm GeV}$. Once again, the bottom red region is
forbidden by BBN because of an excessive inflaton energy density during BBN.
In summary, we observe that, saturating the bound on GWs from BBN, together
with the right DM abundance and successful leptogenesis requires
$\xi_{h}\simeq 0.5$, $M_{N_{2}}\simeq 3\times 10^{11}$ GeV and
$M_{N_{1}}\simeq 10^{6}$ GeV. As discussed above, this parameter space can be
extended, allowing larger values $\\{M_{N_{1}},M_{N_{2}}\\}$ if one considers
stronger non-minimal gravitational couplings by $\xi_{h}\gtrsim 0.5$, thus
allowing a larger reheating temperature [cf. Eq. (3.51)].

Figure 5: Contours of $M_{N_{2}}$ corresponding to the observed baryon
asymmetry [cf. Eq. (3.24)] in the $(k,T_{\text{RH}})$ plane. The red-shaded
region correspond to the lower bound on $T_{\text{RH}}$ from BBN, and the
green region leads to underproduction of $Y_{B}$ due to the kinematic
suppression in inflaton scattering when $M_{N_{2}}$ approaches $m_{\phi}$.

Combining our preceding analyses, it is possible, for a given $V(\phi)$, to
constrain the ($M_{N_{1}}$, $M_{N_{2}}$, $\xi_{h}$) parameter space by
requiring leptogenesis, DM production and reheating to have a common
gravitational origin. Indeed, for a given $k$ and DM mass $M_{N_{1}}$, the
temperature $T_{\text{RH}}$ can be determined by the relic abundance
constraint. In turn, $T_{\text{RH}}$ determines the value of $\xi_{h}$ needed
to reheat the Universe, as well as the value of $M_{N_{2}}$ which gives the
desired baryon asymmetry through leptogenesis. To illustrate our result, we
project the viable parameter space in the ($M_{N_{1}}$,$M_{N_{2}}$) plane in
Fig. 6 for different values of $\xi_{h}$, allowing $k$ to vary within
$k\in[6,\,20]$. In each coloured line segment, gravitational interactions are
responsible for the observed DM relic abundance, the baryon asymmetry and
reheating. Different coloured slanted line segments in this figure correspond
to different choices of the non-minimal coupling $\xi_{h}$, with $\xi_{h}=0$
being ruled out from overproduction of GWs. The maximum possible value for
$\xi_{h}$ is around 13.5, above which the mass $M_{N_{2}}$ necessary to
reproduce the observed baryon asymmetry gets too close to $m_{\phi}$ and
kinematic suppression becomes significant, as can be seen from Fig. 5. Note
that for each $\xi_{h}$, the allowed parameter space satisfying all the
constraints, is rather restricted. This is better seen from the right panel
figure, where we have zoomed in to the $\xi_{h}=1$ case. Interestingly, this
shows that the viable parameter space is approximately independent of $k$,
while $k=6$ and $8$ are excluded by BBN as can be seen from the left panel of
Fig. 2.

Figure 6: Viable parameter space in the ($M_{N_{1}},M_{N_{2}}$) plane for
which gravitational interactions are responsible for the observed DM relic
abundance (in $N_{1}$), the baryon asymmetry (produced from $N_{2}$ decays),
and reheating for $k\in[6,\,20]$. In the left panel, different colours
correspond to $\xi_{h}=\\{0,\,1,\,10\\}$ diagonally from bottom left (red) to
top right (blue). The vertical black segment indicates the range in
$M_{N_{2}}$ for $M_{N_{1}}=4$ PeV for the range in $k$ considered, where the
connection to the IceCube high-energy neutrino excess will be discussed in the
next subsection. In the right panel, we magnify the parameter space for a
fixed non-minimal coupling $\xi_{h}=1$. The dots correspond to even values of
$k$ as indicated.

###  4.2 The case for a decaying gravitational DM & IceCube events

Until now, we have assumed that the DM candidate, $N_{1}$, is absolutely
stable. If it is not, and $N_{1}$ has non-zero Yukawa components, $y_{1i}$,
$N_{1}$ can decay to SM final states. In this case, one necessary (but not
sufficient) constraint on the DM mass and Yukawa coupling arises from the
requirement of having a lifetime larger than the age of the Universe
$\tau_{N_{1}}\gtrsim\tau_{\text{univ}}\simeq 4.35\times 10^{17}$ s. On the
other hand, the IceCube detector has reported the detection of three PeV
neutrinos, a roughly $3\sigma$ excess above expected background rates [136,
137, 138]. The three highest energy events correspond to deposited energies of
1.04 PeV, 1.14 PeV and 2.0 PeV. Although the origin of these very high energy
events is still unclear, it has been shown in Refs. [139, 140, 141, 142, 143,
144, 145, 146, 147, 148, 149, 150, 151, 152, 42] that such events could be
sourced from the decays of superheavy DM particles. The neutrino energy
spectrum presents a high-energy cutoff at $m_{\text{DM}}/2$ [140, 141] if two
body decays including one neutrino are present. The total excess can be
interpreted as high energy neutrinos resulting from the decay of $N_{1}$ with
$\tau_{N_{1}}\approx 10^{28}$ s for both normal and inverted hierarchies [140,
153]. Given that the maximum energy of the IceCube events has been measured to
be about 2 PeV, the mass of the DM particle is constrained to be $\simeq$ 4
PeV. Moreover, the IceCube spectrum sets a lower bound on the DM lifetime of
the order of $10^{28}$ s [141, 153], which is approximately model-independent
and orders of magnitude larger than the lifetime of the Universe. Thus,
satisfying this bound automatically makes $N_{1}$ a nearly stable relic, and
hence a good DM candidate. For $N_{1}\to\ell\,H$ decay, we find

$\displaystyle\tau\equiv\Gamma_{N_{1}\to\ell\,H}^{-1}\simeq\left(\frac{y_{N_{1}}^{2}M_{N_{1}}}{8\,\pi}\right)^{-1}\simeq
10^{28}~{}\text{s}\,\left(\frac{4\times
10^{-29}}{y_{N_{1}}}\right)^{2}\,\left(\frac{1\,\text{PeV}}{M_{N_{1}}}\right)\,;$
(4.1)

that is, the Yukawa coupling $y_{N_{1}}$ must be highly suppressed.

On the other hand, in order to satisfy the observed DM abundance via the
freeze-in mechanism in the early Universe through inverse decay: $\nu\,,H\to
N_{1}$ involving the same Yukawa, one needs [32, 150]

$\displaystyle\Omega_{N_{1}}\,h^{2}\simeq
0.12\,\left(\frac{y_{N_{1}}}{1.2\times
10^{-12}}\right)^{2}\,\left(\frac{M_{N_{1}}}{1\,\text{PeV}}\right)\,.$ (4.2)

This means that the Yukawa required to interpret the PeV IceCube event from
the decay of $N_{1}$, $y_{N_{1}}\sim 10^{-29}$, is far too small for the
thermal bath to populate the Universe from the freeze in mechanism. Thus, if
we are restricted to dimension-four interactions involving RHN and the SM, it
is not possible to simultaneously explain the DM relic density and the IceCube
events. Alternatively, we may consider higher dimensional operators [150],
modified gravity/cosmology [154] or some different production mechanism for DM
[155].

Our minimalistic framework contains a natural avenue to reconcile both the DM
abundance and IceCube events, through the gravitational production of decaying
PeV neutrinos in the early Universe. However, as discussed in Sec. 3.5, the
case with minimal gravitational interactions is excluded by BBN for an
excessive amount of gravitational waves as dark radiation. Thus we then need
to go (slightly) beyond the minimal setup and include non-minimal
gravitational interactions. We show in Fig. 7 contours for
$\Omega_{N_{1}}h^{2}=0.12$ for $M_{N_{1}}=4$ PeV in the ($k$, $T_{\text{RH}}$)
plane. The orange (dashed, dotted) lines correspond to the two dominant
gravitational scattering processes involving the (inflaton, thermal particles)
as discussed in the previous subsection. Note however that gravitational
thermal production requires a high reheating temperature and is not compatible
with the observed baryon asymmetry as can be understood from Fig. 5. In
contrast, at lower $T_{\text{RH}}$, the correct relic density can be produced
from inflaton scattering with a lower value of $\xi_{h}\approx 2.5$. In the
left panel of Fig. 6, we show, by the black vertical line segment, the range
in $M_{N_{2}}$ obtained from varying $k$ while fixing $M_{N_{1}}=4$ PeV. Note
that, since $N_{1}$ is a long-lived stable relic, it does not contribute to
the generation of the baryon asymmetry as its decay takes place below the
electroweak phase transition. In addition, because the Yukawa coupling of
$N_{1}$ is extremely small, its interactions which violate lepton number will
not be in equilibrium, and hence will not wash out any of the asymmetry
produced by $N_{2}$. We summarize our analysis in the tables presented in the
conclusion.

Figure 7: Contours of fixed relic density, $\Omega_{N_{1}}h^{2}=0.12$ for
$M_{N_{1}}=4$ PeV. The upper dotted contour corresponds to production from
gravitational scattering in the thermal bath (and requires a large value of
$\xi_{h}$) and the lower dashed contour corresponds to production from
inflaton scattering (and requires a relatively low value of $\xi_{h}$) Between
the two contours $\Omega_{N_{1}}h^{2}<0.12$ for $M_{N_{1}}=4$ PeV.

##  5 Dark matter & leptogenesis with a Majoron

We have seen in the previous section that our result was particularly
constrained because of the strong dependence of the production of the RHN on
its mass $M_{N_{i}}$, limiting our allowed region to masses above a PeV. In
this section, we consider an alternative mechanism. We include an additional
complex scalar field, $\Phi$ containing the Majoron [156, 99, 157, 158, 159,
160, 145, 161, 162, 163], that acts as an intermediate state in the
interactions of the inflaton and RHNs. This interaction is depicted in Fig. 8.

Figure 8: Feynman diagram for the gravitational production of an on-shell
scalar $S$ coupled to the heavy neutrinos.

The relevant Lagrangian of this extension can be written as

$\displaystyle{\cal
L}_{\Phi}=(-y_{R}^{i}\,\Phi\,\overline{N_{i}^{c}}\,N_{i}+\text{h.c.})+\frac{1}{2}\mu_{\Phi}^{2}\,\Phi^{2}-\frac{1}{4}\lambda_{\Phi}\,\Phi^{4}\,.$
(5.1)

After symmetry breaking, the real part of $\Phi$ acquires a non-zero vacuum
expectation value, around which one can expand the field as:
$\Phi=\frac{1}{\sqrt{2}}(S+v_{S})e^{iJ/v_{S}}$, and $J$ is the Majoron. This
expectation value is the origin of the RHN Majorana masses,
$M_{N_{i}}=y_{R}^{i}\,v_{S}/\sqrt{2}$. Then $m_{S}=\mu_{\Phi}<m_{\phi}$ and
the gravitational production rate of the real scalar, $S$ is101010As shown in
Refs. [68, 67], for the case of a scalar field, the gravitational thermal
production is always negligible with respect to inflaton scattering.

$R^{\phi^{k}}_{S}=\frac{2\times\rho_{\phi}^{2}}{16\pi
M_{P}^{4}}\Sigma_{S}^{k}\,,$ (5.2)

where the factor of two accounts for the fact we produce two scalar particles
per scattering, with [65, 68]

$\Sigma_{S}^{k}=\sum_{n=1}^{\infty}|{\cal
P}^{k}_{2n}|^{2}\left[1+\frac{2\mu_{\Phi}^{2}}{E_{2n}^{2}}\right]^{2}\sqrt{1-\frac{4\mu_{\Phi}^{2}}{E_{2n}^{2}}}\,.$
(5.3)

Since each scalar decays into 2 right-handed neutrinos, we obtain for the
density of $N_{i}$ after integration of the Boltzmann equation [68],

$\displaystyle n_{N_{i}}^{S\phi^{k}}(a_{\rm
RH})\simeq\text{Br}_{i}\times\frac{\sqrt{3}\rho_{\rm RH}^{3/2}}{4\pi
M_{P}^{3}}\frac{k+2}{6k-6}\left(\frac{\rho_{\rm end}}{\rho_{\rm
RH}}\right)^{1-\frac{1}{k}}\Sigma^{k}_{S},$ (5.4)

where we assumed $a_{\rm RH}\gg a_{\rm end}$, and here
$\text{Br}_{i}=\frac{(y_{R}^{i})^{2}}{\sum(y_{R}^{i})^{2}}$ so
$\text{Br}_{i}=\frac{M_{N_{i}}^{2}}{M_{N_{1}}^{2}+M_{N_{2}}^{2}+M_{N_{3}}^{2}}$
if $N_{1,2,3}$ are all lighter than $S$. The relic abundance of $N_{1}$ is
then given by

$\displaystyle\frac{\Omega_{N_{1}}^{S\phi^{k}}h^{2}}{0.12}$
$\displaystyle\simeq$ $\displaystyle\text{Br}_{1}\times\left(\frac{\rho_{\rm
end}}{10^{64}{\rm GeV}^{4}}\right)^{1-\frac{1}{k}}\left(\frac{10^{40}{\rm
GeV}^{4}}{\rho_{\rm
RH}}\right)^{\frac{1}{4}-\frac{1}{k}}\left(\frac{k+2}{6k-6}\right)$ (5.5)
$\displaystyle\times$
$\displaystyle\Sigma_{S}^{k}\times\frac{M_{N_{1}}}{2.5\times
10^{\frac{24}{k}-8}{\rm GeV}}\,,$

whereas the baryon asymmetry follows from Eq. (3.24). Note that so long as
$M_{N_{i}}\ll\mu_{\Phi}\ll m_{\phi}$, the resulting dark matter abundance and
baryon asymmetry will be independent of $m_{S}$.

We show in Figs. 9 and 10 respectively, the parameter space allowed by the
relic abundance and the baryogenesis constraint in the ($k\,,T_{\text{RH}}$)
plane. Comparing Fig. 9 and the dashed lines (from the inflaton scattering) in
the bottom right panel of Fig. 4, we notice that the mass of the dark matter
respecting Planck constraint is much lower, if the branching fraction to
$N_{1}$ is large. For $\text{Br}_{1}=1$, the difference is about 8 orders of
magnitude, and around 6 orders of magnitude for $\text{Br}_{1}=10^{-2}$. The
reason is easy to understand: the production rate through $S$ is boosted in
comparison with the direct production, by a factor

$\displaystyle\frac{R^{S\phi^{k}}_{N_{1}}}{R^{\phi^{k}}_{N_{1}}}\simeq\text{Br}_{1}\frac{m_{\phi}^{2}}{M_{N_{1}}^{2}}\,.$
(5.6)

For smaller branching fraction, the density of $N_{1}$ through this channel is
suppressed and the effect is milder and proportional to $\text{Br}_{1}$, as
one can see in Fig. 9 right panel.

Figure 9: Contours of observed relic abundance assuming $\text{Br}_{1}=1$
(left) and $\text{Br}_{1}=10^{-2}$ (right) for different choices of the DM
mass, considering only Majoron contribution. The purple-shaded region is
disallowed from the warm DM limit (see text).

Figure 10: Contours of $N_{N_{2}}$ corresponding to the observed baryon
asymmetry for $\text{Br}_{2}=1$ (left) and $\text{Br}_{2}=10^{-2}$ (right) in
the $(k,T_{\text{RH}})$ plane. Here only the contribution due to the
intermediate scalar is included. The green-shaded region is kinematically
inaccessible due to $M_{N_{2}}>m_{\phi}/2$ [cf. Fig. 8].

Given that the required mass, $M_{N_{1}}$, can be much lower when we couple
the RHNs to the scalar $S$, and $N_{1}$ is produced relativistically, $N_{1}$
dark matter may still be warm around the time of CMB decoupling. We derive the
warmness constraint by redshifting the $N_{1}$ initial momentum of order
$m_{\phi}$ at $T_{\text{max}}$ to the temperature $T\simeq 1$ eV and require
that the velocity at $T\simeq 1$ eV is less than $2\times 10^{-4}$. This bound
on the velocity is obtained from translating the limit on the warm dark matter
mass from the Lyman-$\alpha$ forest [164] in the case where the abundance is
generated thermally. The current warmness constraint is shown by the purple
region, while the future sensitivity using cosmic 21-cm lines [165] is shown
by the purple dotted curve. In summary, this mechanism interestingly allows
for electroweak scale fermionic dark matter produced gravitationally, which is
not possible by the direct scattering of the inflaton. We show in Fig. 10 the
parameter space allowed to obtain a sufficient amount of baryon asymmetry for
the set of branching ratios $\text{Br}_{2}=1$ and $10^{-2}$. Comparing Figs. 5
and 10 left, we note that for $M_{2}=10^{13}$ GeV the situation is similar to
the direct production because no real enhancement
$\propto\frac{m_{\phi}^{2}}{M_{N_{2}}^{2}}$ exists. However, for
$M_{N_{2}}=10^{11}$ GeV and large values of $k$, $T_{\text{RH}}$ should be
about 3 orders of magnitude larger to obtain the same asymmetry. The reason is
that for a large value of $k$, $Y_{B}\propto\frac{1}{T_{\text{RH}}}$ when $S$
is produced (combining Eqs. (5.4) and (3.24)), and
$\propto\frac{M_{N_{2}}^{2}}{m_{\phi}^{2}T_{\text{RH}}}$ when $N_{2}$ is
produced directly. In other words, $T_{\text{RH}}$ should be compensated by a
factor $\frac{m_{\phi}^{2}}{M_{N_{2}}^{2}}$ to avoid an excessive asymmetry.
As in the case of dark matter, lowering the branching ratio dilutes the effect
as one can see in the right panel of Fig. 10.

Finally, we can combine all of the preceding results, adding the possibility
for a gravitational reheating with non-minimal coupling. We illustrate this in
Fig. 11, which is the analogue of Fig. 6 but with the scalar $S$ as an
intermediate state. For demonstration purposes, here we suppose
$M_{N_{3}}>m_{\phi}/2$ so that $N_{3}$ is not produced by the inflaton or $S$,
resulting in $\text{Br}_{2}=1-\text{Br}_{1}=1-(M_{N_{1}}/M_{N_{2}})^{2}\simeq
1$. As the branching ratios are completely determined by the masses
$M_{N_{1}}$ and $M_{N_{2}}$, for a fixed $k$ and a fixed $\xi_{h}$, there will
be again only one point in the $(M_{N_{1}},M_{N_{2}})$ plane that could
simultaneously obtain the CMB-determined DM relic abundance and the observed
baryon asymmetry. Each color segment in Fig. 11 assumes a fixed $\xi_{h}$ and
allows values of $k\in[6,20]$ that are consistent with the BBN bound on
$T_{\text{RH}}$. The black dot indicates the $M_{N_{1,2}}$ masses, independent
of $k$, required to explain the IceCube high-energy neutrino excess. The green
region is inaccessible because $m_{S}>M_{N_{2}}>m_{\phi}/2$ forbids the
production of $S$ via $\phi$ scattering.111111Note that we have not included
the effects of $\mu_{\Phi}$ in $\Sigma^{k}_{S}$. These start to play a role
for large $\xi_{h}\gtrsim 10^{2}$ when $2M_{N_{2}}$ becomes close to
$m_{\phi}$, since we require $2m_{N_{2}}<m_{S}<m_{\phi}$.  Compared to Fig. 6,
we see that the effect of $S$ as an intermediate state expands the parameter
space of allowed dark matter density and baryon asymmetry. Most notably, the
parameter space opens up towards lower masses, and allows large values of
$\xi_{h}$.

Figure 11: Parameter space satisfying the right dark matter relic abundance
and baryon asymmetry, considering the production through $S$. The line colors
correspond to different values of $\xi_{h}$, with
$\xi_{h}=\\{0,1,10,10^{2}\\}$ from bottom to top, and $\xi_{h}=0$ corresponds
to minimal graviton exchange. Each colored line segment shows the variation of
the predicted masses with $k\in\left[6,\,20\right]$. The black dot marks the
parameter point that can also explain the IceCube high-energy neutrino excess.

##  6 Conclusions

In this paper, we have shown that there exists the possibility that
inflationary reheating, dark matter and the baryon asymmetry can be generated
solely gravitational interactions. The baryon asymmetry is produced through
the decay of a right-handed neutrino $M_{N_{2}}$ (leading first to a non-zero
lepton asymmetry). For minimal gravitational interactions, $\xi_{h}=0$, a
large amount of dark radiation is created in the form of gravitational waves
and is inconsistent with BBN. Thus, we allow for a non-minimal gravitational
coupling $\xi_{h}RH^{2}$ where $H$ the Standard Model Higgs field to enhance
reheating, so that the ratio of gravitational wave energy density to the
radiation is decreased. The lowest $\xi_{h}$ consistent with BBN is around
0.5. The range of the parameter space is 2-8 PeV for the dark matter mass
$M_{N_{1}}$ and $0.3-1.7\times 10^{12}$ GeV for the mass of the lepton number
violating decaying RHN, $M_{N_{2}}$. The range corresponds to a scan over
$\xi_{h}$. Our solution restricts $0.5\lesssim\xi_{h}\lesssim 13.5$ and
$T_{\text{RH}}<3\times 10^{5}$ GeV, where the maximum reheating temperature is
attained with $\xi_{h}\simeq 4.7$. We summarize our results for fixing
different values of $\xi_{h}$ or $M_{N_{1}}$ in the upper sections of Table 2
and Table 3 labeled “Direct gravitational production.” Primordial
gravitational waves generated during inflation allow a large parameter space
with $T_{\text{RH}}\lesssim 5\times 10^{6}$ GeV to be probed in proposed
gravitational wave detectors such as BBO, DECIGO, CE and ET.

Direct gravitational production
---
$\xi_{h}$ | $T_{\text{RH}}\,\text{[GeV]}$ | $M_{N_{1}}\,\text{[PeV]}$ | $M_{N_{2}}\,\text{[GeV]}$
1 | $\\{5.6-1.8\times 10^{4}\\}$ | $\\{2.49-2.54\\}$ | $\\{4.5-4.8\\}\times 10^{11}$
$\\{2.5-2.7\\}$ | $\\{0.11-9.8\times 10^{4}\\}$ | $4.0\star$ | $\\{7.3-9.2\\}\times 10^{11}$
10 | $8.1$ | $8.1$ | $1.7\times 10^{12}$
Gravitational production via $S$
$\xi_{h}$ | $T_{\text{RH}}\,\text{[GeV]}$ | $M_{N_{1}}$ | $M_{N_{2}}\,\text{[GeV]}$
1 | $\\{5.6-1.8\times 10^{4}\\}$ | $\\{7.9-12\\}~{}{\rm TeV}$ | $\\{1.4-2.1\\}\times 10^{9}$
$10$ | $\\{8.1-1.3\times 10^{6}\\}$ | $\\{220-360\\}~{}{\rm TeV}$ | $\\{4.0-6.5\\}\times 10^{10}$
$\\{50-68\\}$ | $\\{2.6\times 10^{3}-2.6\times 10^{7}\\}$ | $4.0~{}{\rm PeV}\star$ | $7.1\times 10^{11}$
$100$ | $\\{8.1\times 10^{3}-9.8\times 10^{7}\\}$ | $\\{7.1-11\\}~{}{\rm PeV}$ | $\\{1.3-2.0\\}\times 10^{12}$
Table 2: Ranges of $T_{\text{RH}}$, DM $(M_{N_{1}})$ and RHN $(M_{N_{2}})$
masses, over which baryon asymmetry and DM relic abundance are simultaneously
satisfied via gravitational yield for different choices of the non-minimal
coupling $\xi_{h}$; $\xi_{h}=0$ corresponds to minimal gravity, which is not
shown since it is excluded by BBN for excessive GWs. The upper section assumes
direction production of $N_{1,2}$, while the lower assumes production via
Majoron’s CP-even partner $S$ as an intermediate state. Here we allow
$k\in[6,20]$ and omit points whose low $T_{\text{RH}}$ values are excluded by
BBN. (For the direct gravitational production, a single value of $k=8$ is
allowed for $\xi_{h}=10$, while no points are allowed for $\xi_{h}\gtrsim
13.5$.) The $\star$ entry corresponds to the DM mass that can explain the
IceCube high-energy neutrino events. Direct gravitational production
---
$\xi_{h}$ | $T_{\text{RH}}\,\text{[GeV]}$ | $M_{N_{1}}\,\text{[PeV]}$ | $M_{N_{2}}\,\text{[GeV]}$
1 |  $0.0084$ (excluded) | – | –
2.5 | 0.11 | $4.0\star$ | $7.3\times 10^{11}$
10 | $8.1$ | $8.1$ | $1.7\times 10^{12}$
Gravitational production via $S$
$\xi_{h}$ | $T_{\text{RH}}\,\text{[GeV]}$ | $M_{N_{1}}$ | $M_{N_{2}}\,\text{[GeV]}$
1 | 0.0084 (excluded) | – | –
10 | $8.1$ | $220~{}{\rm TeV}$ | $4.0\times 10^{10}$
68 | $2.6\times 10^{3}$ | $4.0\star$ | $7.1\times 10^{11}$
100 | $8.1\times 10^{3}$ | $7.1~{}{\rm PeV}$ | $1.3\times 10^{12}$
Table 3: Same as Table 2 but for a fixed $k=8$ as a benchmark.

We also showed that $N_{1}$, if unstable, can explain the recent IceCube PeV
events through its decay into SM neutrinos. In this case, if we want to
accommodate simultaneously the correct DM relic abundance, the observed baryon
asymmetry, gravitational reheating and the IceCube events, the value of
$\xi_{h}$ is fixed for a given $k$. We show this result in the second row of
the upper section of each table where the assumed value of $m_{N_{1}}=4$ PeV
is marked by $\star$.

Finally, we proposed a new scenario where the RHN and the dark matter are
produced through an intermediate scalar state $S$, the CP-even partner of the
Majoron. In this case, the gravitational production of the scalar is not
helicity suppressed by the mass of the final state fermions. As a result, the
mass ranges for $N_{1}$ and $N_{2}$ are increased. For
$0.5\lesssim\xi_{h}\lesssim 100$, the mass range for $N_{1}$ is $4$ TeV to
$11$ PeV, and the range for $N_{2}$ is $7\times 10^{8}-2\times 10^{12}$ GeV.
Finally, the IceCube events can be explained by appropriate choices of
$\xi_{h}$ and $T_{\text{RH}}$ for each $k$, whereas $M_{N_{2}}$ is predicted
to be $7.1\times 10^{11}$ GeV for all $k$ when $M_{N_{1}}$ is fixed to 4 PeV.
We summarize our overall results for fixing different values of $\xi_{h}$ or
$M_{N_{1}}$ in the lower sections of Table 2 and Table 3 labeled
“Gravitational production via $S$”.

## Acknowledgements

We would like to thank Essodjolo Kpatcha and Jong-Hyun Yoon for useful
discussions. We thank Pedro Schwaller for bringing to our attention the
constraint from the gravitational wave as dark radiation. BB received funding
from the Patrimonio Autónomo - Fondo Nacional de Financiamiento para la
Ciencia, la Tecnología y la Innovación Francisco José de Caldas (MinCiencias -
Colombia) grant 80740-465-2020. This project has received funding /support
from the European Union’s Horizon 2020 research and innovation programme under
the Marie Sklodowska-Curie grant agreement No 860881-HIDDeN, and the IN2P3
Master Projet UCMN. BB would also like to acknowledge the hospitality at the
IJCLab during INVISIBLES’22 school+workshop, during which this work was
initiated and fruitful discussions with Ashmita Das. The work of R.C. and
K.A.O. was supported in part by DOE grant DE-SC0011842 at the University of
Minnesota.

## References

  * [1] H. Poincaré, _The Milky Way and the Theory of Gases_ , _Popular Astronomy, vol. 14, pp.475-488_ 14 (1906) 475.
  * [2] F. Zwicky, _Die Rotverschiebung von extragalaktischen Nebeln_ , _Helv. Phys. Acta_ 6 (1933) 110.
  * [3] Planck collaboration, _Planck 2018 results. VI. Cosmological parameters_ , 1807.06209.
  * [4] J.E. Gunn, B.W. Lee, I. Lerche, D.N. Schramm and G. Steigman, _Some Astrophysical Consequences of the Existence of a Heavy Stable Neutral Lepton_ , _Astrophys. J._ 223 (1978) 1015.
  * [5] P. Hut, _Limits on Masses and Number of Neutral Weakly Interacting Particles_ , _Phys. Lett. B_ 69 (1977) 85.
  * [6] B.W. Lee and S. Weinberg, _Cosmological Lower Bound on Heavy Neutrino Masses_ , _Phys. Rev. Lett._ 39 (1977) 165.
  * [7] J.R. Ellis, J.S. Hagelin, D.V. Nanopoulos, K.A. Olive and M. Srednicki, _Supersymmetric Relics from the Big Bang_ , _Nucl. Phys. B_ 238 (1984) 453.
  * [8] V. Silveira and A. Zee, _SCALAR PHANTOMS_ , _Phys. Lett. B_ 161 (1985) 136.
  * [9] J. McDonald, _Gauge singlet scalars as cold dark matter_ , _Phys. Rev._ D50 (1994) 3637 [hep-ph/0702143].
  * [10] C.P. Burgess, M. Pospelov and T. ter Veldhuis, _The Minimal model of nonbaryonic dark matter: A Singlet scalar_ , _Nucl. Phys. B_ 619 (2001) 709 [hep-ph/0011335].
  * [11] H. Davoudiasl, R. Kitano, T. Li and H. Murayama, _The New minimal standard model_ , _Phys. Lett. B_ 609 (2005) 117 [hep-ph/0405097].
  * [12] J.M. Cline, K. Kainulainen, P. Scott and C. Weniger, _Update on scalar singlet dark matter_ , _Phys. Rev. D_ 88 (2013) 055025 [1306.4710].
  * [13] H. Han and S. Zheng, _New Constraints on Higgs-portal Scalar Dark Matter_ , _JHEP_ 12 (2015) 044 [1509.01765].
  * [14] A. Djouadi, O. Lebedev, Y. Mambrini and J. Quevillon, _Implications of LHC searches for Higgs–portal dark matter_ , _Phys. Lett. B_ 709 (2012) 65 [1112.3299].
  * [15] O. Lebedev, H.M. Lee and Y. Mambrini, _Vector Higgs-portal dark matter and the invisible Higgs_ , _Phys. Lett. B_ 707 (2012) 570 [1111.4482].
  * [16] Y. Mambrini, _Higgs searches and singlet scalar dark matter: Combined constraints from XENON 100 and the LHC_ , _Phys. Rev. D_ 84 (2011) 115017 [1108.0671].
  * [17] A. Djouadi, A. Falkowski, Y. Mambrini and J. Quevillon, _Direct Detection of Higgs-Portal Dark Matter at the LHC_ , _Eur. Phys. J. C_ 73 (2013) 2455 [1205.3169].
  * [18] J.A. Casas, D.G. Cerdeño, J.M. Moreno and J. Quilis, _Reopening the Higgs portal for single scalar dark matter_ , _JHEP_ 05 (2017) 036 [1701.08134].
  * [19] G. Arcadi, Y. Mambrini and F. Richard, _Z-portal dark matter_ , _JCAP_ 03 (2015) 018 [1411.2985].
  * [20] M. Escudero, A. Berlin, D. Hooper and M.-X. Lin, _Toward (Finally!) Ruling Out Z and Higgs Mediated Dark Matter Models_ , _JCAP_ 12 (2016) 029 [1609.09079].
  * [21] L. Roszkowski, E.M. Sessolo and S. Trojanowski, _WIMP dark matter candidates and searches—current status and future prospects_ , _Rept. Prog. Phys._ 81 (2018) 066201 [1707.06277].
  * [22] G. Arcadi, M. Dutra, P. Ghosh, M. Lindner, Y. Mambrini, M. Pierre et al., _The waning of the WIMP? A review of models, searches, and constraints_ , _Eur. Phys. J. C_ 78 (2018) 203 [1703.07364].
  * [23] Y. Mambrini, _The Kinetic dark-mixing in the light of CoGENT and XENON100_ , _JCAP_ 09 (2010) 022 [1006.3318].
  * [24] O. Lebedev and Y. Mambrini, _Axial dark matter: The case for an invisible Z’_ , _Phys. Lett. B_ 734 (2014) 350 [1403.4837].
  * [25] A. Alves, A. Berlin, S. Profumo and F.S. Queiroz, _Dark Matter Complementarity and the Z ′ Portal_, _Phys. Rev. D_ 92 (2015) 083004 [1501.03490].
  * [26] A. Alves, G. Arcadi, Y. Mambrini, S. Profumo and F.S. Queiroz, _Augury of darkness: the low-mass dark Z’ portal_ , _JHEP_ 04 (2017) 164 [1612.07282].
  * [27] G. Arcadi, P. Ghosh, Y. Mambrini, M. Pierre and F.S. Queiroz, _$Z^{\prime}$ portal to Chern-Simons Dark Matter_, _JCAP_ 1711 (2017) 020 [1706.04198].
  * [28] D.V. Nanopoulos, K.A. Olive and M. Srednicki, _After Primordial Inflation_ , _Phys. Lett. B_ 127 (1983) 30.
  * [29] M.Y. Khlopov and A.D. Linde, _Is It Easy to Save the Gravitino?_ , _Phys. Lett. B_ 138 (1984) 265.
  * [30] K.A. Olive, D.N. Schramm and M. Srednicki, _Gravitinos as the Cold Dark Matter in an $\Omega$ = 1 Universe_, _Nucl. Phys. B_ 255 (1985) 495.
  * [31] J. McDonald, _Thermally generated gauge singlet scalars as selfinteracting dark matter_ , _Phys.Rev.Lett._ 88 (2002) 091304 [hep-ph/0106249].
  * [32] L.J. Hall, K. Jedamzik, J. March-Russell and S.M. West, _Freeze-In Production of FIMP Dark Matter_ , _JHEP_ 03 (2010) 080 [0911.1120].
  * [33] Y. Mambrini, K.A. Olive, J. Quevillon and B. Zaldivar, _Gauge Coupling Unification and Nonequilibrium Thermal Dark Matter_ , _Phys. Rev. Lett._ 110 (2013) 241306 [1302.4438].
  * [34] N. Bernal, M. Heikinheimo, T. Tenkanen, K. Tuominen and V. Vaskonen, _The Dawn of FIMP Dark Matter: A Review of Models and Constraints_ , _Int. J. Mod. Phys._ A32 (2017) 1730023 [1706.07442].
  * [35] F. Elahi, C. Kolda and J. Unwin, _UltraViolet Freeze-in_ , _JHEP_ 03 (2015) 048 [1410.6157].
  * [36] Y. Mambrini, N. Nagata, K.A. Olive, J. Quevillon and J. Zheng, _Dark matter and gauge coupling unification in nonsupersymmetric SO(10) grand unified models_ , _Phys. Rev. D_ 91 (2015) 095010 [1502.06929].
  * [37] G. Bhattacharyya, M. Dutra, Y. Mambrini and M. Pierre, _Freezing-in dark matter through a heavy invisible Z’_ , _Phys. Rev. D_ 98 (2018) 035038 [1806.00016].
  * [38] D. Chowdhury, E. Dudas, M. Dutra and Y. Mambrini, _Moduli Portal Dark Matter_ , _Phys. Rev. D_ 99 (2019) 095028 [1811.01947].
  * [39] K. Benakli, Y. Chen, E. Dudas and Y. Mambrini, _Minimal model of gravitino dark matter_ , _Phys. Rev. D_ 95 (2017) 095002 [1701.06574].
  * [40] E. Dudas, Y. Mambrini and K. Olive, _Case for an EeV Gravitino_ , _Phys. Rev. Lett._ 119 (2017) 051801 [1704.03008].
  * [41] E. Dudas, T. Gherghetta, Y. Mambrini and K.A. Olive, _Inflation and High-Scale Supersymmetry with an EeV Gravitino_ , _Phys. Rev. D_ 96 (2017) 115032 [1710.07341].
  * [42] E. Dudas, T. Gherghetta, K. Kaneta, Y. Mambrini and K.A. Olive, _Gravitino decay in high scale supersymmetry with R -parity violation_ , _Phys. Rev. D_ 98 (2018) 015030 [1805.07342].
  * [43] K. Kaneta, Y. Mambrini, K.A. Olive and S. Verner, _Inflation and Leptogenesis in High-Scale Supersymmetry_ , _Phys. Rev. D_ 101 (2020) 015002 [1911.02463].
  * [44] N. Bernal, M. Dutra, Y. Mambrini, K. Olive, M. Peloso and M. Pierre, _Spin-2 Portal Dark Matter_ , _Phys. Rev. D_ 97 (2018) 115020 [1803.01866].
  * [45] K. Kaneta, Y. Mambrini and K.A. Olive, _Radiative production of nonthermal dark matter_ , _Phys. Rev. D_ 99 (2019) 063508 [1901.04449].
  * [46] L. Heurtier and F. Huang, _Inflaton portal to a highly decoupled EeV dark matter particle_ , _Phys. Rev. D_ 100 (2019) 043507 [1905.05191].
  * [47] Y. Mambrini, K.A. Olive and J. Zheng, _Post-Inflationary Dark Matter Bremsstrahlung_ , 2208.05859.
  * [48] D.J. Chung, E.W. Kolb and A. Riotto, _Production of massive particles during reheating_ , _Phys. Rev. D_ 60 (1999) 063504 [hep-ph/9809453].
  * [49] G.F. Giudice, E.W. Kolb and A. Riotto, _Largest temperature of the radiation era and its cosmological implications_ , _Phys. Rev. D_ 64 (2001) 023508 [hep-ph/0005123].
  * [50] M.A.G. Garcia, Y. Mambrini, K.A. Olive and M. Peloso, _Enhancement of the Dark Matter Abundance Before Reheating: Applications to Gravitino Dark Matter_ , _Phys. Rev. D_ 96 (2017) 103510 [1709.01549].
  * [51] M.A.G. Garcia, K. Kaneta, Y. Mambrini and K.A. Olive, _Reheating and Post-inflationary Production of Dark Matter_ , _Phys. Rev. D_ 101 (2020) 123507 [2004.08404].
  * [52] M.A.G. Garcia, K. Kaneta, Y. Mambrini and K.A. Olive, _Inflaton Oscillations and Post-Inflationary Reheating_ , _JCAP_ 04 (2021) 012 [2012.10756].
  * [53] B. Barman, N. Bernal, Y. Xu and O. Zapata, _Ultraviolet freeze-in with a time-dependent inflaton decay_ , _JCAP_ 07 (2022) 019 [2202.12906].
  * [54] Y. Ema, R. Jinno, K. Mukaida and K. Nakayama, _Gravitational Effects on Inflaton Decay_ , _JCAP_ 05 (2015) 038 [1502.02475].
  * [55] M. Garny, M. Sandora and M.S. Sloth, _Planckian Interacting Massive Particles as Dark Matter_ , _Phys. Rev. Lett._ 116 (2016) 101302 [1511.03278].
  * [56] M. Garny, A. Palessandro, M. Sandora and M.S. Sloth, _Theory and Phenomenology of Planckian Interacting Massive Particles as Dark Matter_ , _JCAP_ 02 (2018) 027 [1709.09688].
  * [57] Y. Tang and Y.-L. Wu, _Pure Gravitational Dark Matter, Its Mass and Signatures_ , _Phys. Lett. B_ 758 (2016) 402 [1604.04701].
  * [58] Y. Tang and Y.-L. Wu, _On Thermal Gravitational Contribution to Particle Production and Dark Matter_ , _Phys. Lett. B_ 774 (2017) 676 [1708.05138].
  * [59] Y. Ema, R. Jinno, K. Mukaida and K. Nakayama, _Gravitational particle production in oscillating backgrounds and its cosmological implications_ , _Phys. Rev. D_ 94 (2016) 063517 [1604.08898].
  * [60] Y. Ema, K. Nakayama and Y. Tang, _Production of Purely Gravitational Dark Matter_ , _JHEP_ 09 (2018) 135 [1804.07471].
  * [61] Y. Ema, K. Nakayama and Y. Tang, _Production of purely gravitational dark matter: the case of fermion and vector boson_ , _JHEP_ 07 (2019) 060 [1903.10973].
  * [62] M. Chianese, B. Fu and S.F. King, _Impact of Higgs portal on gravity-mediated production of superheavy dark matter_ , _JCAP_ 06 (2020) 019 [2003.07366].
  * [63] M. Chianese, B. Fu and S.F. King, _Interplay between neutrino and gravity portals for FIMP dark matter_ , _JCAP_ 01 (2021) 034 [2009.01847].
  * [64] M. Redi, A. Tesi and H. Tillim, _Gravitational Production of a Conformal Dark Sector_ , _JHEP_ 05 (2021) 010 [2011.10565].
  * [65] Y. Mambrini and K.A. Olive, _Gravitational Production of Dark Matter during Reheating_ , _Phys. Rev. D_ 103 (2021) 115009 [2102.06214].
  * [66] B. Barman and N. Bernal, _Gravitational SIMPs_ , _JCAP_ 06 (2021) 011 [2104.10699].
  * [67] M.R. Haque and D. Maity, _Gravitational dark matter: Free streaming and phase space distribution_ , _Phys. Rev. D_ 106 (2022) 023506 [2112.14668].
  * [68] S. Clery, Y. Mambrini, K.A. Olive and S. Verner, _Gravitational portals in the early Universe_ , _Phys. Rev. D_ 105 (2022) 075005 [2112.15214].
  * [69] S. Clery, Y. Mambrini, K.A. Olive, A. Shkerin and S. Verner, _Gravitational portals with nonminimal couplings_ , _Phys. Rev. D_ 105 (2022) 095042 [2203.02004].
  * [70] A. Ahmed, B. Grzadkowski and A. Socha, _Higgs Boson-Induced Reheating and Dark Matter Production_ , _Symmetry_ 14 (2022) 306.
  * [71] A. Ahmed, B. Grzadkowski and A. Socha, _Higgs boson induced reheating and ultraviolet frozen-in dark matter_ , 2207.11218.
  * [72] M. Fukugita and T. Yanagida, _Baryogenesis Without Grand Unification_ , _Phys. Lett._ B174 (1986) 45.
  * [73] V.A. Kuzmin, V.A. Rubakov and M.E. Shaposhnikov, _On the Anomalous Electroweak Baryon Number Nonconservation in the Early Universe_ , _Phys. Lett._ 155B (1985) 36.
  * [74] W. Buchmuller, P. Di Bari and M. Plumacher, _Cosmic microwave background, matter - antimatter asymmetry and neutrino masses_ , _Nucl. Phys. B_ 643 (2002) 367 [hep-ph/0205349].
  * [75] W. Buchmuller, P. Di Bari and M. Plumacher, _The Neutrino mass window for baryogenesis_ , _Nucl. Phys. B_ 665 (2003) 445 [hep-ph/0302092].
  * [76] P.H. Chankowski and K. Turzynski, _Limits on T(reh) for thermal leptogenesis with hierarchical neutrino masses_ , _Phys. Lett. B_ 570 (2003) 198 [hep-ph/0306059].
  * [77] G. Giudice, A. Notari, M. Raidal, A. Riotto and A. Strumia, _Towards a complete theory of thermal leptogenesis in the SM and MSSM_ , _Nucl. Phys. B_ 685 (2004) 89 [hep-ph/0310123].
  * [78] S. Davidson and A. Ibarra, _A Lower bound on the right-handed neutrino mass from leptogenesis_ , _Phys. Lett._ B535 (2002) 25 [hep-ph/0202239].
  * [79] G. Giudice, M. Peloso, A. Riotto and I. Tkachev, _Production of massive fermions at preheating and leptogenesis_ , _JHEP_ 08 (1999) 014 [hep-ph/9905242].
  * [80] T. Asaka, K. Hamaguchi, M. Kawasaki and T. Yanagida, _Leptogenesis in inflaton decay_ , _Phys. Lett. B_ 464 (1999) 12 [hep-ph/9906366].
  * [81] G. Lazarides and Q. Shafi, _Origin of matter in the inflationary cosmology_ , _Phys. Lett. B_ 258 (1991) 305.
  * [82] B.A. Campbell, S. Davidson and K.A. Olive, _Inflation, neutrino baryogenesis, and (S)neutrino induced baryogenesis_ , _Nucl. Phys. B_ 399 (1993) 111 [hep-ph/9302223].
  * [83] F. Hahn-Woernle and M. Plumacher, _Effects of reheating on leptogenesis_ , _Nucl. Phys. B_ 806 (2009) 68 [0801.3972].
  * [84] R.T. Co, Y. Mambrini and K.A. Olive, _Inflationary Gravitational Leptogenesis_ , 2205.01689.
  * [85] N. Bernal and C.S. Fong, _Dark matter and leptogenesis from gravitational production_ , _JCAP_ 06 (2021) 028 [2103.06896].
  * [86] M.R. Haque and D. Maity, _Gravitational Reheating_ , 2201.02348.
  * [87] D.G. Figueroa and E.H. Tanin, _Inconsistency of an inflationary sector coupled only to Einstein gravity_ , _JCAP_ 10 (2019) 050 [1811.04093].
  * [88] R. Kallosh and A. Linde, _Universality Class in Conformal Inflation_ , _JCAP_ 07 (2013) 002 [1306.5220].
  * [89] S.Y. Choi, J.S. Shim and H.S. Song, _Factorization and polarization in linearized gravity_ , _Phys. Rev. D_ 51 (1995) 2751 [hep-th/9411092].
  * [90] B.R. Holstein, _Graviton Physics_ , _Am. J. Phys._ 74 (2006) 1002 [gr-qc/0607045].
  * [91] J. Ellis, M.A.G. Garcia, D.V. Nanopoulos and K.A. Olive, _Calculations of Inflaton Decays and Reheating: with Applications to No-Scale Inflation Models_ , _JCAP_ 07 (2015) 050 [1505.06986].
  * [92] Planck collaboration, _Planck 2018 results. X. Constraints on inflation_ , _Astron. Astrophys._ 641 (2020) A10 [1807.06211].
  * [93] J. Ellis, M.A.G. Garcia, D.V. Nanopoulos, K.A. Olive and S. Verner, _BICEP/Keck constraints on attractor models of inflation and reheating_ , _Phys. Rev. D_ 105 (2022) 043504 [2112.04466].
  * [94] P. Minkowski, _$\mu\to e\gamma$ at a rate of one out of 109 muon decays?_, _Physics Letters B_ 67 (1977) 421.
  * [95] M. Gell-Mann, P. Ramond and R. Slansky, _Complex Spinors and Unified Theories_ , _Conf. Proc._ C790927 (1979) 315 [1306.4669].
  * [96] T. Yanagida, _HORIZONTAL SYMMETRY AND MASSES OF NEUTRINOS_ , _Conf. Proc._ C7902131 (1979) 95.
  * [97] R.N. Mohapatra and G. Senjanovic, _Neutrino Mass and Spontaneous Parity Violation_ , _Phys. Rev. Lett._ 44 (1980) 912.
  * [98] J. Schechter and J.W.F. Valle, _Neutrino Masses in SU(2) x U(1) Theories_ , _Phys. Rev._ D22 (1980) 2227.
  * [99] J. Schechter and J.W.F. Valle, _Neutrino Decay and Spontaneous Violation of Lepton Number_ , _Phys. Rev. D_ 25 (1982) 774.
  * [100] P.H. Frampton, S.L. Glashow and T. Yanagida, _Cosmological sign of neutrino CP violation_ , _Phys. Lett. B_ 548 (2002) 119 [hep-ph/0208157].
  * [101] M. Raidal and A. Strumia, _Predictions of the most minimal seesaw model_ , _Phys. Lett. B_ 553 (2003) 72 [hep-ph/0210021].
  * [102] A. Ibarra and G.G. Ross, _Neutrino phenomenology: The Case of two right-handed neutrinos_ , _Phys. Lett. B_ 591 (2004) 285 [hep-ph/0312138].
  * [103] T. Rink, _Leptonic CP violation in the minimal type-I seesaw model: Bottom-up phenomenology $\&$ top-down model building_, Master’s thesis, Heidelberg U., 2017.
  * [104] K. Ichikawa, T. Suyama, T. Takahashi and M. Yamaguchi, _Primordial Curvature Fluctuation and Its Non-Gaussianity in Models with Modulated Reheating_ , _Phys. Rev. D_ 78 (2008) 063545 [0807.3988].
  * [105] K. Kainulainen, S. Nurmi, T. Tenkanen, K. Tuominen and V. Vaskonen, _Isocurvature Constraints on Portal Couplings_ , _JCAP_ 06 (2016) 022 [1601.07733].
  * [106] Y. Mambrini, _Particles in the dark Universe_ , _Springer Ed., ISBN 978-3-030-78139-2_ (2021) .
  * [107] K.D. Lozanov and M.A. Amin, _Equation of State and Duration to Radiation Domination after Inflation_ , _Phys. Rev. Lett._ 119 (2017) 061301 [1608.01213].
  * [108] K.D. Lozanov and M.A. Amin, _Self-resonance after inflation: oscillons, transients and radiation domination_ , _Phys. Rev. D_ 97 (2018) 023533 [1710.06851].
  * [109] D. Maity and P. Saha, _(P)reheating after minimal Plateau Inflation and constraints from CMB_ , _JCAP_ 07 (2019) 018 [1811.11173].
  * [110] M.A. Luty, _Baryogenesis via leptogenesis_ , _Phys. Rev. D_ 45 (1992) 455.
  * [111] M. Flanz, E.A. Paschos and U. Sarkar, _Baryogenesis from a lepton asymmetric universe_ , _Phys. Lett. B_ 345 (1995) 248 [hep-ph/9411366].
  * [112] L. Covi, E. Roulet and F. Vissani, _CP violating decays in leptogenesis scenarios_ , _Phys. Lett. B_ 384 (1996) 169 [hep-ph/9605319].
  * [113] W. Buchmuller, P. Di Bari and M. Plumacher, _Leptogenesis for pedestrians_ , _Annals Phys._ 315 (2005) 305 [hep-ph/0401240].
  * [114] S. Davidson, E. Nardi and Y. Nir, _Leptogenesis_ , _Phys. Rept._ 466 (2008) 105 [0802.2962].
  * [115] T. Asaka, H. Nielsen and Y. Takanishi, _Nonthermal leptogenesis from the heavier Majorana neutrinos_ , _Nucl. Phys. B_ 647 (2002) 252 [hep-ph/0207023].
  * [116] V.N. Senoguz and Q. Shafi, _GUT scale inflation, nonthermal leptogenesis, and atmospheric neutrino oscillations_ , _Phys. Lett. B_ 582 (2004) 6 [hep-ph/0309134].
  * [117] B. Barman, D. Borah and R. Roshan, _Nonthermal leptogenesis and UV freeze-in of dark matter: Impact of inflationary reheating_ , _Phys. Rev. D_ 104 (2021) 035022 [2103.01675].
  * [118] T.-H. Yeh, J. Shelton, K.A. Olive and B.D. Fields, _Probing Physics Beyond the Standard Model: Limits from BBN and the CMB Independently and Combined_ , 2207.13133.
  * [119] R.T. Co, D. Dunsky, N. Fernandez, A. Ghalsasi, L.J. Hall, K. Harigaya et al., _Gravitational wave and CMB probes of axion kination_ , _JHEP_ 09 (2022) 116 [2108.09299].
  * [120] K. Saikawa and S. Shirai, _Primordial gravitational waves, precisely: The role of thermodynamics in the Standard Model_ , _JCAP_ 05 (2018) 035 [1803.01038].
  * [121] T. Opferkuch, P. Schwaller and B.A. Stefanek, _Ricci Reheating_ , _JCAP_ 07 (2019) 016 [1905.06823].
  * [122] M. Giovannini, _Gravitational waves constraints on postinflationary phases stiffer than radiation_ , _Phys. Rev. D_ 58 (1998) 083504 [hep-ph/9806329].
  * [123] J. Crowder and N.J. Cornish, _Beyond LISA: Exploring future gravitational wave missions_ , _Phys. Rev. D_ 72 (2005) 083005 [gr-qc/0506015].
  * [124] V. Corbin and N.J. Cornish, _Detecting the cosmic gravitational wave background with the big bang observer_ , _Class. Quant. Grav._ 23 (2006) 2435 [gr-qc/0512039].
  * [125] G.M. Harry, P. Fritschel, D.A. Shaddock, W. Folkner and E.S. Phinney, _Laser interferometry for the big bang observer_ , _Class. Quant. Grav._ 23 (2006) 4887.
  * [126] N. Seto, S. Kawamura and T. Nakamura, _Possibility of direct measurement of the acceleration of the universe using 0.1-Hz band laser interferometer gravitational wave antenna in space_ , _Phys. Rev. Lett._ 87 (2001) 221103 [astro-ph/0108011].
  * [127] S. Kawamura et al., _The Japanese space gravitational wave antenna DECIGO_ , _Class. Quant. Grav._ 23 (2006) S125.
  * [128] K. Yagi and N. Seto, _Detector configuration of DECIGO/BBO and identification of cosmological neutron-star binaries_ , _Phys. Rev. D_ 83 (2011) 044011 [1101.3940].
  * [129] LIGO Scientific collaboration, _Exploring the Sensitivity of Next Generation Gravitational Wave Detectors_ , _Class. Quant. Grav._ 34 (2017) 044001 [1607.08697].
  * [130] D. Reitze et al., _Cosmic Explorer: The U.S. Contribution to Gravitational-Wave Astronomy beyond LIGO_ , _Bull. Am. Astron. Soc._ 51 (2019) 035 [1907.04833].
  * [131] M. Punturo et al., _The Einstein Telescope: A third-generation gravitational wave observatory_ , _Class. Quant. Grav._ 27 (2010) 194002.
  * [132] S. Hild et al., _Sensitivity Studies for Third-Generation Gravitational Wave Observatories_ , _Class. Quant. Grav._ 28 (2011) 094013 [1012.0908].
  * [133] B. Sathyaprakash et al., _Scientific Objectives of Einstein Telescope_ , _Class. Quant. Grav._ 29 (2012) 124013 [1206.0331].
  * [134] M. Maggiore et al., _Science Case for the Einstein Telescope_ , _JCAP_ 03 (2020) 050 [1912.02622].
  * [135] K. Schmitz, _New Sensitivity Curves for Gravitational-Wave Experiments_ , 2002.04615.
  * [136] IceCube collaboration, _First observation of PeV-energy neutrinos with IceCube_ , _Phys. Rev. Lett._ 111 (2013) 021103 [1304.5356].
  * [137] IceCube collaboration, _Evidence for High-Energy Extraterrestrial Neutrinos at the IceCube Detector_ , _Science_ 342 (2013) 1242856 [1311.5238].
  * [138] IceCube collaboration, _Observation of High-Energy Astrophysical Neutrinos in Three Years of IceCube Data_ , _Phys. Rev. Lett._ 113 (2014) 101101 [1405.5303].
  * [139] Y. Bai, R. Lu and J. Salvado, _Geometric Compatibility of IceCube TeV-PeV Neutrino Excess and its Galactic Dark Matter Origin_ , _JHEP_ 01 (2016) 161 [1311.5864].
  * [140] T. Higaki, R. Kitano and R. Sato, _Neutrinoful Universe_ , _JHEP_ 07 (2014) 044 [1405.0013].
  * [141] A. Esmaili, S.K. Kang and P.D. Serpico, _IceCube events and decaying dark matter: hints and constraints_ , _JCAP_ 12 (2014) 054 [1410.5979].
  * [142] A. Bhattacharya, M.H. Reno and I. Sarcevic, _Reconciling neutrino flux from heavy dark matter decay and recent events at IceCube_ , _JHEP_ 06 (2014) 110 [1403.1862].
  * [143] J. Zavala, _Galactic PeV neutrinos from dark matter annihilation_ , _Phys. Rev. D_ 89 (2014) 123516 [1404.2932].
  * [144] C. Rott, K. Kohri and S.C. Park, _Superheavy dark matter and IceCube neutrino signals: Bounds on decaying dark matter_ , _Phys. Rev. D_ 92 (2015) 023529 [1408.4575].
  * [145] E. Dudas, Y. Mambrini and K.A. Olive, _Monochromatic neutrinos generated by dark matter and the seesaw mechanism_ , _Phys. Rev. D_ 91 (2015) 075001 [1412.3459].
  * [146] K. Murase, R. Laha, S. Ando and M. Ahlers, _Testing the Dark Matter Scenario for PeV Neutrinos Observed in IceCube_ , _Phys. Rev. Lett._ 115 (2015) 071301 [1503.04663].
  * [147] L.A. Anchordoqui, V. Barger, H. Goldberg, X. Huang, D. Marfatia, L.H.M. da Silva et al., _IceCube neutrinos, decaying dark matter, and the Hubble constant_ , _Phys. Rev. D_ 92 (2015) 061301 [1506.08788].
  * [148] P. Ko and Y. Tang, _IceCube Events from Heavy DM decays through the Right-handed Neutrino Portal_ , _Phys. Lett. B_ 751 (2015) 81 [1508.02500].
  * [149] S.B. Roland, B. Shakya and J.D. Wells, _PeV neutrinos and a 3.5 keV x-ray line from a PeV-scale supersymmetric neutrino sector_ , _Phys. Rev. D_ 92 (2015) 095018 [1506.08195].
  * [150] M. Chianese and A. Merle, _A Consistent Theory of Decaying Dark Matter Connecting IceCube to the Sesame Street_ , _JCAP_ 04 (2017) 017 [1607.05283].
  * [151] M. Re Fiorentin, V. Niro and N. Fornengo, _A consistent model for leptogenesis, dark matter and the IceCube signal_ , _JHEP_ 11 (2016) 022 [1606.04445].
  * [152] M.A.G. Garcia, Y. Mambrini, K.A. Olive and S. Verner, _Case for decaying spin- 3/2 dark matter_ , _Phys. Rev. D_ 102 (2020) 083533 [2006.03325].
  * [153] C.A. Argüelles, D. Delgado, A. Friedlander, A. Kheirandish, I. Safa, A.C. Vincent et al., _Dark Matter decay to neutrinos_ , 2210.01303.
  * [154] P. Jizba and G. Lambiase, _Tsallis cosmology and its applications in dark matter physics with focus on IceCube high-energy neutrino data_ , 2206.12910.
  * [155] M. Ahmadvand, _Filtered asymmetric dark matter during the Peccei-Quinn phase transition_ , _JHEP_ 10 (2021) 109 [2108.00958].
  * [156] Y. Chikashige, R.N. Mohapatra and R.D. Peccei, _Are There Real Goldstone Bosons Associated with Broken Lepton Number?_ , _Phys. Lett. B_ 98 (1981) 265.
  * [157] V. Berezinsky and J.W.F. Valle, _The KeV majoron as a dark matter particle_ , _Phys. Lett. B_ 318 (1993) 360 [hep-ph/9309214].
  * [158] M. Lattanzi and J.W.F. Valle, _Decaying warm dark matter and neutrino masses_ , _Phys. Rev. Lett._ 99 (2007) 121301 [0705.2406].
  * [159] F. Bazzocchi, M. Lattanzi, S. Riemer-Sørensen and J.W.F. Valle, _X-ray photons from late-decaying majoron dark matter_ , _JCAP_ 08 (2008) 013 [0805.2372].
  * [160] M. Lattanzi, S. Riemer-Sorensen, M. Tortola and J.W.F. Valle, _Updated CMB and x- and $\gamma$-ray constraints on Majoron dark matter_, _Phys. Rev. D_ 88 (2013) 063528 [1303.4685].
  * [161] M. Lattanzi, R.A. Lineros and M. Taoso, _Connecting neutrino physics with dark matter_ , _New J. Phys._ 16 (2014) 125012 [1406.0004].
  * [162] J. Gehrlein and M. Pierre, _A testable hidden-sector model for Dark Matter and neutrino masses_ , _JHEP_ 02 (2020) 068 [1912.06661].
  * [163] E. Dudas, L. Heurtier, Y. Mambrini, K.A. Olive and M. Pierre, _Model of metastable EeV dark matter_ , _Phys. Rev. D_ 101 (2020) 115029 [2003.02846].
  * [164] V. Iršič et al., _New Constraints on the free-streaming of warm dark matter from intermediate and small scale Lyman- $\alpha$ forest data_, _Phys. Rev. D_ 96 (2017) 023522 [1702.01764].
  * [165] M. Sitwell, A. Mesinger, Y.-Z. Ma and K. Sigurdson, _The Imprint of Warm Dark Matter on the Cosmological 21-cm Signal_ , _Mon. Not. Roy. Astron. Soc._ 438 (2014) 2664 [1310.0029].