# Gravitational Portals with Non-Minimal Couplings

Simon Clérya<EMAIL_ADDRESS>Yann Mambrinia,b
<EMAIL_ADDRESS>Keith A. Olivec<EMAIL_ADDRESS>Andrey
Shkerinc<EMAIL_ADDRESS>Sarunas Vernerc<EMAIL_ADDRESS>a Université Paris-
Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France  b CERN, Theoretical Physics
Department, Geneva, Switzerland  c William I. Fine Theoretical Physics
Institute, School of Physics and Astronomy, University of Minnesota,
Minneapolis, MN 55455, USA

###### Abstract

We consider the effects of non-minimal couplings to curvature of the form
$\xi_{S}S^{2}R$, for three types of scalars: the Higgs boson, the inflaton,
and a scalar dark matter candidate. We compute the abundance of dark matter
produced by these non-minimal couplings to gravity and compare to similar
results with minimal couplings. We also compute the contribution to the
radiation bath during reheating. The main effect is a potential augmentation
of the maximum temperature during reheating. A model independent limit of
$\mathcal{O}(10^{12})$ GeV is obtained. For couplings
$\xi_{S}\gtrsim\mathcal{O}(1)$, these dominate over minimal gravitational
interactions.

††preprint: UMN–TH–4116/22††preprint: FTPI–MINN–22/07††preprint: CERN-
TH-2022-025

##  I Introduction

Promoting a field theory Lagrangian from a Lorentz-invariant one to a
generally-covariant one necessarily leads to an interaction between the fields
of the theory and the gravitational field. In the case of a scalar field, $S$,
the natural generalization of this minimal interaction scenario is to
introduce a non-minimal coupling term of the form

$\propto\xi_{S}S^{2}R\;.$ (1)

Here $R$ is the Ricci scalar and $\xi_{S}$ is a non-minimal coupling constant.
This non-minimal coupling to gravity proved to be useful in many applications
to cosmology. Examples include Higgs inflation Bezrukov:2007ep ;
Lebedev:2021xey , where $S$ is associated with the Higgs field degree of
freedom $h$ — the only scalar degree of freedom in the Standard Model,
preheating Ema:2016dny , where $S$ is associated with the inflaton field
$\phi$, and non-perturbative production of dark matter nonminprod , where $S$
represents the scalar dark matter particle $X$.

In the general case, when the fields $\phi$, $h$, and $X$ are all different,
the question arises as to what extent they must interact with each other in
order to successfully reheat the Universe and generate the right amount of
dark matter. Recent studies have shown that interactions via gravity alone, to
which the fields are coupled minimally, is enough for these purposes. Indeed,
the perturbative gravitational production of dark matter through graviton
exchange can play a dominant role during reheating with processes involving
the inflaton MO ; CMOV ; Barman:2021ugy  as well as thermal bath particles
CMOV ; Haque:2021mab . Further, the minimal gravitational coupling can lead to
the completion of the reheating process for certain types of the inflationary
potential, $V(\phi)\sim\phi^{k}$ with $k>2$ CMOV ; Haque:2022kez . Thus,
gravity is strong enough to mediate perturbative channels of reheating and
dark matter production.

The purpose of this work is to study how the inclusion of the non-minimal
coupling terms of the form (1) affect the gravitational production of dark
matter and radiation during reheating. Note that the presence of these terms
is unavoidable: if there were no such couplings at tree level, they would
still be generated by quantum corrections Callan:1970ze . We study particle
production in the processes $hh\rightarrow XX$, $\phi\phi\rightarrow hh$, and
$\phi\phi\rightarrow XX$ which are induced by the non-minimal couplings. Here
$\phi$ represents the inflaton background oscillating around its minimum after
the end of inflation gravprod . Since the scalar fields couple directly to the
curvature scalar $R$, the oscillating background causes the effective masses
of the fields to change non-adiabatically and leads to particle production.
This regime of particle creation has been considered in several different
contexts, including gravitational production of scalar gravscalar ; ema ,
fermion gravferm , and vector dark matter gravvector .

Our main interest is to compare the (dark) matter production channels induced
by the non-minimal couplings with the production via the s-channel graviton
exchange that sets minimal possible production rates. We will see for which
values of the couplings the rates are enhanced, and what are the consequences
on the dark matter density or the temperature attained during reheating.
Throughout the work we adopt the Starobinsky inflationary potential staro ,
although our results are largely independent of the particular form of the
potential. As for the potentials for the fields $h$ and $X$, we take them to
be renormalizable polynomials. We also assume no direct interaction between
$\phi$, $h$, and $X$.

Working in the perturbative regime implies that the non-minimal couplings must
satisfy $|\xi_{S}|\ll M_{P}^{2}/\langle S\rangle^{2}$, where $\langle
S\rangle$ is the vacuum expectation value of $S=\phi,h,X$. The value of
$\xi_{h}$ is constrained from collider experiments as $|\xi_{h}|\lesssim
10^{15}$ higgscons1 .111Note that in the case of Higgs inflation, $\xi_{h}$ is
fixed from CMB measurements Bezrukov:2007ep . Furthermore, the lower bound on
$\xi_{h}$ comes from the fact that the Standard Model electroweak vacuum may
not be absolutely stable HiggsStab . To prevent the vacuum decay due to
quantum fluctuations during inflation fluc , the effective mass of the Higgs
field induced by the non-minimal coupling must be large enough; this gives
$\xi_{h}\gtrsim 10^{-1}$ higgstree ; higgsloop  (see also Markkanen:2018pdo
).222This estimate assumes no new physics interfering the RG running of the
Higgs self-coupling constant until inflationary energy scales.

The paper is organized as follows: The framework for our computation is
presented in Section II. We discuss non-minimal gravitational couplings of the
inflaton, the Higgs boson, and a dark matter scalar in detail. We calculate
the dark matter production rates either from scattering in the thermal bath or
from oscillations in the inflaton condensate. We compare similar processes
obtained from the minimal gravitational particle production. We choose the
Starobinsky model of inflation and discuss the reheating epoch when the
inflaton begins oscillating. In Section III we discuss the resulting abundance
of dark matter produced from the thermal bath and directly from scattering of
the inflaton condensate. We also compute the effects of the non-minimal
couplings on the maximum temperature attained during reheating. We then
compare different processes in Section IV, before summarizing our results in
Section V.

##  II The framework

###  II.1 Scalar-gravity Lagrangian

The theory we consider comprises 3 scalar fields non-minimally coupled to
gravity: the inflaton $\phi$, the Higgs field333We consider the Higgs boson as
a surrogate for any additional scalars with Standard Model couplings. $H$, for
which we adopt the Unitary gauge, $H=(0,h)^{T}/\sqrt{2}$, and the dark matter
candidate $X$. The relevant part of the action takes the form444The metric
signature is chosen as $(+,-,-,-)$.

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

with the conformal factor $\Omega^{2}$ given by

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

Here $M_{P}=2.4\times 10^{18}\,\rm{GeV}$ is the reduced Planck mass, and the
tilde used in Eq. (2) indicates that the theory is considered in the Jordan
frame. For the scalar field Lagrangians we have

$\mathcal{L}_{S}\;=\frac{1}{2}\tilde{g}^{\mu\nu}\partial_{\mu}S\partial_{\nu}S-V_{S}\,,~{}~{}~{}S=\phi,h,X\,.$
(4)

Next, we specify the scalar field potentials. For a model of inflation, we
choose the well-motivated Starobinsky model for which staro

$V_{\phi}\;=\;\frac{3}{4}m_{\phi}^{2}M_{P}^{2}\left(1-e^{-\sqrt{\frac{2}{3}}\frac{\phi}{M_{P}}}\right)^{2}\,.$
(5)

In what follows, we work in the perturbative regime with $\phi\ll M_{P}$,
hence the potential is approximated as

$V_{\phi}\;\simeq\;\frac{1}{2}m_{\phi}^{2}\phi^{2}\,.$ (6)

The inflaton mass, $m_{\phi}$, is fixed by the amplitude of scalar
perturbations inferred from CMB measurements Planck ; for the potential (5)
this gives $m_{\phi}=3\times 10^{13}$ GeV building .

The potential for the Higgs field is taken as follows

$V_{h}\;=\;\frac{1}{2}m_{h}^{2}h^{2}+\frac{1}{4}\lambda_{h}h^{4}\,.$ (7)

Here $m_{h}$ and $\lambda_{h}$ are the Higgs mass and quartic coupling,
correspondingly. Note that both parameters undergo the renormalization group
(RG) running. In what follows we take a weak scale mass, which is a good
approximation at the time of reheating and our results are insensitive to
$\lambda_{h}$. Finally, the dark matter potential is simply given by

$V_{X}\;=\;\frac{1}{2}m_{X}^{2}X^{2}\,.$ (8)

To study the reheating in the theory (2), it is convenient to remove the non-
minimal couplings by performing the redefinition of the metric field. Leaving
the details to Appendix A, we write the action (2) in the Einstein frame,

$\begin{split}\mathcal{S}\;=\;\int
d^{4}x\sqrt{-g}&\left[-\frac{M_{P}^{2}}{2}R+\frac{1}{2}K^{ij}g^{\mu\nu}\partial_{\mu}S_{i}\partial_{\nu}S_{j}\right.\\\
&\qquad\qquad\qquad\left.-\frac{V_{\phi}+V_{h}+V_{X}}{\Omega^{4}}\right]\,.\end{split}$
(9)

Here the indices $i,j$ enumerate the fields $\phi,h,X$, and the kinetic
function is given by

$K^{ij}\;=\;6\frac{\partial\log\Omega}{\partial
S_{i}}\frac{\partial\log\Omega}{\partial
S_{j}}+\frac{\delta^{ij}}{\Omega^{2}}\,.$ (10)

Note that the scalar field kinetic term is not canonical. In general, it is
impossible to make a field redefinition that would bring it to the canonical
form, unless all three non-minimal couplings vanish.555 Such a redefinition
exists if the three-dimensional manifold spanned by the fields $\phi$, $h$ and
$X$ is flat. One can show that it is not the case if at least one of the
couplings is non-zero. For the theory (9) to be well-defined, the kinetic
function (10) must be positive-definite. Computing the eigenvalues, one
arrives at the condition

$\Omega^{2}>0\,,$ (11)

