aainstitutetext: New York University Abu Dhabi
PO Box 129188, Saadiyat Island, Abu Dhabi, United Arab
Emiratesbbinstitutetext: PRISMA+ Cluster of Excellence and Mainz Institute for
Theoretical Physics
Johannes Gutenberg University, 55099 Mainz, Germanyccinstitutetext: Physik
Department, TUM School of Natural Sciences, TU München,
James-Franck-Straße, D-85748 Garching, Germany

# Graviton- and Inflaton-mediated Dark Matter Production
after Large Field Polynomial Inflation

Nicolás Bernal b Julia Harz b,c Martin A. Mojahed b and Yong Xu
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Polynomial inflation is a simple cosmological scenario, which fits the cosmic
microwave background data well. It provides testable predictions for the
tensor-to-scalar ratio and the running of the spectral index. In this work, we
investigate the production of Dirac dark matter (DM) within the framework of
large-field polynomial inflation. We study all relevant production channels
including $i$) non-thermal production through inflaton decays and scatterings,
and $ii$) thermal production from scattering of standard model particles
mediated by inflatons and gravitons. In contrast to small-field polynomial
inflation, where inflaton decay dominates DM production, we find that
graviton-mediated processes can be dominant in the large-field scenario. In
particular, graviton-mediated scatterings can account for the observed relic
abundance even when DM is up to one order of magnitude heavier than the
inflaton. For DM lighter than the inflaton, we demonstrate that the interplay
between graviton- and inflaton-mediated production channels give rise to non-
trivial relations between the DM mass and the reheating temperature required
to account for the DM relic abundance.

MITP-24-056

##  1 Introduction

Cosmic inflation is an established paradigm for solving the horizon and
flatness problem of standard cosmology Starobinsky:1980te ; Guth:1980zm ;
Linde:1981mu ; Albrecht:1982wi . The simplest inflationary scenario is the so-
called slow roll (SR) single-field inflation, where a scalar inflaton field
$\phi$ slowly rolls down its potential Martin:2013tda . For the minimal SR
model with a monomial potential $V(\phi)\propto\phi^{p}$, the most recent
cosmic microwave background (CMB) experiments Planck:2018jri ; BICEP2:2018kqh
have ruled out scenarios with $p\geq 1/2$, since such potentials are too steep
and yield too large tensor-to-scalar ratios $r$. Moreover, according to the
latest Planck data Planck:2018jri , the most favored single-field inflation
models within the Einstein gravity framework are those with concave potentials
$V^{\prime\prime}(\phi)<0$.

As a simple extension of the monomial paradigm, one can consider a
renormalizable polynomial potential with a concave shape that contains a
quadratic, cubic, and a quartic term Hodges:1989dw ; Allahverdi:2006iq ;
Destri:2007pv ; Nakayama:2013jka ; Nakayama:2013txa ; Kallosh:2014xwa ;
Li:2014zfa ; Aslanyan:2015hmi ; Gao:2015yha ; Musoke:2017frr . The cosmic
inflationary dynamics of this model, as well as the parameter space that is in
agreement with CMB data, have recently been studied in the context of small-
and large-field excursions (Refs. Drees:2021wgd  and Drees:2022aea ,
respectively). After the end of inflation, the energy density contained in the
inflaton field has to be transferred to standard-model (SM) degrees of freedom
to generate the SM thermal plasma, for instance via inflaton decays into SM
particles. Interestingly, in the small-field realization of the model, the
reheating temperature is typically smaller than the inflaton mass, while in
the large-field realization it can be much higher.

It is well known that physics beyond the standard model (BSM) is required to
explain both the existence of dark matter (DM) and the dynamical generation of
the baryon asymmetry of the Universe. In the context of polynomial inflation,
a rather large parameter space was found for both thermal and non-thermal
baryogenesis via leptogenesis Xu:2022qpx ; Drees:2024hok . DM can be produced
by several processes, including non-thermal production from inflaton decays
and scatterings, and thermal production from inflaton- and graviton-mediated
scatterings of SM particles in the thermal plasma. In the context of small-
field polynomial inflation, the generation of fermionic DM was investigated in
Ref. Bernal:2021qrl  where it was shown that inflaton decays dominate over all
other production channels due to radiative constraints on the reheating
temperature $T_{\text{rh}}$, where the latter can be at most
$T_{\text{rh}}\sim\mathcal{O}(10^{11})$ GeV.

In this work, we extend the previous work Bernal:2021qrl  by considering a
large-field setup. In this case, radiative corrections to the inflaton
potential do not place strong constraints on reheating, leading to reheating
temperatures as high as $T_{\text{rh}}\sim\mathcal{O}(10^{14})$ GeV
Drees:2022aea . In contrast to the previously studied small-field case, we
find that DM production channels mediated by gravity or the inflaton can
become important. Here we present a detailed account of the interplay of DM
production via gravity-mediated processes and inflation decays and
scatterings. Although the current analysis is focused on polynomial inflation,
our formalism can be applied to other large-field inflation models such as
Starobinsky inflation Starobinsky:1980te , Higgs inflation Bezrukov:2007ep  or
attractor models Kallosh:2013yoa .

The outline of this paper is as follows. We introduce the model setup in
Section 2. The dynamics of inflation and (p)reheating are discussed in Section
3. It will set the stage for the main focus of this work, the production of
DM, which is investigated in Section 4. Here, we analyze all possible
production channels and the corresponding parameter space to successfully
explain the DM abundance. We summarize our findings in Section 5.

##  2 Model Setup

Inspired by Occam’s razor, and the necessity of new physics explaining the
primordial evolution of the Universe and DM, we extend the SM of particle
physics with two new states, which account for the inflaton and DM. The action
$S$ of the model can be written as follows Bernal:2021qrl

$S=\int d^{4}x\,\sqrt{-g}\left(\mathcal{L}_{\rm EH}+\mathcal{L}_{\rm
SM}+\mathcal{L}_{\phi}+\mathcal{L}_{\text{DM}}+\mathcal{L}_{H\phi}\right),$
(1)

where $g$ is the determinant of the Friedmann-Lemaître-Robertson-Walker metric
defined by $g_{\mu\nu}=\text{diag}(1,-a^{2},-a^{2},-a^{2})$, and $a$ is the
cosmic scale factor. $\mathcal{L}_{\rm EH}$ denotes the Einstein-Hilbert
Lagrangian density given by

$\mathcal{L}_{\rm EH}=\frac{M_{P}^{2}}{2}\,R\,,$ (2)

where $R$ denotes the Ricci scalar and $M_{P}\simeq 2.4\times 10^{18}$ GeV is
the reduced Planck mass. $\mathcal{L}_{\rm SM}$ is the usual SM Lagrangian and
$\mathcal{L}_{\phi}$ encodes the inflaton dynamics

$\mathcal{L}_{\phi}=\frac{1}{2}\partial_{\mu}\phi\partial^{\mu}\phi-V(\phi)\,.$
(3)

The inflaton $\phi$ is a real scalar field, and we consider the most general
renormalizable potential111A linear term can always be removed through a field
redefinition.

$V(\phi)=b\,\phi^{2}+c\,\phi^{3}+d\,\phi^{4}.$ (4)

Additionally, we assume that DM is a Dirac fermion $\chi$ with the Lagrangian
density222The results presented in this paper would be qualitatively similar
for the case of a DM Majorana fermion. See Ref. Izumine:2024hyp  for a recent
study of scalar DM production during reheating.

$\mathcal{L}_{\text{DM}}=-y_{\chi}\,\phi\,\overline{\chi}\,\chi+m_{\chi}\,\overline{\chi}\,\chi\,,$
(5)

where $y_{\chi}$ denotes the Yukawa coupling to the inflaton and $m_{\chi}$
denotes the DM mass. The reheating of the Universe is controlled by
interactions of $\phi$ with the SM Higgs doublet $H$

$\mathcal{L}_{H\phi}=-\mu\,\phi\,H^{\dagger}H-\frac{1}{2}\lambda_{\phi
H}\,\phi^{2}\,H^{\dagger}H\,.$ (6)

In this setup, gravitational interactions are also relevant. Expanding the
metric in Eq. (1) around a flat Minkowski background $\eta_{\mu\nu}$ as
$g_{\mu\nu}\simeq\eta_{\mu\nu}+(2/M_{P})\,h_{\mu\nu}$ gives rise to couplings
between the energy-momentum tensor $T_{i}^{\mu\nu}$ for all matter fields $i$
and the graviton $h_{\mu\nu}$ as Choi:1994ax

$\sqrt{-g}\,\mathcal{L}\supset\frac{1}{M_{P}}\,h_{\mu\nu}\,\sum_{i}T_{i}^{\mu\nu},$
(7)

with $i=\\{\phi,\,\chi,\,\text{SM}\\}$. These gravitational interactions can
play a crucial role in particle production in the early Universe, as we will
discuss later.

##  3 Inflation and (P)reheating

In the following, we briefly review the dynamics of polynomial inflation and
the cosmic (p)reheating taking place after inflation.

###  3.1 Polynomial Inflation