which is satisfied automatically for positive values of the couplings. Note
that the negative couplings are also allowed for certain scalar field
magnitudes.

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}}\;,~{}~{}\frac{|\xi_{X}|X^{2}}{M_{P}^{2}}\ll
1\,.$ (12)

We can expand the kinetic and potential terms in the action (9) in powers of
$M_{P}^{-2}$. We obtain a canonical kinetic term for the scalar fields and
deduce the leading-order interactions induced by the non-minimal couplings.
The latter can be brought to the form

$\mathcal{L}_{\rm{non-min.}}\;=\;-\sigma_{hX}^{\xi}h^{2}X^{2}-\sigma_{\phi
X}^{\xi}\phi^{2}X^{2}-\sigma_{\phi h}^{\xi}\phi^{2}h^{2}\,,$ (13)

where the $\sigma_{ij}^{\xi}$ are functions of the couplings $\xi_{i}$,
$\xi_{j}$, the masses $m_{i}$, $m_{j}$, and the Mandelstam variables; see
Appendix A for details.

The small-field approximation (12) implies the bound $\sqrt{|\xi_{S}|}\lesssim
M_{P}/\langle S\rangle$ with $S=\phi,h,X$. Since the inflaton value at the end
of inflation is $\phi_{\rm end}\sim M_{P}$ and afterwards
$\langle\phi^{2}\rangle\sim a^{-3}$, where $a$ is the cosmological scale
factor, then $|\xi_{\phi}|\lesssim(a/a_{\rm end})^{3}$. In particular, at the
onset of inflaton oscillations

$|\xi_{\phi}|\lesssim 1\,.$ (14)

Note that since our calculations involve the effective couplings $\sigma_{\phi
X}^{\xi}$ ($\sigma_{\phi h}^{\xi}$), which depend both on $\xi_{\phi}$ and
$\xi_{X}$ ($\xi_{h}$), the relatively small value of $|\xi_{\phi}|$ can, in
principle, be compensated by a large value of the other couplings.

In Fig. 1, we show the scattering processes obtained from the Lagrangian (13).
These contribute to reheating (when $h$ is in the final state) and dark matter
production (when $X$ is in the final state).

Figure 1: Feynman diagram for the 4-point interactions between the inflaton
$\phi$, the dark matter scalar candidate $X$, and the Higgs boson $h$, given
by the Lagrangian (13).

Finally, in evaluating the cosmological parameters, it is important to stay
within the validity of the low-energy theory. The cutoff of the theory can be
estimated as (see, e.g., higgscons2 )

$\Lambda\sim\frac{M_{P}}{\text{max}_{i}\>|\xi_{i}|}\,.$ (15)

In particular, the temperature of reheating must not exceed $\Lambda$.

###  II.2 Graviton exchange

Let us first consider the case of vanishing $\xi_{\phi,\,h,\,X}$, i.e., the
case of the minimal coupling of the scalar fields to gravity MO ; CMOV ; ema ;
Garny:2015sjg ; Tang:2017hvq ; Chianese:2020yjo ; Redi:2020ffc . It was argued
in MO ; CMOV  that the interaction between the dark and visible sectors
induced by gravity leads to unavoidable contributions to reheating and dark
matter production, in the thermal bath or via the scattering of the inflaton
condensate, through the graviton exchange processes shown in Fig. 2. It is
therefore important to compare the minimal gravitational particle production
to similar processes obtained from the Lagrangian in Eq. (13) with non-minimal
couplings.

Figure 2: Feynman diagram for the (dark) matter production through the
gravitational scattering of the inflaton or the Higgs boson from the thermal
bath.

To study the universal gravitational interactions in minimally coupled
gravity, we expand the space-time metric around flat space using
$g_{\mu\nu}\simeq\eta_{\mu\nu}+2h_{\mu\nu}/M_{P}$, where $h_{\mu\nu}$ is the
canonically-normalized perturbation. The gravitational interactions are
characterized by the following Lagrangian,

${\cal L}_{\rm
min.}=-\frac{1}{M_{P}}h_{\mu\nu}\left(T^{\mu\nu}_{h}+T^{\mu\nu}_{\phi}+T^{\mu\nu}_{X}\right)\,,$
(16)

where the stress-energy tensor is given by

$T^{\mu\nu}_{S}\;=\;\partial^{\mu}S\partial^{\nu}S-g^{\mu\nu}\left[\frac{1}{2}\partial^{\alpha}S\partial_{\alpha}S-V_{S}\right]\,.$
(17)

Note that in this work, we consider only the Higgs field in the visible
sector. Generalization to the complete spectrum of the Standard Model is
straightforward, and we leave it for future work.

For models with minimally coupled gravity, the processes
$\phi/h(p_{1})+\phi/h(p_{2})\rightarrow{h}/X(p_{3})+{h}/X(p_{4})$ can be
parametrized by

$\mathcal{M}^{00}\propto
M_{\mu\nu}^{0}\Pi^{\mu\nu\rho\sigma}M_{\rho\sigma}^{0}\;,$ (18)

where the graviton propagator for the canonically-normalized field
$h_{\mu\nu}$ with exchange momentum $k=p_{1}+p_{2}$ is given by

$\Pi^{\mu\nu\rho\sigma}(k)=\frac{\eta^{\mu\rho}\eta^{\nu\sigma}+\eta^{\mu\sigma}\eta^{\nu\rho}-\eta^{\mu\nu}\eta^{\rho\sigma}}{2k^{2}}\,,$
(19)

and the partial amplitude, $M_{\mu\nu}^{0}$, is given by

$\displaystyle
M_{\mu\nu}^{0}=\frac{1}{2}\left[p_{1\mu}p_{2\nu}+p_{1\nu}p_{2\mu}-\eta_{\mu\nu}p_{1}\cdot
p_{2}-\eta_{\mu\nu}V_{S}^{\prime\prime}\right]\,,$ (20)

with analogous expression for the final state in terms of outgoing momenta
$p_{3,4}$ and the final state potential. In Fig. 2 we show the s-channel
graviton exchange scattering obtained from the Lagrangian (16) for the
production of dark matter from either the Higgs field or the inflaton
condensate as well as the reheating process (the production of Higgs bosons
from the inflaton condensate).

###  II.3 Production rates

In this work, we consider three processes:

  1. A.

The production of dark matter from the scattering of thermal Higgs bosons
(assuming reheating is produced by inflaton decay). In this case, the dark
matter is populated via a freeze-in mechanism throughout the reheating period.

  2. B.

The production of dark matter from direct excitations of the inflaton
condensate. This process, which can be viewed as gravitational inflaton
scattering, is independent of the presence of a thermal bath.

  3. C.

The creation of a radiative bath at the start of reheating arising from the
Higgs boson production through gravitational inflaton scattering. Since such a
process is unavoidable in minimally coupled gravity, it is interesting to know
when such a process becomes dominant in models with non-minimal couplings
$\xi_{i}$.

The thermal dark matter production rate $R(T)$ for the process $hh\rightarrow
XX$ can be calculated from666We include the symmetry factors associated with
identical initial and final states in the definition of
$|\overline{{\cal{M}}}|^{2}$, and a factor of 2 is explicitly included in the
definition of the rate to account for the production of 2 identical particles.
gravitino

$R(T)=\frac{2\times N_{h}}{1024\pi^{6}}\int
f_{1}f_{2}E_{1}\mathop{}\\!\mathrm{d}E_{1}E_{2}\mathop{}\\!\mathrm{d}E_{2}\mathop{}\\!\mathrm{d}\cos\theta_{12}\int|\overline{{\cal
M}}|^{2}\mathop{}\\!\mathrm{d}\Omega_{13}\,,$ (21)

where $E_{i}$ is the energy of particle $i=1,2$, $\theta_{13}$ and
$\theta_{12}$ are the angles formed by momenta ${\bf{p}}_{1,3}$ and
${\bf{p}}_{1,2}$, respectively. $N_{h}=4$ is the number of internal degrees of
freedom for 1 complex Higgs doublet, $|\overline{{\cal M}}|^{2}$ is the matrix
amplitude squared with all symmetry factors included. This accounts for the
explicit factor of 2 in the numerator of Eq. (21). The thermal distribution
function of the incoming Higgs particles is given by the Bose-Einstein
distribution

$f_{i}=\frac{1}{e^{E_{i}/T}-1}\,.$ (22)

The rate for minimal gravitational interactions from Eq. (16) was derived in
CMOV ; Bernal:2018qlk . The rate we use here differs in two respects. As noted
earlier, we only include Higgs scalars in the initial state whereas in CMOV ;
Bernal:2018qlk , all Standard Model particle initial states were included.
Secondly, we keep terms depending on the dark matter mass which had not
previously been taken into account. This allows us to consider dark matter
masses approaching the inflaton mass and/or the reheating temperature.

For minimal (non-minimal) gravitational interactions, we find that the thermal
dark matter production rate can be expressed as

$\displaystyle R_{X}^{T,\,(\xi)}(T)$ $\displaystyle=$
$\displaystyle\beta_{1}^{(\xi)}\frac{T^{8}}{M_{P}^{4}}+\beta_{2}^{(\xi)}\frac{m_{X}^{2}T^{6}}{M_{P}^{4}}+\beta_{3}^{(\xi)}\frac{m_{X}^{4}T^{4}}{M_{P}^{4}}\,,$
(23)

where the coefficients $\beta_{1,\,2,\,3}^{(\xi)}$ are given in Appendix B by
Eqs. (84-86) (Eqs. (80-82)). The ratio of the non-minimal to minimal rate is
shown in Fig. 3. However, we note that when $\xi_{i}\sim\mathcal{O}(1)$ both
rates are comparable and interference effects become significant. The full
coefficients $\beta_{1,\,2,\,3}$ including interference are given by Eqs.
(87-89) from Appendix B. We leave the comparison of the effects on dark matter
production from the two rates for the next section.

Figure 3: Contours of the ratio of the dark matter production rates from the
thermal bath based on non-minimal gravitational interactions to those based on
minimal interactions. The ratio is displayed in the $(\xi_{h},\xi_{X})$ plane.
Note that as discussed in the Introduction, negative values of $\xi_{h}$ may
require new physics (such as supersymmetry) to stabilize the Higgs vacuum.

The rate for dark matter produced from inflaton oscillations of the inflaton
condensate for a potential of the form $V=\lambda\phi^{k}$ were considered in
detail in GKMO2 ; CMOV . The time-dependent inflaton can be written as
$\phi(t)=\phi_{0}(t)\mathcal{P}(t)$, where $\phi_{0}(t)$ is the time-dependent
amplitude that includes the effects of redshift and $\mathcal{P}(t)$ describes
the periodicity of the oscillation. The dark matter production rate is
calculated by writing the potential in terms of the Fourier modes of the
oscillations Ichikawa:2008ne ; Kainulainen:2016vzv ; GKMO2 ; CMOV

$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}\,.$
(24)

For $k=2$ (the only case considered here), the frequency of oscillation is
simply, $\omega=m_{\phi}$.

The rate generated by non-minimal couplings can be readily calculated using
the Lagrangian (13), which leads to

$R_{X}^{\phi,\,\xi}\;=\;\frac{2\times\sigma_{\phi
X}^{\xi~{}2}}{\pi}\frac{\rho_{\phi}^{2}}{m_{\phi}^{4}}\Sigma_{0}^{k}\,,$ (25)

where

$\Sigma_{0}^{k}=\sum_{n=1}^{\infty}|{\cal
P}^{k}_{n}|^{2}\sqrt{1-\frac{4m_{X}^{2}}{E_{n}^{2}}}\,,$ (26)

and $E_{n}=n\omega$ is the energy of the $n$-th inflaton oscillation mode. For
$k=2$, only the second Fourier mode in the sum contributes, with
$\sum|\mathcal{P}^{2}_{n}|^{2}=\frac{1}{16}$. Thus, the rate becomes

$R_{X}^{\phi,\,\xi}\;=\;\frac{2\times\sigma_{\phi
X}^{\xi~{}2}}{16\pi}\frac{\rho_{\phi}^{2}}{m_{\phi}^{4}}\sqrt{1-\frac{m_{X}^{2}}{m_{\phi}^{2}}}\,,$
(27)

where $\rho_{\phi}$ is the energy density of the inflaton and the interaction
term $\sigma_{\phi X}^{\xi}$ is given in Appendix A by Eq. (75).

It was shown in MO  that the dark matter production rate through the exchange
of a graviton, computed from the partial amplitude (18), is

$R^{\phi}_{X}=\frac{2\times\rho_{\phi}^{2}}{256\pi
M_{P}^{4}}\left(1+\frac{m_{X}^{2}}{2m^{2}_{\phi}}\right)^{2}\sqrt{1-\frac{m_{X}^{2}}{m_{\phi}^{2}}}\,,$
(28)

which can be written in the same form as (27) by defining an effective
coupling $\sigma_{\phi X}$

$\sigma_{\phi
X}\;=\;-\frac{m_{\phi}^{2}}{4M_{P}^{2}}\left(1+\frac{m_{X}^{2}}{2m_{\phi}^{2}}\right)\,.$
(29)

A comparison of the non-minimal to minimal rates for the production of dark
matter from inflaton scattering is shown in Fig. 4.

Figure 4: Contours of the ratio of the dark matter production rates from
oscillations in the inflaton condensate based on non-minimal gravitational
interactions to those based on minimal interactions. The ratio is displayed in
the $(\xi_{\phi},\xi_{X})$ plane.

For the production of Higgs bosons through inflaton condensate scattering, we
follow a similar procedure, and from the Lagrangian (13) we find

$R_{h}^{\phi,\,\xi}\;\simeq\;N_{h}\frac{2\times\sigma_{\phi
h}^{\xi~{}2}}{16\pi}\frac{\rho_{\phi}^{2}}{m_{\phi}^{4}}\,,$ (30)

where we assumed that $m_{h}\ll m_{\phi}$, $N_{h}=4$ is the number of internal
degrees of freedom for 1 complex Higgs doublet, and $\sigma_{\phi h}^{\xi}$ is
given in Appendix A by Eq. (76).

On the other hand, it was argued in CMOV  that the scattering
$\phi\phi\rightarrow hh$ through the graviton exchange can also be
parameterized by an effective coupling

$\mathcal{L}_{h}=-\sigma_{\phi h}\phi^{2}h^{2}\,,$ (31)

with

$\sigma_{\phi h}\;=\;-\frac{m_{\phi}^{2}}{4M_{P}^{2}}\,,$ (32)

and the rate $R_{h}^{\phi}$ is given by the analogous expression to (30) with
$\sigma_{\phi h}^{\xi}$ replaced by $\sigma_{\phi h}$.

The full four-point coupling of course is given by the sum $\sigma_{\phi
h/X}^{\xi}+\sigma_{\phi h/X}$. However, except for values where the two are
similar, which occurs when $12\xi^{2}+5\xi\simeq\frac{1}{2}$ (assuming
$m_{X}\ll m_{\phi}$ and taking the $\xi_{i}$ to be equal to $\xi$), either the
minimal or the non-minimal contribution dominates. Thus, for the most part, we
will consider separately the minimal and non-minimal contributions. Note that
for two values of $\xi$ ($\xi\sim-1/2$ and 1/12) destructive interference
could occur causing the entire rate to vanish (at the tree level).

##  III Particle Production with a Non-Minimal Coupling

Given the rates $R_{i}^{j}$ calculated in the previous section, we compute the
evolution for the gravitational (minimal and non-minimal) contribution to the
reheating processes and the dark matter density for the three reactions
outlined above.

###  III.1 $h~{}h\rightarrow X~{}X$

The gravitational scattering of thermal Higgs bosons leads to the production
of massive scalar dark matter particles $X$. The dark matter number density
$n_{X}$ can be calculated from the classical Boltzmann equation

$\frac{dn_{X}}{dt}+3Hn_{X}=R^{T}_{X}\,,$ (33)

where $H=\frac{\dot{a}}{a}$ is the Hubble parameter and the right-hand side of
the equation represents the dark matter production rate. It is more practical
to rewrite the above equation in terms of the scale factor $a$ rather than the
parameters $t$ or $T$.

We proceed by introducing the comoving number density $Y_{X}=na^{3}$ and
rewriting the Boltzmann equation as

$\frac{dY_{X}}{da}=\frac{a^{2}R^{T}_{X}(a)}{H(a)}\,.$ (34)

Since the production rate (23) is a function of the temperature of the thermal
bath, it is necessary to determine the relation between $T$ and $a$ in order
to solve the Boltzmann equation as a function of the scale factor $a$. For the
Starobinsky potential in Eq. (5), at the end of inflation, the inflaton starts
oscillating about a quadratic minimum, and we find the following energy
conservation equations777For the inflaton scattering with
$V(\phi)\sim\phi^{k}$, where $k>2$, see Bernal:2019mhf ; GKMO1 ; GKMO2 ; GKMOV
; Haque:2021mab ; Haque:2022kez ; Ahmed:2021fvt .

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

where $\rho_{\phi}$ and $\rho_{R}$ are the energy density of the inflaton and
radiation, respectively, $\Gamma_{\phi}$ is the inflaton decay rate, and for a
quadratic minimum, we are able to set the equation of state parameter
$w_{\phi}=\frac{P_{\phi}}{\rho_{\phi}}\simeq 0$. We will assume that reheating
occurs due to an effective inflaton coupling to the Standard Model fermions,
given by the interaction Lagrangian

$\mathcal{L}_{\phi-SM}^{y}\;=\;-y\phi\bar{f}f\,,$ (37)

where $y$ is a Yukawa-like coupling, $f$ is a Standard Model fermion, and the
inflaton decay rate is

$\Gamma_{\phi}\;=\;\frac{y^{2}}{8\pi}m_{\phi}\,.$ (38)

If we solve the Friedmann equations (35, 36), we find GKMO1 ; GKMO2 ; CMOV

$\rho_{\phi}(a)=\rho_{\rm end}\left(\frac{a_{\rm end}}{a}\right)^{3}$ (39)

and

$\rho_{R}(a)=\rho_{\rm RH}\left(\frac{a_{\rm
RH}}{a}\right)^{\frac{3}{2}}\frac{1-\left(\frac{a_{\rm
end}}{a}\right)^{\frac{5}{2}}}{1-\left(\frac{a_{\rm end}}{a_{\rm
RH}}\right)^{\frac{5}{2}}}\,,$ (40)

where $a_{\rm{end}}$ is the scale factor at the end of inflation, $\rho_{\rm
end}\equiv\rho_{\phi}(a_{\rm{end}})$ is the inflaton energy density at the end
of inflation when there is no radiation present, $a_{\rm{RH}}$ is the scale
factor at reheating, and $\rho_{\rm
RH}\equiv\rho_{R}(a_{\rm{RH}})=\rho_{\phi}(a_{\rm{RH}})$ is the energy density
at reheating. We note that these equations are strictly valid for $a_{\rm
end}\ll a\ll a_{\rm{RH}}$ and the end of inflation occurs when $\ddot{a}=0$
which corresponds to $\rho_{\rm end}=\frac{3}{2}V(\phi_{\rm{end}})$. For the
Starobinsky potential, $\rho_{\rm end}\simeq 0.175m_{\phi}^{2}M_{P}^{2}$ egno5
.

The radiation energy density can be parameterized as

$\rho_{R}=\frac{g_{T}\pi^{2}}{30}T^{4}\equiv\alpha T^{4}\,,$ (41)

where $g_{T}$ is the number of relativistic degrees of freedom at the
temperature $T$. The maximum temperature is attained when the radiation energy
density reaches its peak at $\rho_{R}(a_{\rm{max}})=\alpha T_{\rm{max}}^{4}$.
It was shown in GKMO1  that the ratio of $a_{\rm max}$ to $a_{\rm end}$ is
given by

$\frac{a_{\rm max}}{a_{\rm end}}=\left(\frac{8}{3}\right)^{\frac{2}{5}}\simeq
1.48\,.$ (42)

Using Eq. (40) we can then express the production rate from gravitational
scattering of thermal particles (23) as a function of the scale factor $a$

$R_{X}^{T,\,(\xi)}(a)\simeq\beta_{1}^{(\xi)}\frac{\rho_{\rm
RH}^{2}}{\alpha^{2}M_{P}^{4}}\left(\frac{a_{\rm
RH}}{a}\right)^{3}\left[\frac{1-\left(\frac{a_{\rm
end}}{a}\right)^{\frac{5}{2}}}{1-\left(\frac{a_{\rm end}}{a_{\rm
RH}}\right)^{\frac{5}{2}}}\right]^{2}\,,$ (43)