In polynomial inflation, the potential takes a form as given in Eq. (4). In
the case where $b=\frac{9}{32}\frac{c^{2}}{d}$, the potential has an
inflection point at $\phi=\phi_{0}\equiv-\frac{3}{8}\frac{c}{d}$, and can be
reparameterized as Drees:2021wgd ; Drees:2022aea

$V(\phi)=d\left[2\phi_{0}^{2}\,\phi^{2}-\frac{8}{3}(1-\beta)\,\phi_{0}\,\phi^{3}+\phi^{4}\right].$
(8)

The parameter $\beta$ controls the flatness of the potential near $\phi_{0}$.
For $\beta=0$, one has an exact inflection point at $\phi_{0}$. A false vacuum
is generated at $\phi>\phi_{0}$ if $\beta<0$, in which case the inflaton may
not roll down to the global minimum of the potential. The latter possibility,
which is physically disfavored as it would not lead to a hot Big Bang, is
avoided in this paper by only considering $0<\beta\ll 1$.

We apply the SR approximation and define the following SR parameters
Lyth:2009zz

$\epsilon_{V}\equiv\frac{M_{P}^{2}}{2}\left(\frac{V^{\prime}}{V}\right)^{2},\quad\eta_{V}\equiv
M_{P}^{2}\,\frac{V^{\prime\prime}}{V}\,,\quad\xi_{V}^{2}\equiv
M_{P}^{4}\,\frac{V^{\prime}\,V^{\prime\prime\prime}}{V^{2}}\,,$ (9)

which must be smaller than one during the accelerated expansion of the
Universe. Inflation ends when the field reaches the value
$\phi=\phi_{\text{end}}$ defined as $\epsilon_{V}(\phi_{\text{end}})=1$.333In
large field inflation, the potential is dominated by the quadratic term at the
end of inflation, which is independent of $\beta$. The total number of
$e$-folds $N_{\star}$ between the time when the CMB pivot scale
$k_{\star}=0.05$ Mpc-1 first crossed the horizon until the end of inflation is
given by Drees:2022aea

$\displaystyle N_{\star}$
$\displaystyle=\int^{\phi_{\star}}_{\phi_{\text{end}}}\frac{1}{\sqrt{2\,\epsilon_{V}(\phi)}}\,\frac{d\phi}{M_{P}}$
$\displaystyle\simeq\frac{1}{24}\Bigg{[}3\left(\frac{\phi}{M_{P}}\right)^{2}-4\,\frac{\phi\,\phi_{0}}{M_{P}^{2}}+15\left(\frac{\phi_{0}}{M_{P}}\right)^{2}$
$\displaystyle\qquad\quad-\left(\frac{\phi_{0}}{M_{P}}\right)^{2}\sqrt{\frac{2}{\beta}}\arctan\left(\frac{\phi_{0}-\phi}{\sqrt{2\beta}\,\phi_{0}}\right)-\left(\frac{\phi_{0}}{M_{P}}\right)^{2}\ln\left(\frac{\phi_{0}-\phi}{M_{P}}\right)^{2}\Bigg{]}\Bigg{|}^{\phi_{\star}}_{\phi_{\text{end}}}.$
(10)

In this work, we consider $50\lesssim N_{\star}\lesssim 65$ Liddle:2003as ;
Drees:2022aea . The power spectrum of the curvature perturbation
$\mathcal{P}_{\zeta}$, the spectral index $n_{s}$, its running $\alpha$,444The
running, defined as $\alpha\equiv\frac{dn_{s}}{d\ln k}$, describes the scale
dependence of the spectral index. and the tensor-to-scalar ratio $r$ are
respectively given by Drees:2021wgd ; Bernal:2021qrl ; Drees:2022aea

$\displaystyle\mathcal{P}_{\zeta}$
$\displaystyle=\frac{1}{24\pi^{2}\,\epsilon_{V}}\,\frac{V}{M_{P}^{4}}\,,$ (11)
$\displaystyle n_{s}$ $\displaystyle=1-6\epsilon_{V}+2\eta_{V}\,,$ (12)
$\displaystyle\alpha$
$\displaystyle=16\,\epsilon_{V}\,\eta_{V}-24\,\epsilon_{V}^{2}-2\,\xi_{V}^{2}\,,$
(13) $\displaystyle r$ $\displaystyle=16\,\epsilon_{V},$ (14)

in the SR formalism. The observables are constrained by the Planck 2018
measurements including baryonic acoustic oscillations at the pivot scale
$k_{\star}=0.05$ Mpc-1 to Planck:2018vyg

$\mathcal{P}_{\zeta}=(2.1\pm 0.1)\times 10^{-9}\,,\quad n_{s}=0.9659\pm
0.0040\,,\quad\alpha=-0.0041\pm 0.0067\,.$ (15)

The most recent constraint on $r$ is obtained from BICEP/Keck 2018 together
with Planck data BICEP:2021xfz :

$r<0.035,\qquad\text{at 95\%~{}C.L.}$ (16)

Note that the model contains four parameters that play a role in inflation,
namely $d$, $\beta$, $\phi_{\star}$, and $\phi_{0}$. The value of
$\phi_{\star}$ is fixed by setting the number of $e$-folds $N_{\star}$, while
the three remaining degrees of freedom must satisfy the constraints in Eqs.
(15) and (16). The parameter space allowed for $\phi_{0}$ is found to be
Drees:2021wgd ; Drees:2022aea

$3\times 10^{-5}~{}M_{P}\lesssim\phi_{0}\lesssim 21.5~{}M_{P}\,.$ (17)

For smaller $\phi_{0}$, it is required that $\phi_{\star}$ be close to
$\phi_{0}$ to fit the CMB data. In this case, even a small loop correction to
the potential from the inflaton couplings could affect the flatness required
by CMB observations Drees:2021wgd . On the other hand, for very large values
of $\phi_{0}$, $\phi_{\star}\ll\phi_{0}$, the first term in Eq. (8) dominates.
This implies that when $\phi_{0}$ is larger than a critical value, the
potential becomes near quadratic, $V(\phi)\simeq
2\,d\,\phi_{0}^{2}\,\phi^{2}$. The latter is in tension with the upper bound
on $r$ from BICEP/Keck 2018 data BICEP:2021xfz , which allows us to derive an
upper bound on $\phi_{0}$.

Having presented the allowed range of $\phi_{0}$ in Eq. (17), we can now
discuss the parameter space of the inflaton mass $m_{\phi}$, which is the
second derivative of the potential in Eq. (8),

$m_{\phi}=2\sqrt{d}\,\phi_{0}\simeq 5.14\times
10^{-8}\left(\frac{\phi_{0}^{2}}{M_{P}}\right),$ (18)

where the last equality was obtained by using $d\simeq 6.61\times
10^{-16}~{}(\phi_{0}/M_{P})^{2}$.555The parameter $d$ controls the overall
amplitude of the power spectrum $\mathcal{P}_{\zeta}$ and was obtained in Ref.
Drees:2021wgd  by fixing $\mathcal{P}_{\zeta}$ to the central value in Eq.
(15). This analytical expression remains a good approximation for
$\phi_{0}\lesssim 5\,M_{P}$. For $\phi_{0}\gtrsim 5\,M_{P}$, the range of $d$
is found to be $10^{-13}\lesssim d\lesssim 10^{-14}$ Drees:2022aea . From the
constraints on $d$ and $\phi_{0}$ presented here and in Eq. (17),
respectively, we find that the inflaton mass can acquire the following range
of possible values

$10^{2}~{}\text{GeV}\lesssim m_{\phi}\lesssim 10^{13}~{}\text{GeV}\,.$ (19)

Future CMB experiments could further constrain the upper limit on the inflaton
mass. The prediction for the running $\alpha$ of the spectral index is
negative and lies in the range Drees:2021wgd ; Drees:2022aea

$\alpha\in[-\mathcal{O}(10^{-2}),\,-\mathcal{O}(10^{-3})]\,.$ (20)

Figure 1: Left panel: Predictions in the $[n_{s},\,r]$ plane for four
benchmark values of $\phi_{0}$ and fixed $\beta$. Right panel: Impact of
$\beta$ on inflationary predictions for fixed $\phi_{0}$. The small (large)
red dot corresponds to $N_{\star}=50$ ($N_{\star}=65$). The blue shaded area
is allowed by the combination of BICEP/Keck and Planck data BICEP:2021xfz ,
while the gray line corresponds to a projected sensitivity of CMB-S4
Abazajian:2019eic .

In the left panel of Fig. 1, we show predictions in the $[n_{s},\,r]$ plane
for four benchmark scenarios for $50\lesssim N_{\star}\lesssim 65$ indicated
by a small and large red dot, respectively: $\phi_{0}=21.5~{}M_{P}$ and
$\beta=1\times 10^{-3}$ (solid red), $\phi_{0}=11~{}M_{P}$ and $\beta=19\times
10^{-3}$ (red dashed), $\phi_{0}=8~{}M_{P}$ and $\beta=6\times 10^{-3}$ (dash-
dotted red) and $\phi_{0}=5~{}M_{P}$ with $\beta=10^{-3}$ (dotted red). Larger
$\phi_{0}$ corresponds to a higher inflationary scale with higher predictions
for $r$. For fixed $\phi_{0}$, a larger $\beta$ implies a steeper potential,
giving rise to a larger prediction for $r$, as shown in the right panel of
Fig. 1. The normalization $d$ of the potential has no impact on these results.
The blue shaded area indicates the parameter space allowed by the combination
of BICEP/Keck 2018 and Planck data BICEP:2021xfz . We note that future CMB
experiments, such as CMB-S4 Abazajian:2019eic  could offer more stringent
constraints on $\phi_{0}$, as shown in gray. As demonstrated above, polynomial
inflation, depending on the parameters $\phi_{0}$, $d$, $\beta$, and
$\phi_{\star}$, can well accommodate recent CMB data. In the large-field
scenario, the potential is dominated by the quadratic term when inflation
ends, which is independent of $\beta$. For the focus of this paper, namely the
DM production during reheating, the relevant parameters are the inflaton mass
and the inflationary scale, which both depend on $\phi_{0}$.

###  3.2 (P)reheating

After the end of cosmic inflation, the inflaton oscillates around the minimum
of its potential and transfers energy into a plasma of lighter degrees of
freedom. The inflaton oscillates around $\phi=0$ during reheating, where the
quadratic term in Eq. (8) dominates, leading to
$V(\phi)\simeq\frac{m_{\phi}^{2}}{2}\phi^{2}$. In our setup, the dominant
channel for energy transfer is the decay of inflatons into pairs of SM Higgs
bosons with mass $m_{H}$. The corresponding decay rate reads

$\Gamma_{\phi}=\frac{1}{8\pi}\,\frac{\mu^{2}}{m_{\phi}}\left[1-\left(\frac{2\,m_{H}}{m_{\phi}}\right)^{2}\right]^{1/2}.$
(21)

The evolution of the inflaton energy density, $\rho_{\phi}$, and the radiation
energy density, $\rho_{R}$, are governed by the following set of Boltzmann
equations

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

where the Hubble expansion rate $\mathcal{H}$ is determined from the Friedmann
equation

$\mathcal{H}^{2}=\frac{\rho_{R}+\rho_{\phi}}{3M_{P}^{2}}\,.$ (24)

An approximate analytical solution for the SM energy density during reheating
is

$\rho_{\phi}(a)\simeq\rho_{\phi}(a_{\text{end}})\left(\frac{a_{\text{end}}}{a}\right)^{3},$
(25)

which is valid as long as $\mathcal{H}\gg\Gamma_{\phi}$. Here,
$a_{\text{end}}$ corresponds to the scale factor at the end of inflation and
the beginning of reheating, $\rho_{R}(a_{\text{end}})=0$ and
$\rho_{\phi}(a_{\text{end}})=3\,M_{P}^{2}\,\mathcal{H}^{2}(a_{\text{end}})$.
The inflationary Hubble scale $\mathcal{H}(a_{\text{end}})$ is given by the
Friedmann equation $\mathcal{H}^{2}(a_{\text{end}})\simeq
V(a_{\text{end}})/(3M_{P}^{2})$. Using Eq. (25), the solution for $\rho_{R}$
during reheating takes the form

$\rho_{R}(a)=\frac{6}{5}M_{P}^{2}\,\Gamma_{\phi}\,\mathcal{H}(a_{\text{end}})\left(\frac{a_{\text{end}}}{a}\right)^{\frac{3}{2}}\left[1-\left(\frac{a_{\text{end}}}{a}\right)^{\frac{5}{2}}\right].$
(26)

While we have verified the validity of the above analytical approximations, we
obtain our numerical results by solving the full set of Boltzmann equations
without relying on any approximations.

The SM radiation energy density is defined as a function of the temperature
$T$ of the SM bath as666As the produced particles interact with the plasma,
they eventually thermalize Allahverdi:2010xz ; Amin:2014eta . In this work, we
assume that the thermalization process occurs instantaneously, which is a good
approximation if the reheating process is not dominated by Planck-suppressed
operators Harigaya:2013vwa . If the instantaneous-thermalization assumption is
violated, then the decay products could follow initial distributions with
lower occupation numbers and harder momenta Harigaya:2013vwa ; Ellis:2015jpg ;
Garcia:2018wtq .

$\rho_{R}(T)=\frac{\pi^{2}}{30}\,g_{\star}(T)\,T^{4},$ (27)

where $g_{\star}(T)$ corresponds to the number of relativistic degrees of
freedom contributing to the SM energy density. The end of reheating is
determined by the onset of SM radiation domination, at a temperature
$T=T_{\text{rh}}$. We define the reheating temperature by the equality
$\mathcal{H}(T_{\text{rh}})\equiv\frac{2}{3}\Gamma_{\phi}$, assuming that the
inflaton decays away at a time scale $t=1/\Gamma_{\phi}$,

$T_{\text{rh}}^{2}=\frac{2}{\pi}\left(\frac{10}{g_{\star}(T_{\text{rh}})}\right)^{\frac{1}{2}}M_{P}\,\Gamma_{\phi}\,.$
(28)

We require $T_{\text{rh}}>T_{\text{BBN}}\simeq 4$ MeV to leave the standard
BBN predictions unaffected by our setup Sarkar:1995dd ; Kawasaki:2000en ;
Hannestad:2004px ; DeBernardis:2008zz ; deSalas:2015glj . Furthermore, the SM
energy density in Eq. (26) has a maximum at
$a=a_{\text{max}}\equiv(8/3)^{2/5}\,a_{\text{end}}$, and the corresponding
maximum temperature $T_{\text{max}}$ Giudice:2000ex  is given by

$T_{\text{max}}^{4}=\frac{60}{\pi^{2}\,g_{\star}(T_{\text{max}})}\left(\frac{3}{8}\right)^{\frac{8}{5}}M_{P}^{2}\,\Gamma_{\phi}\,\mathcal{H}(a_{\text{end}})\,,$
(29)

which is an important parameter in the UV freeze-in process to be discussed in
the next section. Finally, considering that during reheating the temperature
scales as $T(a)\propto a^{-3/8}$ (cf. Eq. (26)), it follows that the Hubble
expansion rate can be expressed as

$\mathcal{H}(T)=\frac{\pi}{3}\,\frac{T^{2}}{M_{P}}\times\begin{dcases}\sqrt{\frac{g_{\star}(T_{\text{rh}})}{10}}\left(\frac{T}{T_{\text{rh}}}\right)^{2}&\text{
for }T_{\text{max}}>T\geq T_{\text{rh}}\,,\\\
\sqrt{\frac{g_{\star}(T)}{10}}&\text{ for }T_{\text{rh}}\geq T\,.\end{dcases}$
(30)

This result will later be used to analytically solve the Boltzmann equation
governing DM production.

Figure 2: Evolution of energy densities (left) and SM temperature (right) for
$T_{\text{rh}}=10^{14}$ GeV (top) and $T_{\text{rh}}=10^{10}$ GeV (bottom)
with for $\phi_{0}=21.5~{}M_{P}$ and $m_{\phi}\simeq 1.5\times 10^{13}$ GeV.

To illustrate the preceding discussion, in Fig. 2 we depict the evolution of
the radiation and inflaton energy densities (left) and the SM temperature
(right), for the largest possible value $\phi_{0}=21.5~{}M_{P}$ and
$m_{\phi}\simeq 1.5\times 10^{13}$ GeV (the latter is obtained from Eq. (18)
using $d\simeq 2\times 10^{-14}$). We have numerically solved the coupled
Boltzmann equations (22) and (23) for two scenarios $T_{\text{rh}}=10^{14}$
GeV (upper panel) and $T_{\text{rh}}=10^{10}$ GeV (lower panel). The former
value is the largest possible reheating temperature allowed in the large field
inflation scenario as we will see in the following, and the latter illustrates
a scenario featuring a reheating temperature lower than the inflaton mass. In
a bosonic reheating process, it is possible to reach temperatures higher than
the inflaton mass due to soft bosons rescattering and forming harder degrees
of freedom Drees:2021wgd . For fixed reheating temperature, the inflaton decay
rate appearing in the Boltzmann equation is obtained via
$\Gamma_{\phi}=3/2\mathcal{H}(T_{\text{rh}})$ with
$\mathcal{H}(T_{\text{rh}})=\frac{\pi}{3}\sqrt{\frac{g_{\star}(T_{\text{rh}})}{10}}\,\frac{T_{\text{rh}}^{2}}{M_{P}}$.
The inflaton does not thermalize during reheating even though the temperature
can be larger than the inflaton mass, since the temperature-independent
inflaton interaction rate $\Gamma_{\phi}$ is always smaller than the Hubble
rate during reheating. Given its relevance to the discussion of the DM
production, we emphasize that during reheating, the maximum temperature
$T_{\text{max}}$ can be much higher than $T_{\text{rh}}$ Giudice:2000ex , as
shown in the lower right panel of Fig. 2.

In general, inflaton couplings to other particles could spoil the flatness of
the potential during inflation. In the polynomial inflationary setup, the
radiative stability of the inflaton potential during inflation requires a
suppressed inflaton-Higgs coupling, which in turn requires Drees:2021wgd