where we assumed that $m_{X}\ll m_{\phi},\,T$, and thus neglected the terms
$\beta_{2,\,3}^{(\xi)}$. If we use
$H\simeq\frac{\sqrt{\rho_{\phi}(a)}}{\sqrt{3}M_{P}}$, which is valid for $a\ll
a_{\rm RH}$, we can rewrite Eq. (34) as

$\frac{dY_{X}^{\xi}}{da}=\frac{\sqrt{3}M_{P}}{\sqrt{\rho_{\rm
RH}}}a^{2}\left(\frac{a}{a_{\rm
RH}}\right)^{\frac{3}{2}}R_{X}^{T,\,(\xi)}(a)\,.$ (44)

We find that the solution to this equation is

$\displaystyle~{}~{}n^{T,\,\xi}_{X}(a_{\rm
RH})=\frac{2\beta_{1}^{\xi}}{\sqrt{3}\alpha^{2}M_{P}^{3}}\frac{\rho_{\rm
RH}^{3/2}}{(1-(a_{\rm end}/a_{\rm RH})^{\frac{5}{2}})^{2}}\times$
$\displaystyle\left(1+3\left(\frac{a_{\rm end}}{a_{\rm
RH}}\right)^{\frac{5}{2}}-\frac{25}{7}\left(\frac{a_{\rm end}}{a_{\rm
RH}}\right)^{\frac{3}{2}}-\frac{3}{7}\left(\frac{a_{\rm end}}{a_{\rm
RH}}\right)^{5}\right)\,,$ (45)

where we integrated Eq. (44) in the interval $a_{\rm end}<a<a_{\rm RH}$.

The relic abundance is given by book

$\Omega_{X}h^{2}=1.6\times 10^{8}\frac{g_{0}}{g_{\rm RH}}\frac{n(T_{\rm
RH})}{T_{\rm RH}^{3}}\frac{m_{X}}{1~{}{\rm GeV}}\,,$ (46)

and if we combine it with Eq. (III.1), we obtain

$\displaystyle\Omega_{X}^{T,\,(\xi)}h^{2}\;=\;\frac{2}{3}\Omega_{k}^{(\xi)}$
$\displaystyle\left[1+3\left(\frac{\rho_{\rm{RH}}}{\rho_{\rm{end}}}\right)^{\frac{5}{6}}-\frac{25}{7}\left(\frac{\rho_{\rm{RH}}}{\rho_{\rm{end}}}\right)^{\frac{1}{2}}\right.$
(47)
$\displaystyle-\left.\frac{3}{7}\left(\frac{\rho_{\rm{RH}}}{\rho_{\rm{end}}}\right)^{\frac{5}{3}}\right]\,,$

with

$\Omega_{k}^{(\xi)}=1.6\times 10^{8}\frac{g_{0}}{g_{\rm
RH}}\frac{\beta_{1}^{(\xi)}\sqrt{3}}{\sqrt{\alpha}}\frac{m_{X}}{1~{}\rm{GeV}}\frac{T_{\rm
RH}^{3}}{M_{P}^{3}}\left[1-\left(\frac{\rho_{\rm RH}}{\rho_{\rm
end}}\right)^{\frac{5}{6}}\right]^{-2},$ (48)

where $g_{0}=43/11$ and we use the Standard Model value $g_{\rm RH}=427/4$.

We observe that $\Omega_{X}^{T,\,\xi}\propto\beta_{1}^{\xi}\,T_{\rm{RH}}^{3}$.
Therefore large values of the couplings $\xi_{h}$ and $\xi_{X}$ would require
a decrease in the reheating temperature. In Section IV we compare the
scattering rates and the dark matter abundances with the minimally coupled
case.

###  III.2 $\phi~{}\phi\rightarrow X~{}X$

Another mode of dark matter production is through the scattering of the
inflaton itself. Whereas the graviton exchange channel was treated with care
in MO ; CMOV , in the case of non-minimal coupling it suffices to replace
$R_{X}^{T,\,\xi}$ in Eq. (44) with the production rate (27),

$\frac{dY_{X}^{\xi}}{da}=\frac{\sqrt{3}M_{P}}{\sqrt{\rho_{\rm
RH}}}a^{2}\left(\frac{a}{a_{\rm
RH}}\right)^{\frac{3}{2}}R_{X}^{\phi,\,\xi}(a)\,,$ (49)

and to integrate between $a_{\rm end}$ and $a_{\rm RH}$, which leads to

$n_{X}^{\phi,\,\xi}(a_{\rm RH})=\frac{\sigma_{\phi
X}^{\xi~{}2}\rho_{\rm{RH}}^{3/2}M_{P}}{4\sqrt{3}\pi
m_{\phi}^{4}}\left[\left(\frac{a_{\rm{RH}}}{a_{\rm{end}}}\right)^{\frac{3}{2}}-1\right]\sqrt{1-\frac{m_{X}^{2}}{m_{\phi}^{2}}}\,.$
(50)

For $a_{\rm RH}\gg a_{\rm end}$, using Eq. (39) we can express
$n_{X}^{\phi,\,\xi}$ as a function of $\rho_{\rm end}$:

$n_{X}^{\phi,\,\xi}(a_{\rm RH})\simeq\frac{\sigma_{\phi X}^{\xi~{}2}\rho_{\rm
RH}\sqrt{\rho_{\rm end}}M_{P}}{4\sqrt{3}\pi
m_{\phi}^{4}}\sqrt{1-\frac{m_{X}^{2}}{m_{\phi}^{2}}}\,,$ (51)

and we find

$\displaystyle\frac{\Omega_{X}^{\phi,\,\xi}h^{2}}{0.12}\simeq\frac{1.3\times
10^{7}\sigma_{\phi
X}^{\xi~{}2}\rho_{\rm{RH}}^{1/4}M_{P}^{2}}{m_{\phi}^{3}}\frac{m_{X}}{1\,\rm{GeV}}\sqrt{1-\frac{m_{X}^{2}}{m_{\phi}^{2}}}\,,$
(52)

where we assumed the Starobinsky value for $\rho_{\rm end}$. The analogous
expression for models with minimally coupled gravity is found by replacing
$\sigma_{\phi X}^{\xi}\rightarrow\sigma_{\phi X}$.

Up to this point we have assumed that the radiation is produced via the direct
inflaton decay to a fermion pair. In the next subsection we discuss an
unavoidable radiation production channel when the inflaton condensate
scattering produces Higgs bosons in models with minimal and non-minimal
coupling to gravity.

###  III.3 $\phi~{}\phi\rightarrow h~{}h$

Gravitational processes that produce dark matter can also populate the thermal
bath in the same way. Even if this Planck-suppressed production mechanism does
not dominate throughout the entire reheating process, it was shown in CMOV
that for $T_{\rm RH}\lesssim 10^{9}$ GeV it is graviton exchange that
dominates the production of the thermal bath at the very beginning of the
reheating, when $\rho_{\phi}\sim\rho_{\rm end}$. In fact, it was shown that
the maximal temperature reached, $T_{\rm max}$, (which can be considered as an
absolute lower bound on $T_{\rm max}$) is $T_{\rm max}\sim 10^{12}$ GeV. It is
therefore natural to determine the value of the couplings ($\xi_{\phi}$,
$\xi_{h}$), for which non-minimal gravitational processes generate the thermal
bath at early times, and the maximal temperature which can be attained by
these processes.

Following the discussion in the previous subsection, to compute the radiation
energy density produced by gravitational couplings we implement the rate
$R_{h}^{\phi,\,\xi}$ (30) into the Friedmann equation (36)

$\frac{d\rho_{R}}{dt}+4H\rho_{R}\;\simeq\;N_{h}\frac{\sigma_{\phi
h}^{\xi~{}2}}{8\pi}\frac{\rho_{\phi}^{2}}{m_{\phi}^{3}}\,,$ (53)

where we took into account that each scattering corresponds to an energy
transfer of $2m_{\phi}$.888Or equivalently that each Higgs quanta carries an
energy $m_{\phi}$. The solution to this equation is

$\rho_{R}\;=\;N_{h}\frac{\sqrt{3}\sigma_{\phi
h}^{\xi~{}2}}{4\pi}\frac{\rho_{\rm
end}^{3/2}M_{P}}{m_{\phi}^{3}}\left[\left(\frac{a_{\rm{end}}}{a}\right)^{4}-\left(\frac{a_{\rm{end}}}{a}\right)^{\frac{9}{2}}\right]\,.$
(54)

Note that the dependence on the scale factor $a$ is very different from that
found in Eq. (40) due to inflaton decay. Indeed, the Higgs bosons produced by
gravitational scattering (minimal as well as non-minimal) are redshifted to a
greater extent because of the high dependence of the rate on their energy due
to the form of the energy-momentum tensor $T_{\mu\nu}^{0}$. Since
$\rho_{R}\propto a^{-4}$ in Eq. (54) (at large $a$) and $\rho_{\phi}\propto
a^{-3}$ in Eq. (39), reheating through this process does not occur (i.e.,
$\rho_{R}$ never comes to dominate the total energy at late times) and
inflaton decay is necessary.999This conclusion is avoided if the inflaton
potential about minimum is approximated by $\phi^{k}$ with a higher power of
$k>4$ CMOV ; Haque:2022kez .

However, as in the case of the reheating from the inflaton decay, the energy
density in Eq. (54) exhibits a maximum when $a=a_{\rm max}=(81/64)a_{\rm
end}$. The maximum radiation density is then,

$\rho_{\rm{max}}^{\xi}\;\simeq\;N_{h}\frac{\sigma_{\phi
h}^{\xi~{}2}}{12\sqrt{3}\pi}\frac{\rho_{\rm
end}^{3/2}M_{P}}{m_{\phi}^{3}}\left(\frac{8}{9}\right)^{8}\,,$ (55)

and from this expression we find that the maximum temperature produced by
gravitational interactions is given by

$\displaystyle T_{\rm{max}}^{\xi}$ $\displaystyle\simeq$ $\displaystyle
6.5\times 10^{11}\left(\frac{|\sigma_{\phi
h}^{\xi}|}{10^{-11}}\right)^{\frac{1}{2}}{\rm GeV}$ $\displaystyle\simeq$
$\displaystyle 1.8\times
10^{12}\sqrt{|\xi|}\left(|5+12\xi|\right)^{\frac{1}{2}}\left(\frac{m_{\phi}}{3\times
10^{13}\,{\rm GeV}}\right){\rm GeV}\,,$