$\left|\left(\frac{\mu}{\phi_{0}}\right)^{2}\ln\left(\frac{\mu}{\phi_{0}}\right)-\left(\frac{\mu}{\phi_{0}}\right)^{2}\right|<64\pi^{2}d\,\beta\,,$
(31)

and translates into an upper bound on the coupling $\mu$. The constraint on
$\mu$ shown in Eq. (31) together with the upper bound on $\phi_{0}=21.5M_{P}$
in Eq. (17) gives rise to an upper bound on the maximal reheating temperature
Drees:2022aea

$T_{\text{rh}}\lesssim 10^{14}~{}\text{GeV}\,.$ (32)

We emphasize that the upper bound on $T_{\text{rh}}$ depends on $\phi_{0}$ and
the couplings involved for the reheating. For smaller $\phi_{0}$, e.g.
$\phi_{0}=M_{P}$, it has been shown that $T_{\text{rh}}\lesssim 10^{11}$ GeV
Drees:2021wgd . Similarly, DM would also contribute (negatively) to the
inflaton potential at loop level, which could potentially spoil inflationary
predictions. To keep radiative corrections under control, $y_{\chi}$ must
satisfy Drees:2021wgd

$\left|y_{\chi}^{4}-3\,y_{\chi}^{4}\,\ln\left(y_{\chi}^{2}\right)\right|<64\pi^{2}\,d\,\beta\,.$
(33)

Assuming the upper bound on $\phi_{0}=21.5M_{P}$ and taking into account the
allowed range of $d$ and $\beta$ within the $2\sigma$ range of the BICEP 2018
data, we obtain

$y_{\chi}\lesssim 10^{-4}\,.$ (34)

We note that within our framework, non-perturbative preheating and pure
gravitational reheating are subdominant to perturbative decay. The trilinear
coupling, $\mu\,\phi\,|H|^{2}$, which induces a tachyonic squared-mass
$m^{2}_{H}\sim\mu\,\phi$ for the field $H$ once $\phi$ becomes negative after
crossing zero, tends to make preheating efficient. However, the quartic self-
interaction term, $\lambda_{H}\,H^{4}$, introduces a positive squared mass
term $m^{2}_{H}\sim\lambda_{H}\left\langle H^{2}\right\rangle$, where
$\left\langle H^{2}\right\rangle$ denotes the variance of $H$. This
backreaction effectively counteracts the tachyonic effects, quickly
terminating preheating and rendering it inefficient Dufaux:2006ee . Finally,
we comment on the possibility of gravitational reheating. It has been
established that rather steep inflaton potentials $V\sim\phi^{p}$ with $p>9$
are required at the reheating stage for successful gravitational reheating
Clery:2021bwz ; Haque:2022kez ; Barman:2022qgt ; Haque:2023zhb .
Interestingly, this bound can be relaxed to $p>4$ by introducing a non-minimal
coupling between gravity and two inflatons Clery:2022wib , or even to $p\geq
2$ by coupling gravity non-minimally to a single inflaton Barman:2023opy . In
our case, the inflaton oscillates around a quadratic potential
$V(\phi)\propto\frac{1}{2}\,m_{\phi}^{2}\,\phi^{2}$ during reheating and is
minimally coupled to gravity, making gravitational reheating inefficient. It
is therefore ignored in the following.

##  4 Dark Matter Production

In the following, we aim to quantify DM production within the large-field
inflationary model during and after reheating. The evolution of the DM number
density $n$ can be tracked using the Boltzmann equation

$\frac{dn}{dt}+3\,\mathcal{H}\,n=\gamma\,,$ (35)

where $\gamma$ denotes the DM production rate density. Taking into account
that during reheating the SM entropy is not conserved, it is convenient to
introduce a comoving number density $N\equiv n\,a^{3}$, with which the
Boltzmann Eq. (35) can be rewritten as

$\frac{dN}{da}=\frac{a^{2}\,\gamma}{\mathcal{H}}\,,$ (36)

or equivalently as

$\frac{dN}{dT}=\begin{dcases}-\frac{8}{\pi}\sqrt{\frac{10}{g_{\star}(T_{\text{rh}})}}\,\frac{M_{P}\,T_{\text{rh}}^{10}}{T^{13}}a^{3}(T_{\text{rh}})\,\gamma&\text{
for }T_{\text{max}}>T\geq T_{\text{rh}}\,,\\\
-\frac{3}{\pi}\sqrt{\frac{10}{g_{\star}(T)}}\,\frac{g_{\star
s}(T_{\text{rh}})}{g_{\star
s}(T)}\,\frac{M_{P}\,T_{\text{rh}}^{3}}{T^{6}}\,a^{3}(T_{\text{rh}})\,\gamma&\text{
for }T_{\text{rh}}\geq T\,,\end{dcases}$ (37)

using Eq. (30). We solve Eq. (37) to obtain the number density. To facilitate
a comparison with observational data, we introduce the DM yield $Y(T)\equiv
n(T)/s(T)$ as a function of the SM entropy density
$s(T)\equiv\frac{2\pi^{2}}{45}\,g_{\star s}\,T^{3}$, where $g_{\star s}(T)$ is
the number of relativistic degrees of freedom that contribute to the SM
entropy. To fit the observed relic abundance, the DM yield has been fixed such
that

$m_{\chi}\,Y_{0}=\Omega_{\chi}h^{2}\,\frac{1}{s_{0}}\,\frac{\rho_{c}}{h^{2}}\simeq
4.3\times 10^{-10}~{}\text{GeV}\,,$ (38)

with $\rho_{c}\simeq 1.05\times 10^{-5}\,h^{2}$ GeV/cm3 the critical energy
density of the Universe, $s_{0}\simeq 2.9\times 10^{3}$ cm-3 the present
entropy density ParticleDataGroup:2020ssz , and $\Omega_{\chi}h^{2}\simeq
0.12$ Planck:2018vyg .

Figure 3: Feynman diagram for fermionic DM production.

The relevant tree-level processes that contribute to DM production are
depicted in Fig. 3, and include

  * •

Inflaton decay (first row),

  * •

Inflaton scattering (second row),

  * •

SM particle scattering (third row),

which we will investigate in the following three subsections, respectively.

###  4.1 Inflaton Decay

The interaction rate density associated with the production of pairs of DM
particles from inflaton decays (through the trilinear interaction in Eq. (5))
is given by

$\gamma=2\,n_{\phi}\,\Gamma_{\phi}\,\text{Br}\,,$ (39)

where $n_{\phi}=\rho_{\phi}/m_{\phi}$ denotes the inflaton number density, and
$\Gamma_{\phi}$ is given in Eq. (21). Note that two DM particles are produced
in each inflaton decay, which is accounted for by the factor 2 in Eq. (39).
Additionally, the branching ratio, Br, for the decay of the inflaton into DM
is given by

$\text{Br}\equiv\frac{\Gamma_{\phi\to\bar{\chi}\chi}}{\Gamma_{\phi\to
H^{\dagger}H}+\Gamma_{\phi\to\bar{\chi}\chi}}\simeq
y_{\chi}^{2}\left(\frac{m_{\phi}}{\mu}\right)^{2}\left(1-\frac{4\,m_{\chi}^{2}}{m_{\phi}^{2}}\right)^{3/2},$
(40)

where we have neglected the Higgs mass. With the interaction rate in Eq. (39),
one can solve Eq. (37) during reheating, i.e. for $T_{\text{max}}\geq T\geq
T_{\text{rh}}$, and obtain

$N(a_{\text{rh}})\equiv N(T_{\text{rh}})\
\simeq\frac{2\pi\,g_{\star}(T_{\text{rh}})}{15}\sqrt{\frac{10}{g_{\star}(T_{\text{rh}})}}\,\text{Br}\,\frac{M_{P}\,T_{\text{rh}}^{2}\,\Gamma_{\phi}}{m_{\phi}}\,a_{\text{rh}}^{3}\,,$
(41)

from which we arrive at the DM yield at present

$Y_{0}=Y_{\text{rh}}=\frac{N(a_{\text{rh}})}{s(T_{\text{rh}})\,a_{\text{rh}}^{3}}\simeq\frac{3}{2}\,\frac{g_{\star}(T_{\text{rh}})}{g_{\star
s}(T_{\text{rh}})}\,\frac{T_{\text{rh}}}{m_{\phi}}\,\text{Br}\simeq\frac{3\,y_{\chi}^{2}}{8\pi^{2}}\,\frac{\sqrt{10\,g_{\star}(T_{\text{rh}})}}{g_{\star
s}(T_{\text{rh}})}\,\frac{M_{P}}{T_{\text{rh}}}\,.$ (42)

Note that DM production from inflaton decay is exponentially suppressed after
reheating and can be neglected.