where we took $\xi_{\phi}=\xi_{h}=\xi$ in the last equality. The analogous
expression for models with minimally coupled gravity is found by replacing
$\sigma_{\phi h}^{\xi}\rightarrow\sigma_{\phi h}$.

To compare the maximum temperature obtained by non-minimal interactions with
respect to minimal gravitational interactions, we can rewrite Eq. (III.3) (now
including minimal interactions in $T_{\rm max}^{\xi}$) as

$T_{\rm{max}}^{\xi}\simeq 1.3\times 10^{12}\left(\frac{|\sigma_{\phi
h}^{\xi}+\sigma_{\phi h}|}{\sigma_{\phi h}}\right)^{\frac{1}{2}}{\rm GeV}\,.$
(57)

The value of $\xi$ for which the maximum temperature generated by the non-
minimal coupling surpasses the one from graviton exchange is shown in Fig. 5
and is determined using

$\sqrt{\frac{|\sigma_{\phi h}^{\xi}|}{|\sigma_{\phi
h}|}}=\sqrt{2|\xi|}\left(|5+12\xi|\right)^{\frac{1}{2}}>1$ (58)

which is satisfied when $\xi>1/12$ or $\xi<-1/2$, as discussed earlier.

Figure 5: The maximum temperature during reheating generated separately by
minimal and non-minimal gravitational scattering of Higgs bosons in the
thermal bath.

As noted above and discussed in CMOV , minimal (and non-minimal) gravitational
interactions for a quadratic inflaton potential do not lead to the completion
of the reheating process, thus requiring additional inflaton interactions for
decay. Although radiation density produced in scattering falls off faster than
that from decay, at early time, the radiation density may in fact dominate and
determine $T_{\rm max}$. To determine when the $\phi~{}\phi\rightarrow h~{}h$
process leads to the maximum temperature, we rewrite Eq. (40) as:

$\rho^{y}_{R}\;=\;\frac{\sqrt{3}y^{2}m_{\phi}M_{P}^{3}}{20\pi}\left(\frac{\rho_{\rm
end}}{M_{P}^{4}}\right)^{\frac{1}{2}}\left[\left(\frac{a_{\rm
end}}{a}\right)^{\frac{3}{2}}-\left(\frac{a_{\rm
end}}{a}\right)^{4}\right]\,.$ (59)

Using Eq. (42), we find that the maximum radiation density produced by the
inflaton decay is given by

$\rho_{\rm
max}^{y}\;=\;\frac{\sqrt{3}y^{2}m_{\phi}M_{P}^{3}}{32\pi}\left(\frac{\rho_{\rm
end}}{M_{P}^{4}}\right)^{\frac{1}{2}}\left(\frac{3}{8}\right)^{\frac{3}{5}}\,.$
(60)

The maximum temperature is therefore determined by (non-minimal) gravitational
interactions when

$y^{2}\lesssim N_{h}\frac{8\rho_{\rm end}\sigma_{\phi
h}^{\xi~{}2}}{9m_{\phi}^{4}}\left(\frac{8}{9}\right)^{8}\left(\frac{8}{3}\right)^{\frac{3}{5}}$
(61)

or

$y\lesssim 1.6~{}\sigma_{\phi h}^{\xi}\sqrt{\frac{\rho_{\rm
end}}{m_{\phi}^{4}}}\simeq 5.4\times 10^{4}~{}\sigma_{\phi
h}^{\xi}\left(\frac{3\times 10^{13}\,\rm{GeV}}{m_{\phi}}\right)\,.$ (62)

This leads to the following reheating temperature:

$\displaystyle T_{\rm{RH}}\;$ $\displaystyle\lesssim$
$\displaystyle\;3.1\times 10^{19}\sigma_{\phi
h}^{\xi}\left(\frac{m_{\phi}}{3\times 10^{13}\,\rm{GeV}}\right)^{-1/2}{\rm
GeV}$ (63) $\displaystyle\lesssim$ $\displaystyle 2.4\times
10^{9}\left(\frac{m_{\phi}}{3\times
10^{13}}\right)^{\frac{3}{2}}\xi(5+12\xi)~{}{\rm GeV}$

where $T_{\rm RH}$ is given by GKMO2

$\rho_{\phi}(a_{\rm RH})=\alpha T_{\rm
RH}^{4}=\frac{12}{25}\Gamma_{\phi}^{2}M_{P}^{2}=\frac{3y^{4}m_{\phi}^{2}M_{P}^{2}}{400\pi^{2}}\,,$
(64)

when the reheating temperature is determined by inflaton decay.

The primary effect of the gravitational scattering processes on reheating is
the augmentation of $T_{\rm max}$ for sufficiently small inflaton decay
coupling, $y$. This can be seen in Fig. 6 where we show the evolution of the
energy density of radiation from scattering and decay as well as the energy
density of the inflaton as a function of $a/a_{\rm end}$ for $\sigma_{\phi
h}^{\xi}=0$ and $\sigma_{\phi h}^{\xi}/\sigma_{\phi h}=100$, respectively.

Figure 6: Evolution of the inflaton density (blue) and the total radiation
density (red), with radiation density produced from inflaton decays (dashed
orange) and $\phi~{}\phi\rightarrow h~{}h$ scattering processes
$\rho_{R}^{\sigma,\,\xi}$ (dotted green) and $\rho_{R}^{\sigma}$ (dash-dotted
purple) with $\sigma_{\phi h}^{\xi}/\sigma_{\phi h}=100$ (or
$\xi_{\phi}=\xi_{h}=\xi\simeq-2.3~{}{or}~{}1.8$), as a function of $a/a_{\rm
end}$ for a Yukawa-like coupling $y=10^{-8}$ and $\rho_{\rm end}\simeq
0.175\,m_{\phi}^{2}M_{P}^{2}$ $\simeq 9\times 10^{62}\,\rm{GeV}^{4}$. The
black dashed lines corresponds to the ratios $a_{\rm{int}}/a_{\rm{end}}\simeq
150$ and $6500$, which agrees with Eq. (66). The numerical solutions are
obtained from Eqs. (35), (36), and (53).

As we saw in Eq. (58), minimal gravitational interactions dominate over non-
minimal interactions when $\sigma_{\phi h}^{\xi}<\sigma_{\phi h}$ or when

$12\xi_{\phi}\xi_{h}+3\xi_{h}+2\xi_{\phi}<\frac{1}{2}\,,$ (65)

when we neglect contributions proportional to the Higgs mass. In this case,
the maximum temperature is determined by gravitational interactions when
$y\lesssim 2.1\times 10^{-6}$ from Eq. (62) using $\sigma_{\phi h}$ from Eq.
(32). The evolution of the energy densities in this case is shown in Fig. 6
with $y=10^{-8}$. However as the energy density of radiation after the maximum
falls faster than $\rho_{\phi}$, reheating in the Universe is determined by
the inflaton decay. For a sufficiently small coupling $y$, the energy density
from the decay dominates the radiation density at $a>a_{\rm int}$, where

$\frac{a_{\rm int}}{a_{\rm end}}\simeq\left(\frac{5\sigma_{\phi
h}^{2}N_{h}\rho_{\rm end}}{y^{2}m_{\phi}^{4}}\right)^{2/5}\simeq
1.6\left(\frac{\sigma_{\phi h}M_{P}}{ym_{\phi}}\right)^{4/5}\,.$ (66)

For $\sigma_{\phi h}=3.8\times 10^{-11}$, $m_{\phi}=3\times 10^{13}$ GeV, and
$y=10^{-8}$ we have $a_{\rm int}\approx 160a_{\rm end}$, as seen in the
figure.

When Eq. (65) is not satisfied, non-minimal interactions may dominate as shown
in the bottom panel of Fig. 6, for $\sigma_{\phi h}^{\xi}=100\sigma_{\phi h}$
and $y=10^{-8}$. The cross-over can be determined from Eq. (66) with the
replacement $\sigma_{\phi h}\rightarrow\sigma_{\phi h}^{\xi}$. In this
example, $a_{\rm int}\approx 6500a_{\rm end}$.

##  IV Results

We now turn to some general results that may be obtained from the framework
described above. Concerning the gravitational production of dark matter from
the thermal bath, the difficulty of populating the Universe via the exchange
of a graviton was already known Bernal:2018qlk ; CMOV . Summing the minimal
and non-minimal contributions in Eq. (47), we find for $\rho_{\rm
RH}\ll\rho_{\rm end}$

$\displaystyle\frac{\Omega_{X}^{T}}{0.12}$ $\displaystyle\simeq$
$\displaystyle\left[1+30f(\xi_{h},\xi_{X})\right]\left(\frac{T_{\rm
RH}}{10^{14}~{}{\rm GeV}}\right)^{3}\left(\frac{m_{X}}{4.0\times 10^{9}~{}{\rm
GeV}}\right)$ (67) $\displaystyle=$
$\displaystyle\left[1+120\xi^{2}(1+6\xi+12\xi^{2})\right]$
$\displaystyle\times\left(\frac{T_{\rm RH}}{10^{14}~{}{\rm
GeV}}\right)^{3}\left(\frac{m_{X}}{4.0\times 10^{9}~{}{\rm GeV}}\right)$

with

$f(\xi_{h},\xi_{X})=\xi_{h}^{2}+2\xi_{h}\xi_{X}+\xi_{X}^{2}+12\xi_{h}\xi_{X}\left(\xi_{h}+\xi_{X}+4\xi_{h}\xi_{X}\right)$

where we assumed $\xi_{h}=\xi_{X}=\xi$ in the last equality, for simplicity.
It is clear that, if we set $\xi=0$, i.e. if we consider only graviton
exchange, the reheating temperature necessary to obtain a reasonable density
respecting the data Planck  is dangerously close to the mass of the inflaton,
even for extremely large dark matter masses. This problem had already been
raised in Bernal:2018qlk  and resolved in MO ; CMOV  by considering the dark
matter produced from the (minimal) gravitational inflaton scattering.