Figure 4: Parameter space that reproduces the observed DM abundance for DM
produced through 2-body decays of the inflaton for different values of the
Yukawa coupling $y_{\chi}$ with $m_{\phi}=1.5\times 10^{13}~{}\text{GeV}$ and
$\phi_{0}=21.5~{}M_{P}$ (thick black lines). The colored bands are in tension
with BBN (orange), Lyman-$\alpha$ (red) or generate radiative corrections from
Higgs (blue) or DM loops (green) to the potential that could spoil the
flatness required for successful inflation. The vertical red dotted line
corresponds to $m_{\chi}=m_{\phi}/2$.

In Fig. 4, we show with thick black lines the parameter space
$[m_{\chi},\,T_{\text{rh}}]$ required to reproduce the observed DM relic
density through inflaton decays, for several benchmark values for the Yukawa
coupling $y_{\chi}$. As an example, we have again fixed
$\phi_{0}=21.5~{}M_{P}$, which corresponds to an inflaton mass $m_{\phi}\simeq
1.5\times 10^{13}$ GeV for $d=2\times 10^{-14}$. Energy conservation dictates
that only DM with masses up to $m_{\phi}/2$ can be produced from inflaton
decays, as illustrated by the vertical red line in Fig. 4. Several constraints
to the model parameters are overlaid in color. Low reheating temperatures with
$T_{\text{rh}}\lesssim 4~{}\text{MeV}$ Sarkar:1995dd ; Kawasaki:2000en ;
Hannestad:2004px ; DeBernardis:2008zz ; deSalas:2015glj  are in tension with
BBN, whereas high reheating temperatures and large Yukawa couplings are
disfavored due to the destabilization of the inflaton potential when including
loop corrections, cf. Eqs. (32) and (34). Finally, small masses are not
allowed due to the Lyman-$\alpha$ constraints: The DM particles produced from
inflaton decays had an initial energy $E_{\chi}=m_{\phi}/2$ at the time of
production. Due to redshifting, the present momentum of the DM particles is
$p_{\chi}(a_{0})=E_{\chi}\times a_{\text{rh}}/a_{0}$ with $a_{0}$ and
$a_{\text{rh}}$ being the scale factors at present and at the end of
reheating, respectively. As a cold DM candidate, it is required that the DM
velocity at present is smaller than that of warm DM constrained by
Lyman-$\alpha$, namely $p_{\chi}(a_{0})/m_{\chi}\lesssim 1.8\times 10^{-8}$
Masina:2020xhk , which leads to Bernal:2021qrl

$\frac{m_{\chi}}{\text{keV}}\gtrsim\frac{m_{\phi}}{T_{\text{rh}}}\,.$ (43)

The Lyman-$\alpha$ constraint and the condition $y_{\chi}<10^{-4}$ can be
combined to give a lower bound on the DM mass

$m_{\chi}\gtrsim 4.8~{}\text{keV}.$ (44)

We note that compared to the previous analysis in a small-field model
Bernal:2021qrl , where $m_{\chi}\lesssim 10^{11}$ GeV, DM can be as large as
$m_{\chi}\simeq m_{\phi}/2\sim 10^{13}$ GeV in the large-field scenario.
However, other production channels, e.g. gravitational production, can also be
efficient in the large-field scenario and will modify the results presented in
Fig. 4, as we demonstrate below.

###  4.2 Inflaton Scattering

In addition to the production from inflaton decays, DM particles are also
generated from inflaton scatterings, as shown in the second row of Fig. 3. The
relevant processes are $s$-channel annihilations mediated by a graviton or an
inflaton (due to the $\phi^{3}$ term in the potential) and $t$-channel
annihilations mediated by DM. Compared to decays, the inflaton-mediated
$s$-channel and DM-mediated $t$-channel annihilations are always suppressed
because of extra couplings appearing in the cross section. We thus only focus
on gravitational annihilation, which dominates the production if the direct
coupling between the inflaton and the DM is sufficiently suppressed.

For the DM production from gravitational annihilation of inflatons, the
interaction-rate density is given by Mambrini:2021zpp ; Bernal:2021kaj ;
Barman:2021ugy ; Barman:2022tzk

$\gamma=\frac{\rho_{\phi}^{2}}{m_{\phi}^{2}}\,\frac{m_{\chi}^{2}}{64\pi
M_{P}^{4}}\left(1-\frac{m_{\chi}^{2}}{m_{\phi}^{2}}\right)^{3/2},$ (45)

where $m_{\chi}$ in the second term arises from spin summation in the final
states. A detailed derivation of Eq. (45) can be found in appendix A. The
corresponding DM yield at present is

$Y_{0}\simeq\frac{g_{\star}}{128\,g_{\star
s}}\sqrt{\frac{g_{\star}}{10}}\left(\frac{m_{\chi}}{m_{\phi}}\right)^{2}\left(\frac{T_{\text{rh}}}{M_{P}}\right)^{3}\left[\left(\frac{T_{\text{max}}}{T_{\text{rh}}}\right)^{4}-1\right]\left(1-\frac{m_{\chi}^{2}}{m_{\phi}^{2}}\right)^{3/2}.$
(46)

Note that the enhancement factor
$\left[(T_{\text{max}}/T_{\text{rh}})^{4}-1\right]$ arises from the
noninstantaneous nature of reheating. In large field polynomial inflation, a
large ratio $T_{\text{max}}/T_{\text{rh}}$ naturally occurs (cf. Fig. 2), and
the corresponding large enhancement of gravitational DM production from
inflaton annihilation can successfully account for the observed relic
abundance, as depicted by the solid black line in Fig. 5. For inflaton masses
of the order $m_{\phi}\sim 10^{13}$ GeV, gravitational production is efficient
even when the reheating temperature is as low as $T_{\text{rh}}\sim 10^{5}$
GeV in the regime where the mass of the DM reaches its kinematic threshold.
Note that when $m_{\chi}=m_{\phi}$, the production is closed from the phase
space. By comparing Eq. (46) and Eq. (42), we find that for a given reheating
temperature $T_{\text{rh}}$, inflaton mass $m_{\phi}$, and DM mass
$m_{\chi}(\lesssim m_{\phi}/2)$, the gravitational annihilation dominates over
the decay if

$y_{\chi}\lesssim\frac{\pi}{4}\,\sqrt{\frac{g_{\star}}{30}}\,\frac{m_{\chi}}{m_{\phi}}\left(\frac{T_{\text{max}}}{M_{P}}\right)^{2}.$
(47)

Figure 5: Parameter space that reproduces the observed DM abundance for DM
produced through gravitational scatterings of inflatons or SM particles with
$m_{\phi}=1.5\times 10^{13}~{}\text{GeV}$ and $\phi_{0}=21.5M_{P}$ (thick
black lines). The colored bands are in tension with BBN (orange),
Lyman-$\alpha$ (red) or generate a radiative-unstable potential (blue). The
vertical red dotted lines correspond to $m_{\chi}=m_{\phi}/2$,
$m_{\chi}=m_{\phi}$, $m_{\chi}=T_{\text{rh}}$ and $m_{\chi}=T_{\text{max}}$.

So far we have focused on production channels for DM with mass
$m_{\chi}\lesssim m_{\phi}$. As mentioned previously, the temperature during
reheating in the large-field scenario can be larger than the inflaton mass.
This means that it is possible to source DM with a mass as large as
$T_{\text{max}}$ through the scattering of SM particles from the plasma, as we
demonstrate in the next section.

###  4.3 SM Particles Scattering

In the following, we explore the DM production from scatterings of SM
particles (last row of Fig. 3), mediated by either the graviton or inflaton.
These processes allow the production of DM with $m_{\chi}\sim T_{\text{max}}$
and feature a UV freeze-in behavior due to the suppression of a heavy scale.

####  4.3.1 Graviton Mediation

For gravitational DM production from the annihilation of SM particles in the
thermal plasma, the interaction rate is given by Garny:2015sjg ; Tang:2017hvq
; Garny:2017kha ; Bernal:2018qlk

$\gamma\simeq\alpha_{f}\frac{T^{8}}{M_{P}^{4}}\,,$ (48)

where $\alpha_{f}\simeq 1.1\times 10^{-3}$ is obtained for fermionic DM by
summing over the contributions from all SM degrees of freedom in the thermal
plasma Garny:2017kha . The corresponding DM yield at present is given by

$Y_{0}\simeq\begin{dcases}\frac{45\,\alpha_{f}}{2\pi^{3}g_{\star
s}}\sqrt{\frac{10}{g_{\star}}}\frac{T_{\text{rh}}^{3}}{M_{P}^{3}}\left\\{\left[1-\left(\frac{m_{\chi}}{T_{\text{rh}}}\right)^{3}\right]+2\left[1-\left(\frac{T_{\text{rh}}}{T_{\text{max}}}\right)^{4}\right]\right\\}\simeq\frac{135\,\alpha_{f}}{2\pi^{3}g_{\star
s}}\sqrt{\frac{10}{g_{\star}}}\frac{T_{\text{rh}}^{3}}{M_{P}^{3}}\\\ \hskip
239.00298pt\text{for }m_{\chi}\lesssim T_{\text{rh}}\,,\\\
\frac{45\,\alpha_{f}}{\pi^{3}g_{\star
s}}\sqrt{\frac{10}{g_{\star}}}\,\frac{T_{\text{rh}}^{7}}{M_{P}^{3}\,m_{\chi}^{4}}\left[1-\left(\frac{m_{\chi}}{T_{\text{max}}}\right)^{4}\right]\qquad\qquad\qquad\text{for
}T_{\text{rh}}\lesssim m_{\chi}\lesssim T_{\text{max}}\,.\end{dcases}$ (49)