On the other hand, from Eq. (67) we see that there is another solution to this
tension if one allows for non-minimal gravitational couplings. Indeed, it is
easy to see that for values of $\xi_{i}\gtrsim 0.1$
($f(\xi_{h},\xi_{X})\gtrsim{\frac{1}{30}}$), non-minimal gravitational
production dominates over graviton exchange. In this case, it becomes easier
to obtain the correct dark matter density for more reasonable values of
$T_{\rm RH}$ and/or $m_{X}$. For example, for a common value
$\xi=\xi_{h}=\xi_{X}=1$, a temperature of $T_{\rm RH}\sim 1.2\times 10^{13}$
GeV, thus slightly below the inflaton mass, is sufficient to produce an EeV
dark matter candidate, whereas for $\xi=1000$, $T_{\rm RH}\sim 10^{11}$ GeV
will saturate the relic density for a 2.6 TeV dark matter mass. We show this
result in Fig. 7 where we plot the reheating temperature needed to satisfy the
relic density constraint as function of $m_{X}$ for different value of $\xi$.
For each value of $\xi$, the relic density exceeds $\Omega_{X}h^{2}=0.12$
above the corresponding curve. As one can see, the line for $\xi=0$ is in the
upper corner of the figure at high values of $T_{\rm RH}$ and $m_{X}$ and
these drop significantly at higher values of $\xi$.

Figure 7: Region of parameter space respecting the relic density constraint
$\Omega_{X}h^{2}=0.12$ in the plane ($m_{X}$,$T_{\rm RH}$) for different
values of $\xi=\xi_{h}=\xi_{X}$ and $\rho_{\rm end}\simeq
0.175\,m_{\phi}^{2}M_{P}^{2}$ in the case of gravitational production from the
thermal bath $h~{}h\rightarrow X~{}X$. Both minimal and non-minimal
contributions are taken into account.

As was shown in MO ; CMOV , another possibility to avoid the necessity of high
reheating temperatures and/or dark matter masses is the production of matter
from the oscillations within the inflaton condensate when the energy stored in
the condensate is much larger than the reheating temperature. A simple
comparison between Eqs. (47) and (52) shows that the production of dark matter
via inflaton scattering when $\xi_{i}\neq 0$ generally dominates over the
production of dark matter from the thermal bath:

$\displaystyle\frac{\Omega_{X}^{\phi,\,\xi}}{\Omega_{X}^{T,\,\xi}}$
$\displaystyle\simeq$ $\displaystyle 34\frac{(\sigma_{\phi
X}^{\xi})^{2}}{\beta_{1}^{\xi}}\frac{M_{P}^{5}}{T_{\rm RH}^{2}m_{\phi}^{3}}$
(68) $\displaystyle\simeq$ $\displaystyle 185\frac{M_{P}m_{\phi}}{T_{\rm
RH}^{2}}\frac{(5+12\xi)^{2}}{1+6\xi+12\xi^{2}}\gg 1\,,\ $

where we took $\xi=\xi_{\phi}=\xi_{h}=\xi_{X}$ and $m_{X}\ll m_{\phi}$ in the
last equality. We are therefore able to state that the relic density of dark
matter generated by the non-minimal gravitational scattering of the inflaton
is always much more abundant than that produced by the thermal bath.

Dark matter production from inflaton scattering via minimal graviton exchange
also dominates over minimal gravitational thermal production CMOV . This state
of affairs is anything but surprising. Indeed, the energy available in the
inflaton condensate at the onset of oscillations is much greater than that
available in the thermal bath during the reheating process. As the scattering
cross-sections are themselves highly dependent on the energies through the
energy-momentum tensor, it is quite normal that inflaton scattering is the
dominant process for both minimal and non-minimal gravitational couplings.

Since inflaton scattering dominates in both the minimal and non-minimal
gravitational interactions we can compare the two. We obtain

$\frac{\Omega_{X}^{\phi,\,\xi}}{\Omega_{X}^{\phi}}=\frac{\sigma_{\phi
X}^{\xi~{}2}}{\sigma_{\phi X}^{2}}\simeq 4\xi^{2}(5+12\xi)^{2}\,,$ (69)

and we see again that non-minimal interactions dominate when $\xi>1/12$ or
$<-1/2$.

We show in Fig. 8 the region of the parameter space in the ($m_{X}$, $T_{\rm
RH}$) plane allowed by the relic density constraint, adding all of the minimal
and non minimal gravitational contributions, from inflaton scattering and as
well as Higgs scattering from the thermal bath taking
$\xi_{\phi}=\xi_{h}=\xi_{X}=\xi$. As expected, for $\xi=0$ we recover the
result found in CMOV . As one can see, the difficulty in the gravitational
production from the thermal bath is indeed alleviated as a reheating
temperature $T_{\rm RH}\simeq 10^{11}$ GeV allows for the production of a PeV
scale dark matter candidate. If in addition we introduce the non-minimal
couplings $\xi$, the necessary reheating temperature to fit the Planck data
may be as low as the electroweak scale for a GeV candidate if $\xi\gtrsim
1000$.

Figure 8: Region of parameter space respecting the relic density constraint
$\Omega_{X}h^{2}=0.12$ in the plane ($m_{X}$,$T_{\rm RH}$) for different
values of $\xi_{\phi}=\xi_{h}=\xi_{X}=\xi$ and $\rho_{\rm end}\simeq
0.175\,m_{\phi}^{2}M_{P}^{2}$ in the case of production from gravitational
inflaton scattering $\phi~{}\phi\rightarrow X~{}X$. Both minimal and non-
minimal contributions are taken into account.

Finally, we note that given the dark matter mass and reheating temperature (if
that sector of beyond the Standard Model physics were known), the contours in
Fig. 8 allow us to place an upper bound on the non-minimal couplings, $\xi$.
We can rewrite Eq. (52) as

$\displaystyle\frac{\Omega_{X}h^{2}}{0.12}$ $\displaystyle=$ $\displaystyle
4.1\times 10^{-7}(12\xi^{2}+5\xi+\frac{1}{2})^{2}\left(\frac{T_{\rm
RH}}{10^{10}{\rm GeV}}\right)$ (70)
$\displaystyle\times\left(\frac{m_{X}}{1{\rm
GeV}}\right)\left(\frac{m_{\phi}}{3\times 10^{13}{\rm GeV}}\right)\,,$

when $m_{X}\ll m_{\phi}$ and $\xi=\xi_{\phi}=\xi_{X}$. Then, for example, if
$m_{X}=1$ TeV, and $T_{\rm RH}=10^{9}$ GeV, we obtain an upper limit of
$|\xi|\lesssim 4$.

##  V Conclusions

In this paper, we have generalized the minimal gravitational interactions in
the early Universe, i.e., the s-channel exchange of a graviton, to include
non-minimal couplings of all scalars to the Ricci curvature $R$. We consider a
scalar sector $S_{i}$ consisting of the inflaton condensate $\phi$, the Higgs
field $H$ and a dark matter candidate $X$, and we have analyzed the impact of
couplings of the type $\xi_{i}S_{i}^{2}R$ on the reheating process and dark
matter production. The latter can be generated by the thermal Higgs scattering
or excitations of the inflaton, both through minimal and non-minimal
gravitational couplings. Whereas the Higgs scattering through the exchange of
a graviton necessitates a very large reheating temperature and/or dark matter
mass in order to fulfill Planck CMB constraints ($T_{\rm RH}\simeq 10^{14}$
GeV with $m_{X}\simeq 10^{9}$ GeV), for $\xi\gtrsim 0.1$, the non-minimal
coupling dominates the process and alleviates the tension. For $\xi\simeq
1000$, a dark matter mass of $\sim 1$ PeV with $T_{\rm RH}\simeq 10^{10}$ GeV
will satisfy the constraint, see Fig. 7. However, thermal production is not
the sole source of dark matter production through gravity. When we include the
contribution (necessarily present) of the inflaton scattering, we showed that
the energy stored in the condensate at the end of inflation compensates
largely the reduced gravitational Planck coupling. These processes yield the
correct relic abundance through minimal graviton exchange for a dark matter
mass of $\sim 10^{8}$ GeV with $T_{\rm RH}\simeq 10^{10}$ GeV, and the
constraint is satisfied for a dark matter mass of $\sim 100$ GeV and $T_{\rm
RH}\gtrsim 10^{4}$ GeV if one adds non-minimal couplings of the order
$\xi\simeq 100$ as we show in Fig. 8. Gravitational inflaton scattering also
affects the reheating process, producing a maximum temperature $\simeq
10^{12}$ GeV with minimal couplings, reaching as large as $T_{\rm
max}^{\xi}\simeq 5|\xi|T_{\rm max}\simeq 10^{14}$ GeV for $\xi=100$ as one can
see in Fig. 5. This result can be re-expressed as an upper limit to $|\xi|$
given values of $m_{X}$ and $T_{\rm RH}$.

We can not over-emphasize that all of our results are unavoidable, in the
sense that they are purely gravitational, and do not rely on physics beyond
the Strandard Mode. The relic density of dark matter, and maximum temperature
of the thermal bath computed here should be considered as lower bounds, that
should be implemented in any extension of the Standard Model, whatever is its
nature.

Note added : During the completion of the manuscript, some overlapping results
were presented in Aoki:2022dzd .

Acknowledgements.  The authors want to thank Emilian Dudas for useful
discussions. This work was made possible by with the support of the Institut
Pascal at Université Paris-Saclay during the Paris-Saclay Astroparticle
Symposium 2021, with the support of the P2IO Laboratory of Excellence (program
“Investisse ments d’avenir” ANR-11-IDEX-0003-01 Paris-Saclay and
ANR-10-LABX-0038), the P2I axis of the Graduate School Physics of Université
Paris-Saclay, as well as IJCLab, CEA, IPhT, APPEC, the IN2P3 master projet
UCMN and EuCAPT ANR-11-IDEX-0003-01 Paris-Saclay and ANR-10-LABX-0038). This
project has received support from the European Union’s Horizon 2020 research
and innovation programme under the Marie Sk  $-$ ${\rm l}$ odowska-Curie grant
agreement No 860881-HIDDeN. The work of K.A.O. and A.S. was supported in part
by DOE grant DE-SC0011842 at the University of Minnesota.

## Appendix

###  A PARTICLE PRODUCTION WITH A NON-MINIMAL COUPLING

The full Jordan frame action we consider is given by Eq. (2). The conformal
transformation to the Einstein frame is given by

$g_{\mu\nu}\;=\;\Omega^{2}\tilde{g}_{\mu\nu}\,,$ (71)

where $g_{\mu\nu}$ is the Einstein frame spacetime metric and the conformal
factor is expressed by Eq. (3). It can readily be shown that the scalar
curvature transforms as (see, e.g., Fujii:2003pa )

$\tilde{R}=\Omega^{2}\left[R+6g^{\mu\nu}\nabla_{\mu}\nabla_{\nu}\ln\Omega-6g^{\mu\nu}\left(\nabla_{\mu}\ln\Omega\right)\left(\nabla_{\nu}\ln\Omega\right)\right]\,.$
(72)

After eliminating the total divergence term, we find the Einstein frame action
(9).

To find the effective interaction terms we assume the small field limit (12)
and expand the conformal factors in the Einstein frame action. We find the
following effective interaction Lagrangian:

$\displaystyle\mathcal{L}_{\rm{eff}}\;$ $\displaystyle=$
$\displaystyle\;-\frac{1}{2}\left(\frac{\xi_{\phi}\phi^{2}}{M_{P}^{2}}+\frac{\xi_{X}X^{2}}{M_{P}^{2}}\right)\partial^{\mu}h\partial_{\mu}h-\frac{1}{2}\left(\frac{\xi_{h}h^{2}}{M_{P}^{2}}+\frac{\xi_{X}X^{2}}{M_{P}^{2}}\right)\partial^{\mu}\phi\partial_{\mu}\phi-\frac{1}{2}\left(\frac{\xi_{\phi}\phi^{2}}{M_{P}^{2}}+\frac{\xi_{h}h^{2}}{M_{P}^{2}}\right)\partial^{\mu}X\partial_{\mu}X$
(73) $\displaystyle+$
$\displaystyle\frac{6\xi_{h}\xi_{X}hX}{M_{P}^{2}}\partial^{\mu}h\partial_{\mu}X+\frac{6\xi_{h}\xi_{\phi}h\phi}{M_{P}^{2}}\partial^{\mu}h\partial_{\mu}\phi+\frac{6\xi_{\phi}\xi_{X}\phi
X}{M_{P}^{2}}\partial^{\mu}\phi\partial_{\mu}X+m_{X}^{2}X^{2}\left(\frac{\xi_{\phi}\phi^{2}}{M_{P}^{2}}+\frac{\xi_{h}h^{2}}{M_{P}^{2}}\right)$
$\displaystyle+$ $\displaystyle
m_{\phi}^{2}\phi^{2}M_{P}^{2}\left(\frac{\xi_{X}X^{2}}{M_{P}^{2}}+\frac{\xi_{h}h^{2}}{M_{P}^{2}}\right)+m_{h}^{2}h^{2}\left(\frac{\xi_{\phi}\phi^{2}}{M_{P}^{2}}+\frac{\xi_{X}X^{2}}{M_{P}^{2}}\right)\,,$

and we can rewrite the above Lagrangian in terms of the effective couplings as
Eq. (13), with

$\displaystyle\sigma_{hX}^{\xi}\;$ $\displaystyle=$
$\displaystyle\;\frac{1}{4M_{P}^{2}}\left[\xi_{h}(2m_{X}^{2}+s)+\xi_{X}(2m_{h}^{2}+s)\right.$
(74) $\displaystyle+$
$\displaystyle\left.\left(12\xi_{X}\xi_{h}(m_{h}^{2}+m_{X}^{2}-t)\right)\right]\,,$
$\sigma_{\phi
X}^{\xi}\;=\;\frac{1}{2M_{P}^{2}}\left[\xi_{\phi}m_{X}^{2}+12\xi_{\phi}\xi_{X}m_{\phi}^{2}+3\xi_{X}m_{\phi}^{2}+2\xi_{\phi}m_{\phi}^{2}\right]\,,$
(75) $\sigma_{\phi
h}^{\xi}\;=\;\frac{1}{2M_{P}^{2}}\left[\xi_{\phi}m_{h}^{2}+12\xi_{\phi}\xi_{h}m_{\phi}^{2}+3\xi_{h}m_{\phi}^{2}+2\xi_{\phi}m_{\phi}^{2}\right]\,,$
(76)

where $s,t$ are the Mandelstam variables. The latter couplings assume an
inflaton condensate in the initial state rather than a thermal Higgs in the
initial state accounting for the lack of symmetry in the three couplings.

###  B THERMAL PRODUCTION

In this appendix we calculate the thermal dark matter production rate
$R_{X}^{T,\,\xi}(T)$ arising from the effective four-point interaction
$\sigma_{hX}h^{2}X^{2}$, where $\sigma_{hX}$ is given by Eq. (74). We also
calculate the production rate $R_{X}^{T}(T)$ for the thermal scattering
processes mediated by gravity alone, ${\rm{SM}}~{}{\rm{SM}}\rightarrow X~{}X$,
that are unavoidable in models with a minimal coupling to gravity
($\xi_{\phi,h,X}=0$) Bernal:2018qlk ; CMOV , and compare the two results.

The production rate $R_{X}^{T,\,\xi}(T)$ can be computed from Eq. (21). The
matrix element squared is given by

$|\overline{{\cal M}}^{hX,\,\xi}|^{2}\;=\;4\sigma_{hX}^{\xi~{}2}\,,$ (77)

where in the limit where the Higgs boson mass is neglected, the Mandelstam
variables $s$ and $t$ are given by

$\displaystyle t\,=$
$\displaystyle\,\dfrac{s}{2}\left(\sqrt{1-\frac{4m_{X}^{2}}{s}}\cos\theta_{13}-1\right)+m_{X}^{2}\,,$
(78) $\displaystyle s\,=$
$\displaystyle\,2E_{1}E_{2}\left(1-\cos\theta_{12}\right)\,.$ (79)

We find the following coefficients for Eq. (23)

$\displaystyle\beta_{1}^{\xi}$
$\displaystyle\;=\;\frac{\pi^{3}}{2700}\left[\xi_{h}^{2}+2\xi_{h}\xi_{X}+\xi_{X}^{2}+12\xi_{h}\xi_{X}\left(\xi_{h}+\xi_{X}+4\xi_{h}\xi_{X}\right)\right]\,,$
(80) $\displaystyle\beta_{2}^{\xi}$
$\displaystyle\;=\;\frac{\zeta(3)^{2}\xi_{h}}{2\pi^{5}}\left[\xi_{h}+\xi_{X}+6\xi_{h}\xi_{X}-12\xi_{h}\xi_{X}^{2}\right]\,,$
(81) $\displaystyle\beta_{3}^{\xi}$
$\displaystyle\;=\;\frac{\xi_{h}^{2}}{576\pi}\,.$ (82)

Similarly, using Eqs. (18)-(20), we find the matrix element squared for
minimally coupled gravity:

$|\overline{{\cal
M}}^{hX}|^{2}\;=\;\frac{1}{4M_{P}^{4}}\frac{\left(t(s+t)-2m_{X}^{2}t+m_{X}^{4}\right)^{2}}{s^{2}},$
(83)

where we have neglected the Higgs field mass. We find the coefficients:

$\displaystyle\beta_{1}$ $\displaystyle\;=\;\frac{\pi^{3}}{81000}\,,$ (84)
$\displaystyle\beta_{2}$
$\displaystyle\;=\;-\frac{\zeta(3)^{2}}{30\pi^{5}}\,,$ (85)
$\displaystyle\beta_{3}$ $\displaystyle\;=\;\frac{1}{4320\pi}\,.$ (86)

Note that when both contributions are kept, and we neglect $m_{h}\ll m_{X}$,
the full coefficients (including interference) are given by

$\displaystyle\beta_{1}^{\xi}$ $\displaystyle=$
$\displaystyle\frac{\pi^{3}}{81000}\left[30\xi_{h}^{2}\left(12\xi_{X}(4\xi_{X}+1)+1\right)\right.$
(87)
$\displaystyle\left.+10\xi_{h}(6\xi_{X}+1)^{2}+10\xi_{X}(3\xi_{X}+1)+1\right]\,,$
$\displaystyle\beta_{2}^{\xi}$ $\displaystyle=$
$\displaystyle-\frac{\zeta(3)^{2}}{60\pi^{5}}\left[2+10\xi_{X}\right.$ (88)
$\displaystyle\left.+5\xi_{h}\left(1+6\xi_{X}+6\xi_{h}\left(6\xi_{X}(2\xi_{X}-1)-1\right)\right)\right]\,,$
$\displaystyle\beta_{3}^{\xi}$ $\displaystyle=$
$\displaystyle\frac{1}{8640\pi}\left[2+5\xi_{h}\left(32\xi_{h}-2\right)\right]\,.$
(89)

which reduces to Eqs. (84-86) when all $\xi_{i}=0$.