DM lighter than $T_{\text{rh}}$ is produced both during and after reheating,
while DM heavier than $T_{\text{rh}}$ (and lighter than $T_{\text{max}}$) is
produced only during reheating. For $m_{\chi}\lesssim T_{\text{rh}}$, the
gravitational scattering of SM particles is subdominant compared to the
channels studied in the two preceding subsections. However, in contrast to the
channels considered in Sections 4.1 and 4.2, gravitational scattering of SM
particles can account for the observed relic abundance even if
$m_{\chi}\gtrsim m_{\phi}$, as shown by the solid black line in the upper
right corner of Fig. 5. Consequently, DM as heavy as $m_{\chi}\simeq
T_{\text{max}}\simeq 4\times 10^{14}$ GeV can be produced in the large-field
polynomial inflation scenario, which is three orders of magnitude larger than
the upper limit in the small-field scenario.

####  4.3.2 Inflaton Mediation

Finally, we focus on the inflaton-mediated UV freeze-in production of DM from
scatterings of SM particles, as depicted in the second diagram of the last row
of Fig. 3. The annihilation cross section $\sigma$ for this process is given
by

$\sigma(\mathfrak{s})\simeq\frac{1}{8\pi}\frac{y^{2}_{\chi}\,\mu^{2}}{(\mathfrak{s}-m_{\phi}^{2})^{2}+m_{\phi}^{2}\,\Gamma_{\phi}^{2}}\left(1-\frac{4m_{\chi}^{2}}{\mathfrak{s}}\right)^{3/2},$
(50)

where $\sqrt{\mathfrak{s}}$ corresponds to the center-of-mass energy.777We
take the Higgs to be massless before the onset of the electroweak phase
transition. The thermal mass of the Higgs is dominated by
$\sim\frac{1}{4}\,y_{t}\,T$ with $y_{t}$ being the quark-top Yukawa coupling
Quiros:1999jp , which is negligible compared to the momentum of the thermal
Higgs. Since the SM temperature can be higher than the inflaton mass, the DM
can be produced resonantly during reheating Barman:2024mqo . The corresponding
interaction rate density is

$\gamma(T)\simeq\frac{T}{8\pi^{4}}\int_{4\,m_{\chi}^{2}}^{\infty}d\mathfrak{s}\,\mathfrak{s}^{3/2}\,\sigma(\mathfrak{s})\,K_{1}\left(\frac{\sqrt{\mathfrak{s}}}{T}\right),$
(51)

and is well approximated by Bernal:2020fvw ; Bernal:2020yqg

$\gamma(T)\simeq\frac{y^{2}_{\chi}\,\mu^{2}}{64\pi^{4}}\times\begin{cases}\frac{32}{\pi}\frac{T^{6}}{m_{\phi}^{4}}&\quad\text{
for }T\ll m_{\phi}\,,\\\
\frac{m_{\phi}^{2}\,T}{\Gamma_{\phi}}\,K_{1}\left(\frac{m_{\phi}}{T}\right)\left[1-\left(\frac{2m_{\chi}}{m_{\phi}}\right)^{2}\right]^{3/2}&\quad\text{
for }T\gtrsim m_{\phi}\,,\end{cases}$ (52)

where $K_{1}$ denotes the modified Bessel function of the first kind. The
analytical approximations shown in Eq. (52) have been validated against the
full numerical result. Finally, as a consistency check, we note that DM does
not thermalize, as $\gamma(T)\ll\mathcal{H}(T)\,n_{\text{eq}}(T)$ with
$n_{\text{eq}}(T)\simeq\frac{4}{\pi^{2}}\,T^{3}$ for $T\gg m_{\chi}$.

The inflaton mass is typically larger than the reheating temperature due to
the radiative stability condition Drees:2021wgd ; Bernal:2021qrl . In the case
where $m_{\phi}\gg T_{\text{rh}}$, which almost always holds for small-field
inflation and sometimes also for large-field inflation, we obtain the
following analytic approximation for the DM yield produced from inflaton-
mediated Higgs decays using the first line of Eq. (52),888Equation (53)
refines the result reported in Ref. Bernal:2021qrl  by a factor $195/135$
accounting for the DM production during reheating, which was not considered
previously.

$Y_{0}\simeq\frac{195\,y_{\chi}^{2}}{4\pi^{8}\,g_{\star
s}}\,\sqrt{\frac{10}{g_{\star}}}\,\frac{\mu^{2}\,T_{\text{rh}}\,M_{P}}{m_{\phi}^{4}}\,\simeq\frac{195\,y_{\chi}^{2}}{g_{\star
s}\,\pi^{6}}\left(\frac{T_{\text{rh}}}{m_{\phi}}\right)^{3}\,.$ (53)

This result is parametrically suppressed by a factor

$\frac{T_{\text{rh}}}{M_{P}}\left(\frac{T_{\text{rh}}}{m_{\phi}}\right)^{3}\ll
1$ (54)

compared to the contribution from inflaton decays (cf. Eq. (42)), and hence
always negligible.

However, in large-field scenarios, it is also possible to have
$T_{\text{rh}}\gtrsim m_{\phi}$. In this case, we find that the ratio of the
DM yield between annihilation and decay is proportional to a factor
$T_{\text{rh}}/(g_{\star}\,m_{\phi})$. This result is obtained by solving Eq.
(37) using the analytical expressions in Eq. (52).999Obtaining a full
analytical result for the yield $Y_{0}$, analogous to Eq. (53), is
considerably more complex for the present case due to the different scales
involved, hence we do not present a corresponding result. Numerically, we find
that in the most extreme case of large-field polynomial inflation consistent
with CMB data, where $T_{\text{rh}}\sim 10^{14}$ GeV, inflaton-mediated Higgs
annihilation provides a percent level correction to the DM relic abundance
produced from inflaton decays.

###  4.4 Combined analysis

In the following, we combine all previously studied processes and investigate
their interplay for different DM masses $m_{\chi}$ and reheating temperatures
$T_{\mathrm{rh}}$. Figure 6 shows, with thick black lines, the parameter space
that reproduces the entire observed DM abundance considering all processes
previously discussed. From left to right, we show the results for different
values of the Yukawa coupling $y_{\chi}$ with $\phi_{0}=21.5~{}M_{P}$ (first
column), $\phi_{0}=5~{}M_{P}$ (second column) and $\phi_{0}=0.8~{}M_{P}$
(third column). The two red dotted vertical lines correspond to
$m_{\chi}=m_{\phi}$ and $m_{\chi}=m_{\phi}/2$.

Figure 6: Parameter space (thick black lines) that reproduce the whole
observed DM abundance taking into account all processes discussed in the text,
for different values of the Yukawa coupling $y_{\chi}$ and
$\phi_{0}=21.5~{}M_{P}$ (first column), $\phi_{0}=5~{}M_{P}$ (second column)
and $\phi_{0}=0.8~{}M_{P}$ (third column). The colored bands are in tension
with BBN (orange), Lyman-$\alpha$ (red) or generate a radiative-unstable
potential (blue). The two red dotted vertical lines correspond to
$m_{\chi}=m_{\phi}$ and $m_{\chi}=m_{\phi}/2$.

We start by considering the largest possible value $\phi_{0}=21.5~{}M_{P}$,
corresponding to $m_{\phi}=1.5\times 10^{13}~{}\text{GeV}$, shown in the four
panels to the left. The upper panel shows the viable parameter space for
$y_{\chi}=10^{-8}$, and the lower panels show how the parameter space changes
as the value of $y_{\chi}$ is lowered. For $y_{\chi}$ larger than a threshold
value, the viable parameter space is continuous in $T_{\text{rh}}$ and
discontinuous in $m_{\chi}$, as seen in the three upper panels. As $y_{\chi}$
is lowered below a threshold value, the parameter space becomes continuous in
$m_{\chi}$ and discontinuous in $T_{\text{rh}}$. We find numerically that the
threshold value is $y_{\chi}\approx 10^{-17}$, which is in agreement with the
analytical threshold value obtained in Eq. (47). For $y_{\chi}\gtrsim
10^{-17}$, the branch to the left corresponds to light DM
($m_{\chi}<m_{\phi}/2$), which is produced via both inflaton decays and
graviton-mediated annihilations of inflatons and SM particles. The branch to
the right corresponds to heavier DM masses ($m_{\phi}/2\lesssim
m_{\chi}<T_{\text{max}}$) and is independent of the value of the Yukawa
coupling $y_{\chi}$ for $T_{\text{rh}}\gtrsim 10^{5}$ GeV, where DM production
is dominated by gravitational processes.101010For $T_{\text{rh}}\lesssim
10^{5}$ GeV, and DM masses very close to the kinematic threshold, DM
production is dominated by inflaton decays. DM is overproduced between the two
branches, while it is underproduced in the two complementary regions of
parameter space. For $y_{\chi}\lesssim 10^{-17}$, there is a continuous range
of viable DM masses that could account for the whole observed relic abundance,
$10^{9}~{}\text{GeV}\lesssim m_{\chi}<T_{\text{max}}$, where both upper and
lower limits come from the intersection of gravitational production from
scattering of SM particles and the constraint from radiative stability. Here,
two branches occur for different values of $T_{\text{rh}}$: If
$T_{\text{rh}}\gtrsim 10^{5}$ GeV, DM is produced by gravitational processes,
whereas if $T_{\text{rh}}\lesssim 10^{5}$ GeV it is generated by the decay of
the inflaton. Note that the black lines correspond to the parameter space
where the correct relic abundance is reached, while the white region within
them corresponds to overproduction. The remaining part leads to
underproduction.

As $\phi_{0}$ decreases, the maximum allowed value for the reheating
temperature $T_{\text{rh}}$ decreases, resulting in suppression of the
gravitational production channels. Consequently, the interplay between
graviton-mediated processes and inflaton decay becomes less prominent and
eventually vanishes for smaller values of $\phi_{0}$, as shown in the second
and third columns of Figure 6. In practice, we find that for $m_{\phi}\lesssim
10^{11}$ GeV or equivalently $\phi_{0}\lesssim 0.8\,M_{P}$, gravitational
channels can no longer account for the DM relic abundance in any part of the
parameter space that is consistent with the radiative bound on
$T_{\text{rh}}$, as seen in the third column of Fig. 6.

##  5 Conclusions

Polynomial inflation is a simple cosmological framework that is consistent
with the cosmic microwave background data and provides testable predictions
for the tensor-to-scalar ratio and the running of the spectral index
Drees:2021wgd ; Drees:2022aea . However, it is desirable that a consistent
cosmological model also explains other observations, including the dark matter
(DM) abundance. Within the polynomial inflation framework, it was shown in
Ref. Bernal:2021qrl  that for small field values, the direct decay of the
inflaton dominates over all other production mechanisms due to the upper bound
on the reheating temperature. However, the maximum reheating temperature can
be up to three orders of magnitude higher in large-field scenarios compared to
the small-field case. Consequently, DM-production channels that are suppressed
in the small-field setup become important in the large-field setup.

We investigated all relevant production channels, including inflaton decays
and inflaton- and graviton-mediated scatterings; see Fig. 3 for the relevant
Feynman diagrams. We find that gravitational inflaton annihilations dominate
over inflaton decays, for inflaton-DM couplings smaller than the threshold
value presented in Eq. (47). When considering all production channels
simultaneously, the parameter space $(m_{\chi},\,T_{\text{rh}})$ where the
whole DM relic abundance is correctly accounted for is generally separated
into two disconnected branches. For $y_{\chi}$ larger than a threshold value,
the two branches are separated by a range of values for $m_{\chi}$ where the
correct relic abundance cannot be accounted for. For $y_{\chi}$ smaller than
the threshold, the two branches are separated by a range of values for
$T_{\text{rh}}$ where the relic abundance cannot be correctly reproduced, see
Fig. 6 for details. Finally, we find that graviton-mediated SM particle
annihilation can account for the whole relic abundance for DM masses as heavy
as the maximum temperature $T_{\text{max}}\sim 10^{14}~{}\text{GeV}$, which is
three orders of magnitude higher than the heaviest DM mass that can be
accommodated in the small field setup.

In summary, the current work offers a comprehensive investigation of DM
production after large-field polynomial inflation. The presented formalism can
be readily extended to other inflationary setups.

###### Acknowledgements.

NB received funding from the Spanish FEDER / MCIU-AEI under the grant
FPA2017-84543-P. MAM acknowledges support from the DFG Collaborative Research
Centre “Neutrinos and Dark Matter in Astro- and Particle Physics” (SFB 1258).
JH, MAM, and YX acknowledge support from the Cluster of Excellence “Precision
Physics, Fundamental Interactions, and Structure of Matter” (PRISMA+ EXC
2118/1) funded by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) within the German Excellence Strategy (Project No. 390831469).

##  Appendix A Gravitational Scattering

In this appendix, we present a detailed derivation of the interaction-rate
density for DM production from graviton-mediated inflaton annihilation,
corresponding to the first diagram in the second row of Fig. 3. We label the
particle momenta as $\phi(p_{1})\,\phi(p_{2})\to
h_{\mu\nu}(q)\to\chi(p_{3})\,\bar{\chi}(p_{4})$. The relevant vertices are
given by Choi:1994ax

$-\frac{i}{M_{P}}\left[p_{1\mu}\,p_{2\nu}+p_{1\nu}\,p_{2\mu}-\eta_{\mu\nu}\left(p_{1}\cdot
p_{2}+m_{\phi}^{2}\right)\right]$ (55)

for $\phi\phi h_{\mu\nu}$, and

$-\frac{i}{4M_{P}}\left[\left(p_{3}-p_{4}\right)_{\mu}\gamma_{\nu}+\left(p_{3}-p_{4}\right)_{\nu}\gamma_{\mu}-2\eta_{\mu\nu}\left(\not{p}_{3}-\not{p}_{4}-2m_{\chi}\right)\right]$
(56)

for $\bar{\chi}\chi h_{\mu\nu}$. The graviton propagator reads

$\Pi^{\mu\nu\rho\sigma}=\frac{1}{2q^{2}}\left(\eta^{\rho\nu}\eta^{\sigma\mu}+\eta^{\rho\mu}\eta^{\sigma\nu}-\eta^{\rho\sigma}\eta^{\mu\nu}\right).$
(57)

From these Feynman rules, we readily obtain the squared matrix element

$\sum_{s_{3},s_{4}}|\mathcal{M}|^{2}=\frac{1}{2\,M_{P}^{4}}m_{\chi}^{2}(m_{\phi}^{2}-m_{\chi}^{2})\,,$
(58)

where $s_{3},s_{4}$ denotes the spins of the final states. With the matrix
element at hand, we arrive at the final result for the interaction-rate
density

$\gamma=\frac{\rho_{\phi}^{2}}{m_{\phi}^{2}}\,\frac{m_{\chi}^{2}}{64\pi
M_{P}^{4}}\left(1-\frac{m_{\chi}^{2}}{m_{\phi}^{2}}\right)^{3/2}.$ (59)