## References

  * (1) F. L. Bezrukov and M. Shaposhnikov, Phys. Lett. B 659 (2008), 703-706 [arXiv:0710.3755 [hep-th]].
  * (2) O. Lebedev, Prog. Part. Nucl. Phys. 120 (2021), 103881 doi:10.1016/j.ppnp.2021.103881 [arXiv:2104.03342 [hep-ph]].
  * (3) Y. Ema, R. Jinno, K. Mukaida and K. Nakayama, JCAP 02 (2017), 045 [arXiv:1609.05209 [hep-ph]].
  * (4) B. A. Bassett and S. Liberati, Phys. Rev. D 58, 021302 (1998) [erratum: Phys. Rev. D 60, 049902 (1999)] [arXiv:hep-ph/9709417 [hep-ph]]; S. Tsujikawa, K. i. Maeda and T. Torii, Phys. Rev. D 60, 063515 (1999) [arXiv:hep-ph/9901306 [hep-ph]]; T. Markkanen and S. Nurmi, JCAP 02, 008 (2017) [arXiv:1512.07288 [astro-ph.CO]].
  * (5) Y. Mambrini and K. A. Olive, Phys. Rev. D 103 (2021) no.11, 115009 [arXiv:2102.06214 [hep-ph]].
  * (6) S. Clery, Y. Mambrini, K. A. Olive and S. Verner, [arXiv:2112.15214 [hep-ph]];
  * (7) B. Barman and N. Bernal, JCAP 06 (2021), 011 doi:10.1088/1475-7516/2021/06/011 [arXiv:2104.10699 [hep-ph]]; B. Barman, N. Bernal, Y. Xu and Ó. Zapata, [arXiv:2202.12906 [hep-ph]].
  * (8) M. R. Haque and D. Maity, [arXiv:2112.14668 [hep-ph]].
  * (9) M. R. Haque and D. Maity, [arXiv:2201.02348 [hep-ph]].
  * (10) C. G. Callan, Jr., S. R. Coleman and R. Jackiw, Annals Phys. 59, 42-73 (1970) doi:10.1016/0003-4916(70)90394-5
  * (11) L. Parker, Phys. Rev. Lett. 21, 562-564 (1968); L. Ford, Phys. Rev. D 35, 2955 (1987).
  * (12) D. J. H. Chung, E. W. Kolb and A. Riotto, Phys. Rev. Lett. 81, 4048-4051 (1998) [arXiv:hep-ph/9805473 [hep-ph]]; S. Hashiba and J. Yokoyama, JCAP 01, 028 (2019) [arXiv:1809.05410 [gr-qc]]; N. Bernal, A. Chatterjee and A. Paul, JCAP 12, 020 (2018) [arXiv:1809.02338 [hep-ph]]; L. Li, T. Nakama, C. M. Sou, Y. Wang and S. Zhou, JHEP 07, 067 (2019) [arXiv:1903.08842 [astro-ph.CO]]; J. A. R. Cembranos, L. J. Garay and J. M. Sánchez Velázquez, JHEP 06, 084 (2020) [arXiv:1910.13937 [hep-ph]]; N. Herring, D. Boyanovsky and A. R. Zentner, Phys. Rev. D 101, no.8, 083516 (2020) [arXiv:1912.10859 [gr-qc]]; E. Babichev, D. Gorbunov, S. Ramazanov and L. Reverberi, JCAP 09 (2020), 059 [arXiv:2006.02225 [hep-ph]]; S. Ling and A. J. Long, Phys. Rev. D 103, no.10, 103532 (2021) [arXiv:2101.11621 [astro-ph.CO]].
  * (13) Y. Ema, R. Jinno, K. Mukaida and K. Nakayama, JCAP 05, 038 (2015) [arXiv:1502.02475 [hep-ph]]; Y. Ema, R. Jinno, K. Mukaida and K. Nakayama, Phys. Rev. D 94, no.6, 063517 (2016) [arXiv:1604.08898 [hep-ph]]; Y. Ema, K. Nakayama and Y. Tang, JHEP 09, 135 (2018) [arXiv:1804.07471 [hep-ph]]; Y. Ema, K. Mukaida and K. Nakayama, [arXiv:2112.12774 [hep-ph]].
  * (14) Y. Ema, K. Nakayama and Y. Tang, JHEP 07, 060 (2019) [arXiv:1903.10973 [hep-ph]]. N. Herring and D. Boyanovsky, Phys. Rev. D 101, no.12, 123522 (2020) [arXiv:2005.00391 [astro-ph.CO]];
  * (15) A. Ahmed, B. Grzadkowski and A. Socha, JHEP 08, 059 (2020) [arXiv:2005.01766 [hep-ph]]; C. Gross, S. Karamitsos, G. Landini and A. Strumia, JHEP 03, 174 (2021) [arXiv:2012.12087 [hep-ph]]; B. Barman, N. Bernal, A. Das and R. Roshan, JCAP 01, no.01, 047 (2022) [arXiv:2108.13447 [hep-ph]].
  * (16) A. A. Starobinsky, Phys. Lett. B 91, 99 (1980).
  * (17) M. Atkins and X. Calmet, Phys. Rev. Lett. 110, no.5, 051301 (2013) [arXiv:1211.0281 [hep-ph]]; Z. Z. Xianyu, J. Ren and H. J. He, Phys. Rev. D 88, 096013 (2013) [arXiv:1305.0251 [hep-ph]].
  * (18) R. A. Flores and M. Sher, Phys. Rev. D 27, 1679 (1983); M. Sher, Phys. Rept. 179, 273-418 (1989); G. Isidori, G. Ridolfi and A. Strumia, Nucl. Phys. B 609, 387-409 (2001) [arXiv:hep-ph/0104016 [hep-ph]]; F. Bezrukov, M. Y. Kalmykov, B. A. Kniehl and M. Shaposhnikov, JHEP 10, 140 (2012) [arXiv:1205.2893 [hep-ph]]; G. Degrassi, S. Di Vita, J. Elias-Miro, J. R. Espinosa, G. F. Giudice, G. Isidori and A. Strumia, JHEP 08, 098 (2012) [arXiv:1205.6497 [hep-ph]]; D. Buttazzo, G. Degrassi, P. P. Giardino, G. F. Giudice, F. Sala, A. Salvio and A. Strumia, JHEP 12, 089 (2013) [arXiv:1307.3536 [hep-ph]]; A. V. Bednyakov, B. A. Kniehl, A. F. Pikelner and O. L. Veretin, Phys. Rev. Lett. 115, no.20, 201802 (2015) [arXiv:1507.08833 [hep-ph]].
  * (19) A. Vilenkin, Phys. Lett. 115B, 91 (1982); A. Vilenkin and L. H. Ford, Phys. Rev. D 26, 1231 (1982); A. D. Linde, Phys. Lett. 116B, 335 (1982); A. A. Starobinsky, Phys. Lett. 117B, 175 (1982); A. Vilenkin, Nucl. Phys. B 226, 527 (1983); A. D. Linde, Phys. Lett. 160B, 243 (1985); K. Enqvist, K. W. Ng and K. A. Olive, Nucl. Phys. B 303, 713 (1988).
  * (20) J. R. Espinosa, G. F. Giudice and A. Riotto, JCAP 05, 002 (2008) [arXiv:0710.2484 [hep-ph]].
  * (21) M. Herranen, T. Markkanen, S. Nurmi and A. Rajantie, Phys. Rev. Lett. 113, no.21, 211102 (2014) [arXiv:1407.3141 [hep-ph]]. T. Markkanen, S. Nurmi, A. Rajantie and S. Stopyra, JHEP 06, 040 (2018) [arXiv:1804.02020 [hep-ph]].
  * (22) T. Markkanen, A. Rajantie and S. Stopyra, Front. Astron. Space Sci. 5, 40 (2018) [arXiv:1809.06923 [astro-ph.CO]].
  * (23) N. Aghanim et al. [Planck], Astron. Astrophys. 641, A6 (2020) [arXiv:1807.06209 [astro-ph.CO]].
  * (24) See e.g., J. Ellis, M. A. G. García, N. Nagata, D. V. Nanopoulos, K. A. Olive and S. Verner, Int. J. Mod. Phys. D 29, no.16, 2030011 (2020) [arXiv:2009.01709 [hep-ph]].
  * (25) F. Bezrukov, A. Magnin, M. Shaposhnikov and S. Sibiryakov, JHEP 01, 016 (2011) [arXiv:1008.5157 [hep-ph]].
  * (26) M. Garny, M. Sandora and M. S. Sloth, Phys. Rev. Lett. 116, no.10, 101302 (2016) [arXiv:1511.03278 [hep-ph]]; M. Garny, A. Palessandro, M. Sandora and M. S. Sloth, JCAP 02, 027 (2018) [arXiv:1709.09688 [hep-ph]].
  * (27) Y. Tang and Y. L. Wu, Phys. Lett. B 774, 676-681 (2017) [arXiv:1708.05138 [hep-ph]].
  * (28) M. Chianese, B. Fu and S. F. King, JCAP 06, 019 (2020) [arXiv:2003.07366 [hep-ph]]; M. Chianese, B. Fu and S. F. King, JCAP 01, 034 (2021) [arXiv:2009.01847 [hep-ph]].
  * (29) M. Redi, A. Tesi and H. Tillim, [arXiv:2011.10565 [hep-ph]].
  * (30) K. Benakli, Y. Chen, E. Dudas and Y. Mambrini, Phys. Rev. D 95 (2017) no.9, 095002 [arXiv:1701.06574 [hep-ph]]; E. Dudas, Y. Mambrini and K. Olive, Phys. Rev. Lett. 119 (2017) no.5, 051801 [arXiv:1704.03008 [hep-ph]].
  * (31) N. Bernal, M. Dutra, Y. Mambrini, K. Olive, M. Peloso and M. Pierre, Phys. Rev. D 97 (2018) no.11, 115020 [arXiv:1803.01866 [hep-ph]].
  * (32) M. A. G. Garcia, K. Kaneta, Y. Mambrini and K. A. Olive, JCAP 04, 012 (2021) [arXiv:2012.10756 [hep-ph]].
  * (33) K. Ichikawa, T. Suyama, T. Takahashi and M. Yamaguchi, Phys. Rev. D 78 (2008), 063545 [arXiv:0807.3988 [astro-ph]].
  * (34) K. Kainulainen, S. Nurmi, T. Tenkanen, K. Tuominen and V. Vaskonen, JCAP 06 (2016), 022 [arXiv:1601.07733 [astro-ph.CO]].
  * (35) N. Bernal, F. Elahi, C. Maldonado and J. Unwin, JCAP 11 (2019), 026 [arXiv:1909.07992 [hep-ph]].
  * (36) M. A. G. Garcia, K. Kaneta, Y. Mambrini and K. A. Olive, Phys. Rev. D 101 (2020) no.12, 123507 [arXiv:2004.08404 [hep-ph]].
  * (37) M. A. G. Garcia, K. Kaneta, Y. Mambrini, K. A. Olive and S. Verner, [arXiv:2109.13280 [hep-ph]].
  * (38) A. Ahmed, B. Grządkowski and A. Socha, [arXiv:2111.06065 [hep-ph]].
  * (39) J. Ellis, M. A. G. Garcia, D. V. Nanopoulos and K. A. Olive, JCAP 07, 050 (2015) [arXiv:1505.06986 [hep-ph]].
  * (40) Y. Mambrini, “Particles in the dark Universe", Springer Ed., ISBN 978-3-030-78139-2
  * (41) S. Aoki, H. M. Lee, A. G. Menkara and K. Yamashita, [arXiv:2202.13063 [hep-ph]].
  * (42) Y. Fujii and K. Maeda, “The scalar-tensor theory of gravitation,” Cambridge University Press, 2003.