## References

  * (1) A.A. Starobinsky, _A New Type of Isotropic Cosmological Models Without Singularity_ , _Phys. Lett. B_ 91 (1980) 99.
  * (2) A.H. Guth, _The Inflationary Universe: A Possible Solution to the Horizon and Flatness Problems_ , _Phys. Rev. D_ 23 (1981) 347.
  * (3) A.D. Linde, _A New Inflationary Universe Scenario: A Possible Solution of the Horizon, Flatness, Homogeneity, Isotropy and Primordial Monopole Problems_ , _Phys. Lett. B_ 108 (1982) 389.
  * (4) A. Albrecht and P.J. Steinhardt, _Cosmology for Grand Unified Theories with Radiatively Induced Symmetry Breaking_ , _Phys. Rev. Lett._ 48 (1982) 1220.
  * (5) J. Martin, C. Ringeval and V. Vennin, _Encyclopædia Inflationaris_ , _Phys. Dark Univ._ 5-6 (2014) 75 [1303.3787].
  * (6) Planck collaboration, _Planck 2018 results. X. Constraints on inflation_ , _Astron. Astrophys._ 641 (2020) A10 [1807.06211].
  * (7) BICEP2, Keck Array collaboration, _BICEP2 / Keck Array x: Constraints on Primordial Gravitational Waves using Planck, WMAP, and New BICEP2/Keck Observations through the 2015 Season_ , _Phys. Rev. Lett._ 121 (2018) 221301 [1810.05216].
  * (8) H.M. Hodges, G.R. Blumenthal, L.A. Kofman and J.R. Primack, _Nonstandard Primordial Fluctuations From a Polynomial Inflaton Potential_ , _Nucl. Phys. B_ 335 (1990) 197.
  * (9) R. Allahverdi, K. Enqvist, J. García-Bellido and A. Mazumdar, _Gauge invariant MSSM inflaton_ , _Phys. Rev. Lett._ 97 (2006) 191304 [hep-ph/0605035].
  * (10) C. Destri, H.J. de Vega and N.G. Sanchez, _MCMC analysis of WMAP3 and SDSS data points to broken symmetry inflaton potentials and provides a lower bound on the tensor to scalar ratio_ , _Phys. Rev. D_ 77 (2008) 043509 [astro-ph/0703417].
  * (11) K. Nakayama, F. Takahashi and T.T. Yanagida, _Polynomial Chaotic Inflation in the Planck Era_ , _Phys. Lett. B_ 725 (2013) 111 [1303.7315].
  * (12) K. Nakayama, F. Takahashi and T.T. Yanagida, _Polynomial Chaotic Inflation in Supergravity_ , _JCAP_ 08 (2013) 038 [1305.5099].
  * (13) R. Kallosh, A. Linde and A. Westphal, _Chaotic Inflation in Supergravity after Planck and BICEP2_ , _Phys. Rev. D_ 90 (2014) 023534 [1405.0270].
  * (14) T. Li, Z. Sun, C. Tian and L. Wu, _The Renormalizable Three-Term Polynomial Inflation with Large Tensor-to-Scalar Ratio_ , _Eur. Phys. J. C_ 75 (2015) 301 [1407.8063].
  * (15) G. Aslanyan, L.C. Price, J. Adams, T. Bringmann, H.A. Clark, R. Easther et al., _Ultracompact minihalos as probes of inflationary cosmology_ , _Phys. Rev. Lett._ 117 (2016) 141102 [1512.04597].
  * (16) T.-J. Gao and Z.-K. Guo, _Inflection point inflation and dark energy in supergravity_ , _Phys. Rev. D_ 91 (2015) 123502 [1503.05643].
  * (17) N. Musoke and R. Easther, _Expectations for Inflationary Observables: Simple or Natural?_ , _JCAP_ 12 (2017) 032 [1709.01192].
  * (18) M. Drees and Y. Xu, _Small field polynomial inflation: reheating, radiative stability and lower bound_ , _JCAP_ 09 (2021) 012 [2104.03977].
  * (19) M. Drees and Y. Xu, _Large field polynomial inflation: parameter space, predictions and (double) eternal nature_ , _JCAP_ 12 (2022) 005 [2209.07545].
  * (20) Y. Xu, _Polynomial Inflation and Its Aftermath_ , Ph.D. thesis, U. Bonn (main), 2022.
  * (21) M. Drees and Y. Xu, _Parameter space of leptogenesis in polynomial inflation_ , _JCAP_ 04 (2024) 036 [2401.02485].
  * (22) N. Bernal and Y. Xu, _Polynomial inflation and dark matter_ , _Eur. Phys. J. C_ 81 (2021) 877 [2106.03950].
  * (23) F.L. Bezrukov and M. Shaposhnikov, _The Standard Model Higgs boson as the inflaton_ , _Phys. Lett. B_ 659 (2008) 703 [0710.3755].
  * (24) R. Kallosh, A. Linde and D. Roest, _Superconformal Inflationary $\alpha$-Attractors_, _JHEP_ 11 (2013) 198 [1311.0472].
  * (25) S. Izumine and K. Nakayama, _Effects of gravitational particle production on Higgs portal dark matter_ , 2403.05199.
  * (26) 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].
  * (27) D.H. Lyth and A.R. Liddle, _The primordial density perturbation: Cosmology, inflation and the origin of structure_ (2009).
  * (28) A.R. Liddle and S.M. Leach, _How long before the end of inflation were observable perturbations produced?_ , _Phys. Rev. D_ 68 (2003) 103503 [astro-ph/0305263].
  * (29) Planck collaboration, _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6 [1807.06209].
  * (30) BICEP, Keck collaboration, _Improved Constraints on Primordial Gravitational Waves using Planck, WMAP, and BICEP/Keck Observations through the 2018 Observing Season_ , _Phys. Rev. Lett._ 127 (2021) 151301 [2110.00483].
  * (31) K. Abazajian et al., _CMB-S4 Science Case, Reference Design, and Project Plan_ , 1907.04473.
  * (32) R. Allahverdi, R. Brandenberger, F.-Y. Cyr-Racine and A. Mazumdar, _Reheating in Inflationary Cosmology: Theory and Applications_ , _Ann. Rev. Nucl. Part. Sci._ 60 (2010) 27 [1001.2600].
  * (33) M.A. Amin, M.P. Hertzberg, D.I. Kaiser and J. Karouby, _Nonperturbative Dynamics Of Reheating After Inflation: A Review_ , _Int. J. Mod. Phys. D_ 24 (2014) 1530003 [1410.3808].
  * (34) K. Harigaya and K. Mukaida, _Thermalization after/during Reheating_ , _JHEP_ 05 (2014) 006 [1312.3097].
  * (35) J. Ellis, M.A.G. García, D.V. Nanopoulos, K.A. Olive and M. Peloso, _Post-Inflationary Gravitino Production Revisited_ , _JCAP_ 03 (2016) 008 [1512.05701].
  * (36) M.A.G. García and M.A. Amin, _Prethermalization production of dark matter_ , _Phys. Rev. D_ 98 (2018) 103504 [1806.01865].
  * (37) S. Sarkar, _Big bang nucleosynthesis and physics beyond the standard model_ , _Rept. Prog. Phys._ 59 (1996) 1493 [hep-ph/9602260].
  * (38) M. Kawasaki, K. Kohri and N. Sugiyama, _MeV scale reheating temperature and thermalization of neutrino background_ , _Phys. Rev. D_ 62 (2000) 023506 [astro-ph/0002127].
  * (39) S. Hannestad, _What is the lowest possible reheating temperature?_ , _Phys. Rev. D_ 70 (2004) 043506 [astro-ph/0403291].
  * (40) F. De Bernardis, L. Pagano and A. Melchiorri, _New constraints on the reheating temperature of the universe after WMAP-5_ , _Astropart. Phys._ 30 (2008) 192.
  * (41) P.F. de Salas, M. Lattanzi, G. Mangano, G. Miele, S. Pastor and O. Pisanti, _Bounds on very low reheating scenarios after Planck_ , _Phys. Rev. D_ 92 (2015) 123534 [1511.00672].
  * (42) 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].
  * (43) J.F. Dufaux, G.N. Felder, L. Kofman, M. Peloso and D. Podolsky, _Preheating with trilinear interactions: Tachyonic resonance_ , _JCAP_ 07 (2006) 006 [hep-ph/0602144].
  * (44) S. Cléry, Y. Mambrini, K.A. Olive and S. Verner, _Gravitational portals in the early Universe_ , _Phys. Rev. D_ 105 (2022) 075005 [2112.15214].
  * (45) M.R. Haque and D. Maity, _Gravitational reheating_ , _Phys. Rev. D_ 107 (2023) 043531 [2201.02348].
  * (46) B. Barman, S. Cléry, R.T. Co, Y. Mambrini and K.A. Olive, _Gravity as a portal to reheating, leptogenesis and dark matter_ , _JHEP_ 12 (2022) 072 [2210.05716].
  * (47) M.R. Haque, D. Maity and R. Mondal, _Gravitational neutrino reheating_ , _Phys. Rev. D_ 109 (2024) 063543 [2311.07684].
  * (48) S. Cléry, Y. Mambrini, K.A. Olive, A. Shkerin and S. Verner, _Gravitational portals with nonminimal couplings_ , _Phys. Rev. D_ 105 (2022) 095042 [2203.02004].
  * (49) B. Barman, N. Bernal and J. Rubio, _Rescuing gravitational-reheating in chaotic inflation_ , _JCAP_ 05 (2024) 072 [2310.06039].
  * (50) Particle Data Group collaboration, _Review of Particle Physics_ , _PTEP_ 2020 (2020) 083C01.
  * (51) I. Masina, _Dark matter and dark radiation from evaporating primordial black holes_ , _Eur. Phys. J. Plus_ 135 (2020) 552 [2004.04740].
  * (52) Y. Mambrini and K.A. Olive, _Gravitational Production of Dark Matter during Reheating_ , _Phys. Rev. D_ 103 (2021) 115009 [2102.06214].
  * (53) N. Bernal and C.S. Fong, _Dark matter and leptogenesis from gravitational production_ , _JCAP_ 06 (2021) 028 [2103.06896].
  * (54) B. Barman and N. Bernal, _Gravitational SIMPs_ , _JCAP_ 06 (2021) 011 [2104.10699].
  * (55) B. Barman, N. Bernal, Y. Xu and Ó. Zapata, _Ultraviolet freeze-in with a time-dependent inflaton decay_ , _JCAP_ 07 (2022) 019 [2202.12906].
  * (56) M. Garny, M. Sandora and M.S. Sloth, _Planckian Interacting Massive Particles as Dark Matter_ , _Phys. Rev. Lett._ 116 (2016) 101302 [1511.03278].
  * (57) Y. Tang and Y.-L. Wu, _On Thermal Gravitational Contribution to Particle Production and Dark Matter_ , _Phys. Lett. B_ 774 (2017) 676 [1708.05138].
  * (58) 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].
  * (59) 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].
  * (60) M. Quiros, _Finite temperature field theory and phase transitions_ , in _ICTP Summer School in High-Energy Physics and Cosmology_ , pp. 187–259, 1, 1999 [hep-ph/9901312].
  * (61) B. Barman, N. Bernal and Y. Xu, _Resonant Reheating_ , 2404.16090.
  * (62) N. Bernal, A. Donini, M.G. Folgado and N. Rius, _Kaluza-Klein FIMP Dark Matter in Warped Extra-Dimensions_ , _JHEP_ 09 (2020) 142 [2004.14403].
  * (63) N. Bernal, A. Donini, M.G. Folgado and N. Rius, _FIMP Dark Matter in Clockwork/Linear Dilaton Extra-Dimensions_ , _JHEP_ 04 (2021) 061 [2012.10453].