# Isocurvature Constraints on Scalar Dark Matter Production from the Inflaton

Marcos A. G. Garciaa<EMAIL_ADDRESS>Mathias Pierreb
<EMAIL_ADDRESS>Sarunas Vernerc<EMAIL_ADDRESS>a Departamento de
Física Teórica, Instituto de Física, Universidad Nacional Autónoma de México,
Ciudad de México C.P. 04510, Mexico  b Deutsches Elektronen-Synchrotron DESY,
Notkestr. 85, 22607 Hamburg, Germany  c Institute for Fundamental Theory,
Physics Department, University of Florida, Gainesville, FL 32611, USA

(August 27, 2024)

###### Abstract

We investigate the production of a spectator scalar dark matter field that is
directly coupled to the inflaton during inflation and reheating. We consider
two specific inflationary potentials, namely the Starobinsky and T-model of
inflation, which satisfy the constraints on the scalar tilt, $n_{s}$, and
tensor-to-scalar ratio, $r$, measured by the Planck satellite. Excitation of
light scalar dark matter during inflation may result in large isocurvature
perturbations, which can be avoided by inducing a sizable effective dark
matter mass during the inflationary phase. For purely gravitational
production, the Planck isocurvature constraints require the dark matter mass
to be larger than the Hubble scale at horizon exit, with $m_{\chi}\gtrsim
0.5H_{*}$. For small bare dark matter masses $m_{\chi}\ll H_{*}$, these
constraints translate into a lower bound on the dark matter coupling to the
inflaton. We argue that these constraints can applied to a wide class of
single-field slow-roll inflation models. We also derive isocurvature, dark
matter abundance, and Lyman-$\alpha$ constraints on the direct coupling and
bare dark matter mass.

††preprint: DESY-23-027

##  I Introduction

The nature and origin of dark matter (DM) remains one of the biggest puzzles
in fundamental physics. The stringent limits from direct DM detection
searches, such as XENON1T XENON:2018voc , LUX LUX:2016ggv , PandaX
PandaX-4T:2021bab  and LZ LZ:2022ufs , along with the absence of detection
from indirect DM and collider experiments, create tension with the standard
weakly interacting massive particle (WIMP) paradigm without providing clues
about the composition of the invisible fraction of the universe. This tension
provides a strong incentive to explore alternative models of dark matter
Arcadi:2017kky ; Escudero:2016gzx ; Mambrini .

Among the various proposed dark matter production mechanisms, a particularly
minimalistic model assumes that the dark sector of the universe is populated
independently of the visible sector, but via the same mechanism. That is, the
energy density stored in the degrees of freedom that drive primordial
inflation dissipates during the post-inflationary reheating epoch and
populates the dark sector via the freeze-in mechanism Hall:2009bx ;
Bernal:2017kxu . If the dark sector only interacts directly with the inflaton
field, dark matter would never have been in thermal equilibrium with Standard
Model particles, and this could lead to non-perturbative gravitational
particle production Parker:1968mv ; Parker:1969au ; Parker:1971pt ;
Ford:1986sy ; Ema1 ; Ema2 ; Ema3 ; Herranen:2015ima ; Markkanen:2015xuw .
Therefore, its comoving phase space distribution (PSD) would be fully
determined by the end of reheating.

Our previous study Garcia:2022vwm  (see also Chung:1998zb ; Chung:1998ua ;
Kuzmin:1998kk ; Ling:2021zlj ; Kaneta:2022gug ) explored this minimal out-of-
equilibrium scenario for a scalar DM particle. We determined the relic
abundance and structure formation constraints for a wide range of inflaton-DM
coupling values, covering pure gravitational production, weak direct coupling
(perturbative production), and strong direct coupling (non-perturbative
production) regimes. In the weak coupling regime, we found that Lyman-$\alpha$
forest measurements of the matter power spectrum allow DM masses as small as
$\sim 0.3\,{\rm meV}$. However, purely gravitationally produced dark matter
faces a tension with the Cosmic Microwave Background (CMB) bounds on the
amplitude of the isocurvature power spectrum for DM masses smaller than the
instantaneous Hubble scale at the end of inflation Chung:2004nh ; Chung:2011xd
; Chung:2015pga ; Herring:2019hbe ; Padilla:2019fju ; Ling:2021zlj ;
Redi:2022zkt . The current constraints on the isocurature power spectrum from
${\textit{P}lanck}$ are given by $\beta_{\rm
iso}\;\equiv\;\mathcal{P}_{\mathcal{S}}(k_{*})/(\mathcal{P}_{\mathcal{R}}(k_{*})+\mathcal{P}_{\mathcal{S}}(k_{*}))\;<\;0.038$
at the 95% C.L. with the pivot scale $k_{*}=0.05\,{\rm Mpc}^{-1}$
Planck:2018jri , where $\mathcal{P}_{\mathcal{R}}$ and
$\mathcal{P}_{\mathcal{S}}$ denote the curvature and isocurvature power
spectra, respectively.

As argued in Gordon:2000hv ; Herring:2019hbe , isocurvature modes are always
strongly suppressed relative to adiabatic modes at linear order when there is
no dark matter misalignment. A cursory computation shows that the curvature
perturbation $\mathcal{R}$ is conserved in super-horizon scales, consistent
with the absence of late-time isocurvature Renaux-Petel:2014htw . However,
during inflation, the tachyonic instability of light DM modes can efficiently
source isocurvature modes at second order Chung:2004nh ; Chung:2011xd ;
Ling:2021zlj . In our previous work Garcia:2022vwm , we postponed the detailed
exploration of power spectrum constraints for the direct decay production
mechanism. In this paper, we investigate the constraints on the dark matter
abundance, the Lyman-$\alpha$ forest, and the isocurvature at quadratic order
for this out-of-equilibrium DM production scenario.

This paper is structured as follows: Section II discusses the non-perturbative
gravitational production of scalar dark matter during inflation and reheating.
In Section III, we introduce the phase space distribution of DM and compute
the dark matter relic abundance. Section IV is dedicated to computing the dark
matter isocurvature power spectrum, and in Section V, we combine the relic
abundance, structure formation, and isocurvature constraints. Finally, Section
VI summarizes our findings. In this paper, we use natural units
$k_{B}=\hbar=c=1$ and the metric signature $(+,-,-,-)$.

##  II Inflation and Dark Matter Production

We examine the production of a scalar spectator field, $\chi$, that is
associated with dark matter. We assume that the scalar inflaton field, $\phi$,
can couple to this singlet scalar dark matter particle. For both fields, a
minimal coupling to gravity is assumed111The presence of large non-minimal
couplings to gravity in the inflaton or dark matter sector would significantly
alter the model, and we leave this analysis for future work. and a dimension-4
interaction $\frac{1}{2}\sigma\phi^{2}\chi^{2}$. The scalar field action can
be expressed as

$\displaystyle\mathcal{S}\;=\;\int\mathop{}\\!\mathrm{d}^{4}x\,$
$\displaystyle\sqrt{-g}\Big{[}\frac{1}{2}(\partial_{\mu}\phi)^{2}+\frac{1}{2}(\partial_{\mu}\chi)^{2}-V(\phi)$
$\displaystyle-\frac{1}{2}\left(m_{\chi}^{2}+\sigma\phi^{2}\right)\chi^{2}+\mathcal{L}_{\phi-{\rm
SM}}\Big{]}\,.$ (1)

Here $m_{\chi}$ represents the bare mass of dark matter, $\sigma$ is the
dimensionless inflaton-DM coupling, $\mathcal{L}_{\phi-{\rm SM}}$ is the not
yet specified inflaton-Standard Model interaction, that is responsible for
reheating of the universe, and $V(\phi)$ is the inflaton potential. We
consider two particular models of inflation: the Starobinsky model
Starobinsky:1980te ; Ellis:2013nxa

$V(\phi)\;=\;\frac{3}{2}\lambda
M_{P}^{4}\left[1-e^{-\sqrt{\frac{2}{3}}\frac{\phi}{M_{P}}}\right]^{2}\,,$ (2)

and the T-model Kallosh:2013maa

$V(\phi)\;=\;\lambda
M_{P}^{4}\left[\sqrt{6}\tanh\left(\frac{\phi}{\sqrt{6}M_{P}}\right)\right]^{2}\,,$
(3)

where $M_{P}=1/\sqrt{8\pi G_{N}}\simeq 2.435\times 10^{18}\,{\rm GeV}$ denotes
the reduced Planck mass and the constant $\lambda$ determines the inflaton
mass at its minimum at $\phi=0$, with $m_{\phi}=\sqrt{2\lambda}M_{P}$ for both
models. The value of $\lambda$ is fixed by the measurement of the amplitude of
the primordial curvature power spectrum. The principal CMB observables are
evaluated at the pivot scale $k_{*}=0.05\,{\rm Mpc}^{-1}$ used in the Planck
analysis Planck:2018jri , and they can be expressed in terms of the slow-roll
parameters:

$\displaystyle n_{s}$ $\displaystyle\simeq 1-6\epsilon_{*}+2\eta_{*}\,,\qquad
r\simeq 16\epsilon_{*}\,,$ (4) $\displaystyle A_{S*}$
$\displaystyle\simeq\frac{V_{*}}{24\pi^{2}\epsilon_{*}M_{P}^{4}}\,.$ (5)

Here $V_{*}=V(\phi_{*})$, $n_{s}$ is the scalar tilt, $r$ is the tensor-to-
scalar ratio, $A_{S*}\simeq 2.1\times 10^{-9}$ is the curvature power spectrum
amplitude, and the slow-roll parameters are
$\epsilon=\frac{1}{2}M_{P}^{2}(V_{,\phi}/V)^{2}$ and
$\eta=M_{P}^{2}\left(V_{,\phi\phi}/V\right)$. The normalization of the
potentials (2) and (3) can be approximated as Garcia:2020wiy ; Ellis:2021kad

$\displaystyle\lambda\simeq\frac{12\pi^{2}A_{S*}}{N_{*}^{2}}\;\simeq\;\frac{2.5\times
10^{-7}}{N_{*}^{2}}\leavevmode\nobreak\ (\rm{Starobinsky})\,,$ (6)
$\displaystyle\lambda\simeq\frac{3\pi^{2}A_{S*}}{N_{*}^{2}}\;\simeq\;\frac{6.2\times
10^{-8}}{N_{*}^{2}}\leavevmode\nobreak\ (\rm{T\mbox{-}Model})\,,$ (7)

where $N_{*}$ is the number of $e$-folds after horizon crossing.222Note that
the dimensionless parameter $\lambda$ uniquely determines the full inflaton
potential. For this choice of parameter, the Starobinsky and T-model
potentials coincide at large field values. However, the two potentials differ
at the end of inflation and around the minimum, resulting in different
inflaton masses. For a nominal choice of $N_{*}=55$ $e$-folds, we find
$n_{s}\simeq 0.965\leavevmode\nobreak\ (0.963)$ and $r\simeq
0.004\leavevmode\nobreak\ (0.004)$ for Starobinsky model (T-model) of
inflation, and both models are highly favored by current CMB measurements
Planck:2018jri .333See Ref. Ellis:2021kad  for a dedicated analysis on the
BICEP/Keck/WMAP/Planck data constraints on the Starobinsky and T-models of
inflation.

During inflation, the coupling of the dark matter scalar to gravity and the
inflaton results in the non-adiabatic variation of the dark matter effective
mass, leading to subsequent particle production Parker:1968mv ; Parker:1969au
; Parker:1971pt ; Ford:1986sy . After the end of inflation, the scalar
inflaton field, $\phi$, undergoes coherent oscillations about the minimum of
its potential, and the effective mass of dark matter oscillates accordingly.
This leads to additional post-inflationary production of dark matter.444See
Refs. Opferkuch:2019zbd ; Mambrini:2021zpp ; Haque:2021mab ; Clery:2021bwz ;
Clery:2022wib ; Barman:2021ugy ; Barman:2022tzk ; Haque:2022kez  for the
treatment of gravitational production of dark matter during reheating. As the
expansion rate of the universe rapidly slows down at this time, the inflaton
decays into elementary particles, reheating the universe.

Assuming a vanishing vacuum expectation value (VEV) for $\chi$, the population
of the dark matter relic abundance will be sourced by the growth in the $\chi$
fluctuation modes. For convenience, we use conformal time $\tau$, which is
related to cosmic time $t$ as
$\mathop{}\\!\mathrm{d}t/\mathop{}\\!\mathrm{d}\tau=a$, where $a$ is the scale
factor of the flat Friedmann-Robertson-Walker spacetime (FRW) metric. We
introduce the rescaled dark matter field

$X\;\equiv\;a\chi\,,$ (8)

and variation of the action (1) with respect to the DM field $\chi$ leads to
the following equation of motion

$\left[\frac{\mathop{}\\!\mathrm{d}^{2}}{\mathop{}\\!\mathrm{d}\tau^{2}}-\nabla^{2}-\frac{a^{\prime\prime}}{a}+a^{2}m_{\chi}^{2}+\sigma
a^{2}\phi^{2}\right]X(\tau,\boldsymbol{x})\;=\;0\,,$ (9)

where prime denotes the derivative with respect to the conformal time, $\tau$.
The momentum modes of the field $X$ are determined from its canonically
quantized form

$X(\tau,\boldsymbol{x})=\int\frac{\mathop{}\\!\mathrm{d}^{3}\boldsymbol{k}}{(2\pi)^{3/2}}\,e^{-i\boldsymbol{k}\cdot\boldsymbol{x}}\left[X_{k}(\tau)\hat{a}_{\boldsymbol{k}}+X_{k}^{*}(\tau)\hat{a}^{\dagger}_{-\boldsymbol{k}}\right]\,,$
(10)

where $\boldsymbol{k}$ corresponds to the comoving momentum,
$k=|\boldsymbol{k}|$, and $\hat{a}_{\boldsymbol{k}}$ and
$\hat{a}^{\dagger}_{\boldsymbol{k}}$ are the annihilation and creation
operators, respectively, satisfying the canonical commutation relations
$[\hat{a}_{\boldsymbol{k}},\hat{a}^{\dagger}_{\boldsymbol{k}^{\prime}}]=\delta(\boldsymbol{k}-\boldsymbol{k}^{\prime})$,
$[\hat{a}_{\boldsymbol{k}},\hat{a}_{\boldsymbol{k}^{\prime}}]=[\hat{a}^{\dagger}_{\boldsymbol{k}},\hat{a}^{\dagger}_{\boldsymbol{k}^{\prime}}]=0$.
The canonical commutation relations between the field, $X_{k}$, and its
momentum conjugate, $X_{k}^{\prime}$, are fulfilled by the Wronskian
constraint $X_{k}X^{*\prime}_{k}-X_{k}^{*}X_{k}^{\prime}\;=\;i$. Substituting
Eq. (10) into Eq. (9), we obtain the equation of motion for the mode
functions, $X_{k}$,

$X_{k}^{\prime\prime}+\omega_{k}^{2}X_{k}\;=\;0\,,$ (11)

with the mode frequency given by

$\displaystyle\omega_{k}^{2}\;$
$\displaystyle\equiv\;k^{2}-\frac{a^{\prime\prime}}{a}+a^{2}m_{\chi}^{2}+\sigma
a^{2}\phi^{2}\leavevmode\nobreak\ $
$\displaystyle=\;k^{2}+a^{2}\left(\frac{R}{6}+m_{\chi}^{2}+\sigma\phi^{2}\right)\,,$
(12)

where $R$ denotes the Ricci scalar. The initial conditions for $X_{k}$ are
chosen to be the positive-frequency Bunch-Davies vacuum

$X_{k}(\tau_{0})\;=\;\frac{1}{\sqrt{2\omega_{k}}}\,,\
X_{k}^{\prime}(\tau_{0})\;=\;-\frac{i\omega_{k}}{\sqrt{2\omega_{k}}}\,,$ (13)

which are valid for modes that are deep inside the Hubble horizon, with
$|\tau_{0}\omega_{k}|\gg 1$. The efficiency of DM particle production is
determined by the ratio of couplings $\sigma/\lambda$ Garcia:2021iag ;
Garcia:2022vwm . Eq. (II) can be written equivalently as

$\displaystyle\omega_{k}^{2}$ $\displaystyle=k^{2}+a^{2}m_{\rm{eff}}^{2}\,,$
(14) $\displaystyle{\rm{with}}\leavevmode\nobreak\ \leavevmode\nobreak\
m_{\rm{eff}}^{2}$ $\displaystyle\equiv
m_{\chi}^{2}+\sigma\phi^{2}-\left(\frac{1-3w}{2}\right)H^{2}\,.$ (15)

Here $H=\dot{a}/a$ is the Hubble parameter and $w\equiv P/\rho$ is the
equation of state parameter, where $P$ and $\rho$ are the background pressure
and energy density, respectively. For low momentum (IR) modes, the mode
frequency squared may become negative during inflation if
$k^{2}<a^{2}|m_{\rm{eff}}^{2}|$, leading to a fast mode growth during
inflation due to the tachyonic instability. If the DM bare mass is negligible,
$m_{\chi}\ll H$, this occurs when

$k^{2}\lesssim a^{2}\left(\frac{1-3w}{2}H^{2}-\sigma\phi^{2}\right)\,,$ (16)

which implies the tachyonic growth for the low momentum (IR) modes when

$\frac{\sigma}{\lambda}\;\lesssim\;\frac{\alpha(1-3w)}{2}\left(\frac{M_{P}}{\phi}\right)^{2}\leavevmode\nobreak\
{\rm{with}}\leavevmode\nobreak\ \leavevmode\nobreak\
\alpha\equiv\frac{H^{2}}{\lambda M_{P}^{2}}\,.$ (17)

During inflation, when $t\ll t_{\rm{end}}$, $\alpha\simeq
0.5\leavevmode\nobreak\ (2)$ for the Starobinsky model (T-model); at the end
of inflation, when $t=t_{\rm{end}}$, $\alpha\simeq 0.12\leavevmode\nobreak\
(0.33)$, which implies the tachyonic growth when $\sigma/\lambda\lesssim
0.5\leavevmode\nobreak\ (0.3)$ for the Starobinsky model (T-model). Therefore,
for a smaller effective coupling, the DM phase space distribution is highly
sensitive to the inflationary dynamics.

If the dimensionless inflaton-DM coupling is zero, $\sigma=0$, then the scalar
dark matter is produced solely through gravitational interactions. We find
from Eq. (11) that tachyonic mode growth occurs when

$k^{2}\lesssim a^{2}\left(\frac{1-3w}{2}H^{2}-m_{\chi}^{2}\right)\,.$ (18)

During inflation, $w\simeq-1$, and for low momentum (IR) modes, the tachyonic
mode enhancement arises when $m_{\chi}^{2}\lesssim 2H^{2}$.

It is important to note that the growth of scalar fluctuations for scales that
exit the horizon 50-60 $e$-folds before the end of inflation can efficiently
generate isocurvature perturbations that would leave an imprint in the CMB. In
Sections IV and V, we provide a detailed discussion and calculation of dark
matter isocurvature and constrain the associated parameter space. Our results
are largely independent of the specific model of inflation and can be applied
to a variety of inflationary scenarios.

Figure 1: Phase space distribution of the dark matter field $\chi$ for a range
of effective couplings $\sigma/\lambda$ and bare masses $m_{\chi}$ for the
T-model of inflation. Left panel: the PSD dependence on the inflaton-DM
coupling for a fixed DM mass $m_{\chi}=10^{-2}H_{\rm{end}}$. As the coupling
ratio $\sigma/\lambda$ increases, the long wavelength (IR) modes become more
suppressed, and dominant particle production occurs towards the end of
inflation. Right panel: the PSD dependence on the DM mass for a fixed
inflaton-DM coupling $\sigma/\lambda=10^{-2}$. In the IR regime,
$f_{\chi}\propto q^{-2.9}$, and in the UV regime, $f_{\chi}\propto q^{-9/2}$.
In both panels, vertical lines are displayed corresponding to the rescaled
comoving momenta of the present-day horizon scale, $q_{0}\simeq 5.8\times
10^{-27}$, and the Planck pivot scale, $q_{*}\simeq 1.3\times 10^{-24}$,
assuming $N_{*}=55$. For reference, the horizon scale at the end of inflation
is $q_{\rm end}\simeq 0.4$.

##  III Dark matter phase space distribution

To extract the DM phase space distribution, $f_{\chi}(k,t)$, we numerically
solve the equation of motion (11) for a wide range of momentum modes. This
distribution is determined by the comoving particle occupation number of the
scalar field, $n_{k}$, which can be expressed in a UV-finite form
Kofman:1997yn ; Garcia:2021iag 555This expression can also be equivalently
expressed as
$n_{k}=\frac{\omega_{k}}{2}(|X_{k}^{\prime}|^{2}/\omega_{k}^{2}+|X_{k}|^{2})-\frac{1}{2}$
Kofman:1997yn .

$f_{\chi}(k,t)\;=\;n_{k}\;=\;\frac{1}{2\omega_{k}}\left|\omega_{k}X_{k}-iX^{\prime}_{k}\right|^{2}\,.$
(19)

This quantity can only be interpreted as the occupation number for non-
tachyonic modes, and therefore it must be evaluated well after the end of
inflation. For simplicity, we introduce a dimensionless comoving momentum
rescaled at the end of inflation with respect to the inflaton mass,
$\boldsymbol{q}=\boldsymbol{k}/(m_{\phi}a_{\rm{end}})$. The comoving number
density of the produced DM particles can then be expressed as Garcia:2022vwm

$n_{\chi}\left(\frac{a}{a_{\rm
end}}\right)^{3}\;=\;\frac{m_{\phi}^{3}}{2\pi^{2}}\int\mathop{}\\!\mathrm{d}q\,q^{2}f_{\chi}(q,t)\,.$
(20)

The inflaton mass is given by $m_{\phi}=\sqrt{2\lambda}M_{P}\simeq
4.1H_{\rm{end}}\simeq 3\times 10^{13}\,\rm{GeV}$ for the Starobinsky model and
$m_{\phi}=\sqrt{2\lambda}M_{P}=2.5H_{\rm{end}}\simeq 1.6\times
10^{13}\,\rm{GeV}$ for the T-model.

For concreteness, we assume that the Planck pivot scale $k_{*}$ leaves the
horizon $N_{*}=55$ $e$-folds before the end of inflation and a total duration
of inflation of $N_{\rm tot}=76.5$ $e$-folds. We have chosen this value to
ensure that the present comoving scale was deeply inside the Hubble horizon at
the beginning of inflation, without extending the period of accelerated
expansion too far back in time. While the duration of inflation does not
affect the form of the PSD due to the linearity of the equation of motion for
the momentum modes (11), integrated quantities, such as $n_{\chi}$, generally
depend on the infrared (IR) cutoff in the spectrum, which we discuss in detail
below.

We employ numerical techniques developed in Ref. Ellis:2014opa  to compute the
power spectra in multifield inflation, to obtain the complete PSD down to
wavenumbers below the CMB pivot scale without the need for extrapolation. This
represents an improvement over our previous analysis in Garcia:2022vwm . The
integrand of Eq. (20) is shown in Fig. 1 for a range of DM masses and
couplings. Although we present the results for the T-model of inflation, they
can be translated to other inflationary models by identifying different values
of Hubble parameter at the end of inflation between the models.666This
translation is less relevant for the isocurvature constraints since they
primarily depend on $H_{*}$. We discuss it in detail in Section IV. To
translate the constraints from the T-model to the Starobinsky model, we can
use the relationship

$\left[\sigma\left(\frac{M_{P}}{H_{\rm end}}\right)^{2}\right]_{\rm
T\mbox{-}model}\;=\;\left[\sigma\left(\frac{M_{P}}{H_{\rm
end}}\right)^{2}\right]_{\rm Staro.}\,,$ (21)

which, in terms of the effective dimensionless coupling ratio, corresponds to
$(\sigma/\lambda)_{\rm T\mbox{-}model}\simeq 2.8(\sigma/\lambda)_{\rm
Staro.}$. For the purely gravitational regime, when $\sigma=0$, the
constraints can be directly related by transforming between different values
of $H_{\rm{end}}$.777Note that, for our choice of normalization for the
Starobinsky and T-model potentials, the inflaton mass is different in both
cases, and the comoving momentum $q$ has to be rescaled accordingly.

The left panel of Fig. 1 displays the phase space distributions for a fixed
mass $m_{\chi}=10^{-2}H_{\rm end}$ for several values of the effective
coupling constant ratio $\sigma/\lambda\ll 1$. For the smallest couplings, the
spectrum presents a steep red tilt, corresponding to $f_{\chi}\propto q^{-3}$
in the pure gravitational case ($\sigma=0$). This is consistent with the
$k$-dependence of the squared mode function from Eq. (44), as derived in
Appendix B, which is devoted to analytical approximations. For a comprehensive
treatment of purely gravitational dark matter production, see also Ref.
Garcia:2022vwm ; Kaneta:2022gug ; Ford:2021syk .

The left panel of Fig. 1 also illustrates the suppression of the growth of
long wavelength (IR) modes for $\sigma\neq 0$. This suppression is more
significant for modes that exit the horizon earlier because the inflaton field
value during inflation is larger, and the tachyonic enhancement is weaker (see
Eq. (17)). This observation is qualitatively consistent with the analytical
approximations provided in App. B. Couplings $\sigma/\lambda\gtrsim 10^{-1}$
are not IR-sensitive, and result in particle production that primarily occurs
towards the end of inflation and during reheating. As discussed in App. B, the
phase space distribution in the long wavelength (IR) regime scales as
$f_{\chi}\propto q^{-2\nu}$, where $\nu=\sqrt{9/4-m_{\rm{eff}}^{2}/H^{2}}$ for
real $\nu$, and reduces to $f_{\chi}\propto q^{-3}$ when $m_{\rm{eff}}\ll H$,
which corresponds to $\sigma\langle\phi\rangle^{2}\ll m_{\chi}^{2}$.

For a non-vanishing coupling $\sigma$, the number density $n_{\chi}$ converges
over the full range of momenta $0\leq q<\infty$. However, if we assume a
finite duration of inflation, the integrated PSD is sensitive to the IR
endpoint for $\sigma/\lambda\lesssim 0.1$, which may be associated with the
mode that exited the horizon at the beginning of inflation, or more
conservatively, with the wavenumber of the present horizon scale. We choose to
ignore the contribution of larger scales that have never entered our causal
radius and simply renormalize the background Ling:2021zlj ; Herring:2019hbe ;
Starobinsky:1994bd ; Felder:1999wt .

In all cases, the $q>1$ UV tail of the PSD corresponds to modes excited during
the post-inflationary reheating epoch. In the small coupling regime, one can
use the perturbative Boltzmann formalism to study how the energy of the
inflaton dissipates while it undergoes coherent oscillations around its
minimum. The UV tail of the distribution is given by Garcia:2022vwm ;
Kaneta:2022gug

$\displaystyle f_{\chi}(q,t)$
$\displaystyle\simeq\frac{\sqrt{3}\pi(\sigma-\lambda)^{2}\rho_{\phi}^{3/2}(t_{\rm
end})M_{P}}{16m_{\phi}^{7}}$ $\displaystyle\times
q^{-9/2}e^{-1.56\left(\frac{a_{\rm end}}{a_{\rm
reh}}\right)^{2}q^{2}}\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ (q>1)\,,$ (22)

where $\rho_{\phi}(t_{\rm{end}})$ is the inflaton energy density at the end of
inflation and $a_{\rm{reh}}$ is the scale factor when reheating occurs.888We
define the time of reheating as the moment when the energy density of
radiation starts to dominate the energy budget of the universe.

The right panel of Fig. 1 demonstrates the impact of varying the bare mass of
the DM field for a fixed inflaton-DM coupling of $\sigma/\lambda=10^{-2}$,
near the purely gravitational particle production regime. We limit our focus
to light bare masses, $10^{-3}<m_{\chi}/H_{\rm end}<1$, and the results are
extrapolated to smaller masses. Following the end of inflation, the dark
matter scalar is not yet fully decoupled, and the PSD relaxes to its final
form as reheating proceeds. This relaxation is modulated by the oscillations
of the inflaton, and is dependent on the dark matter mass $m_{\chi}$, with
lighter DM asymptoting to the final value at a smaller rate. This need for a
longer evaluation time is in conflict with the increased precision needed to
enforce the Wronskian constraint for $X_{k}$, and better results are achieved
with extrapolation (see Ref. Garcia:2022vwm  for further details).

##  IV Dark matter isocurvature

In the previous section, we numerically computed the PSD for the scalar dark
matter field $\chi$. The results show that the PSD is red-tilted and has a
large amplitude at CMB scales if $\sigma/\lambda<10^{-1}$. This efficient mode
excitation during inflation is expected to affect the primordial power
spectrum of scalar fluctuations. We note that we assume a zero mean vacuum
expectation value for the DM scalar field ($\langle\chi\rangle=0$, but
$\langle\chi^{2}\rangle\neq 0$).999See Refs. Ema3 ; Tenkanen:2019aij  for
discussions on isocurvature constraints assuming a non-zero VEV of the dark
matter scalar field. For a vanishing background value of $\chi$, the
fluctuations of the dark matter field will not directly source the curvature
perturbation. Therefore, they can be considered to be pure isocurvature
fluctuations in the comoving gauge Chung:2011xd ; Chung:2004nh . Isocurvature
fluctuations have not yet been detected at CMB scales, which imposes an upper
bound on the primordial power spectrum of
$\mathcal{P}_{\mathcal{S}}(k_{*})\lesssim 8.3\times 10^{-11}$ Planck:2018jri .

Figure 2: DM isocurvature power spectrum for different inflaton-DM couplings
with $m_{\chi}/H_{\rm end}=10^{-2}$, with each coupling represented by a
different color. The vertical lines indicate the present horizon scale and the
Planck pivot scale, as shown in Fig. 1. Here we set $N_{\rm tot}=76.5$. The
left panel corresponds to the T-model, while the right panel corresponds to
the Starobinsky model, where the DM coupling $\sigma$ has been rescaled using
Eq. (21).

The suppression of the linear isocurvature power spectrum results from the
geodesic evolution of the background in the scalar field manifold. However,
the growth in the energy density of DM is dominated by the quadratic
fluctuations that source the variance $\langle\chi^{2}\rangle$ Chung:2004nh ;
Chung:2015pga ; Ling:2021zlj ; Redi:2022zkt . The second-order contribution to
the isocurvature power spectrum can be approximated using equation
Liddle:1999pr ; Chung:2004nh ; Ling:2021zlj

$\mathcal{P}_{\mathcal{S}}(k)\;=\;\frac{k^{3}}{2\pi^{2}\rho_{\chi}^{2}}\int\mathop{}\\!\mathrm{d}^{3}\boldsymbol{x}\
\langle\delta\rho_{\chi}(\boldsymbol{x})\delta\rho_{\chi}(0)\rangle
e^{-i\boldsymbol{k}\cdot\boldsymbol{x}}\,,$ (23)

where $\rho_{\chi}$ and $\delta\rho_{\chi}$ denote the DM energy density and
its fluctuation, respectively. Upon Fourier transformation, the spectrum can
be evaluated in terms of the $\chi$ mode functions and their (conformal) time
derivatives Ling:2021zlj . The details of the evaluation of the integral in
equation (23) are provided in Appendix A. Like the PSD, the late-time
isocurvature power spectrum is to be evaluated after DM decoupling. For
$10^{-3}\lesssim m_{\chi}/H_{\rm end}$, it is sufficient to evolve the mode
functions $X_{k}$ a few ($\sim 5-6$) $e$-folds after the end of inflation to
obtain a good approximation of the final asymptotic value of
$\mathcal{P}_{\mathcal{S}}$.

Numerical approach. Following the method detailed in Appendix A, we compute
the isocurvature power spectrum for various DM masses and inflaton-DM
couplings, down to CMB scales and below. Fig. 2 shows the shape of the
isocurvature power spectrum for a broad range of comoving scales with fixed
$m_{\chi}/H_{\rm end}=10^{-2}$ and varying coupling ratio $\sigma/\lambda$.
The scales are normalized to the comoving momentum of the mode that exited the
horizon at the end of inflation. The figure demonstrates that a weak coupling
between the inflaton and the dark sector produces a nearly scale-invariant
isocurvature spectrum, whose amplitude is in strong disagreement with the
Planck constraint. As the coupling strength increases, the amplitude is
suppressed for long wavelength modes, which is consistent with a more
efficient particle production at the onset of reheating.

Fig. 3 illustrates the dependence of $\mathcal{P}_{\mathcal{S}}(k_{*})$ on the
bare DM mass $m_{\chi}$ for three effective couplings $\sigma/\lambda$. We
consider the T-model of inflation in this example, but identical constraints
can be obtained for the Starobinsky model using the rescaling relationship
(21). Importantly, our results are not limited to these specific models of
inflation, and can be applied to different models as well by introducing a
similar rescaling relationship. As expected, the amplitude of the isocurvature
spectrum is smaller for heavy DM since the tachyonic enhancement is weaker
(see Eq. (14)). Conversely, for $m_{\chi}\ll H_{\rm end}$, a clear red-tilted
power-law trend can be observed. Therefore, for a given coupling, we can
extrapolate to smaller masses to obtain the full constraints on the minimally
coupled scalar DM scenario. We find numerically that for $N_{*}=55$, the
isocurvature constraints are satisfied in the purely gravitational regime
when101010We note that the value of the Hubble parameter at the end of
inflation, $H_{\rm{end}}$, depends on the specific inflationary potential
being considered.

$\displaystyle m_{\chi}\gtrsim
1.11\,H_{\rm{end}}\qquad(\text{Starobinsky})\,,\leavevmode\nobreak\ $
$\displaystyle m_{\chi}\gtrsim 1.34\,H_{\rm{end}}\qquad(\text{T-model})\,,$
(24)

and in the small bare mass limit $m_{\chi}/H_{\rm end}\ll 1$, we find

$\displaystyle\frac{\sigma}{\lambda}\;\gtrsim\;0.008\qquad(\text{Starobinsky})\,,\leavevmode\nobreak\
$
$\displaystyle\frac{\sigma}{\lambda}\;\gtrsim\;0.02\qquad(\text{T-model})\,.$
(25)

ignoring a weak DM mass dependence (see Fig. 5 below). The isocurvature
constraints are also weakly sensitive to the value of $N_{*}$, which in turn
implies that they are insensitive to the reheating mechanism. Therefore, these
constraints can be expressed in terms of the Hubble parameter at horizon
crossing, $H_{*}$, given by Eqs. (29) and (31).

Analytical approach. Appendix B explores an analytical derivation of the
isocurvature constraint based on the evaluation of the mode functions after
inflation. We approximate this evaluation using a finite-duration pure de
Sitter era. For the purely gravitational case ($\sigma=0)$, we find that the
limit for the Starobinsky model (T-model) is $m_{\chi}\gtrsim
1.18H_{\rm{end}}\,(1.41H_{\rm{end}})$, and for negligible bare mass, we find
the limit $\sigma/\lambda\gtrsim 0.006\,(0.02)$ for the Starobinsky model
(T-model) (see Eqs. (B) and (B) and their corresponding derivation). These
analytical results are in good agreement with the fully numerical constraints
of Eqs. (IV) and (IV).

However, small discrepancies arise due to several factors. Firstly, some
parameters that we approximate as constant actually vary with time as the
dynamics of the universe departs from a pure de Sitter era due to slow-roll
corrections. In particular, we neglected the explicit time derivative terms
for the mode functions in the integrand of the isocurvature expression of Eq.
(41). Additionally, the field excursion during inflation is quite large, with
$\Delta\phi/M_{p}\equiv(\phi_{*}-\phi_{\text{end}})/M_{p}\sim 4.7\,(5.3)$ for
the Starobinsky model (T-model), and the end of inflation deviates
significantly from a pure de Sitter phase. Typically, for couplings
$\sigma/\lambda\sim\mathcal{O}(10^{-2})$, we find that the parameter
$\nu^{2}=9/4-m_{\rm{eff}}^{2}/H^{2}$ changes sign between the beginning and
end of inflation. This means that the frequency for IR modes switches from a
purely imaginary to a real quantity during inflation, making any estimate
based on a constant $\nu$ inaccurate. Furthermore, the mode functions continue
to evolve after the end of inflation, typically until $a/a_{\text{end}}\simeq
20-30$, as the universe transitions from a quasi-de Sitter to a matter
domination phase.

Infrared sensitivity. As discussed in Appendix B, the analytical approach
allows us to identify the dependence of our results on the infrared cutoff,
$k_{\text{IR}}$. From Eq. (50), we observe that the isocurvature power
spectrum has a logarithmic dependence on the IR cutoff, given by
$\sim\log(k_{*}/k_{\text{IR}})$. However, since it is also exponentially
sensitive to the ratio of the effective dark matter mass with respect to the
Hubble scale during inflation, changing $k_{\text{IR}}$ by orders of magnitude
would only result in a small correction of $\mathcal{O}(1)$ to the
isocurvature bound. Therefore, our results are robust and rather insensitive
to the infrared cutoff.

Figure 3: DM isocurvature power spectrum at the Planck pivot scale
$k_{*}=0.05\,{\rm Mpc}^{-1}$ with $N_{*}=55$ $e$-folds as a function of the DM
mass $m_{\chi}/H_{\rm{end}}$. We show three different inflaton-DM effective
couplings $\sigma/\lambda=0,10^{-2}$, and $10^{-3/2}$.

##  V Relic abundance, structure formation, and isocurvature constraints

We now discuss the implications of the computed phase space distribution (PSD)
and isocurvature power spectrum on the model parameter space. First, we must
ensure that the measured DM relic abundance, with closure fraction
$\Omega_{\rm DM}h^{2}=0.1198$ Planck:2018vyg , can be obtained at some
reheating temperature $T_{\rm BBN}\leq T_{\rm reh}\leq m_{\phi}$ with $T_{\rm
BBN}\sim 1\leavevmode\nobreak\ \text{MeV}$. With the PSD determined, the relic
density can be computed as Garcia:2022vwm

$\Omega_{\rm DM}\;\simeq\;\frac{m_{\chi}n_{\chi}}{\rho_{c}}=\;\frac{1}{6\pi
q_{0}^{3}}\left(\frac{m_{\chi}H_{0}}{M_{P}^{2}}\right)\int_{q_{0}}^{\infty}\mathop{}\\!\mathrm{d}q\,q^{2}f_{\chi}(q)\,.$
(26)

Here $H_{0}=100\,h\,\rm{km\,s^{-1}\,Mpc^{-1}}$ is the present Hubble
parameter, $\rho_{c}=1.05\times 10^{-5}\,h^{2}\,\rm{GeV\,cm^{-3}}$ is the
present critical energy density, with $h\sim 0.67$ Planck:2018vyg , and we
integrate only down to the present comoving scale $q_{0}$. For the coupling
ratio $\sigma/\lambda\ll 1$, saturation of the DM relic density leads to a
functional relation between the reheating temperature and the DM mass. This
connection was studied in detail in Garcia:2022vwm , and we reproduce the
results in Fig. 4. For pure gravitational production, the relic density is
mass-independent, with $\Omega_{\rm DM}h^{2}\simeq 0.12(T_{\rm reh}/34\,{\rm
GeV})$ Ling:2021zlj ; Garcia:2022vwm . As the inflaton-DM coupling is
increased, a functional relation between $T_{\rm reh}$ and $m_{\chi}$ arises,
given by $(T_{\rm reh}/1\,{\rm GeV})\propto(m_{\chi}/1\,{\rm GeV})^{\gamma}$,
with $-1\leq\gamma\leq 0$. For $\sigma/\lambda\gtrsim 10^{-1/2}$, $\Omega_{\rm
DM}\propto T_{\rm reh}m_{\chi}$, corresponding to $\gamma=-1$. Further details
can be found in Garcia:2022vwm .

Figure 4: DM masses and reheating temperatures for which the measured DM relic
abundance is saturated, for $\sigma/\lambda\ll 1$. The gray shaded region
corresponds to the exclusion of light DM masses from the Lyman-$\alpha$
measurement of the matter power spectrum.

For $\sigma/\lambda\lesssim 10^{-1}$, the DM relic abundance can be saturated
with very small DM masses. However, if the DM particles are very light, they
could possess a non-negligible pressure component and depart substantially
from the cold dark matter paradigm. This departure is manifested as a
suppression in the matter power spectrum and the corresponding erasure of
clustered structure overdensities. This cutoff in the spectrum is not present
down to scales $k_{\text{H}}(a=1)\simeq 15\,h\,\text{Mpc}^{-1}$, as determined
by the Lyman-$\alpha$ forest measurement of the distribution of matter. For a
thermalized DM relic particle (warm dark matter, WDM), this constraint is
presented as a lower bound on its mass,
$m_{\text{WDM}}>m_{\text{WDM}}^{\text{Ly}\mbox{-}\alpha}$ Narayanan:2000tp ;
Viel:2005qj ; Viel:2013fqw ; Baur:2015jsy ; Irsic:2017ixq ; Palanque-
Delabrouille:2019iyz ; Garzilli:2019qki .

For out-of-equilibrium DM production mechanisms, the Ly-$\alpha$ bound is
dependent on the details of the production and decoupling of the dark
particles. Notably, one can translate the WDM lower bound to a bound on the
not-equilibrated DM by matching their equation of state parameters. This
matching procedure yields the lower bound Ballesteros:2020adh

$\,m_{\chi}^{\text{Ly}\mbox{-}\alpha}\;=\;m_{\rm
WDM}^{\text{Ly}\mbox{-}\alpha}\left(\frac{T_{\star}}{T_{\rm
WDM,0}}\right)\sqrt{\frac{\langle q^{2}\rangle}{\langle q^{2}\rangle_{\rm
WDM}}}\,,$ (27)

where $T_{\rm WDM,0}$ is the WDM temperature, $T_{\star}$ is the
characteristic energy scale of the produced DM, $T_{\star}=m_{\phi}(a_{\rm
end}/a_{0})$ for the model discussed in this work, and

$\langle
q^{2}\rangle\,\equiv\,\dfrac{\int\mathop{}\\!\mathrm{d}q\,q^{4}f_{\chi}(q)}{\int\mathop{}\\!\mathrm{d}q\,q^{2}f_{\chi}(q)}\,.$
(28)

Given the PSD for $\chi$, displayed in Fig. 1, we can readily examine the
structure formation constraint. However, for weak inflaton-DM coupling,
$\langle q^{2}\rangle$ is dependent on $m_{\chi}$. Therefore, Eq. (27) is a
non-linear equation in $m_{\chi}$ whose solution provides the lower bound
$m_{\chi}^{\text{Ly}\mbox{-}\alpha}$ Garcia:2022vwm . The Lyman-$\alpha$
constraint is shown in Fig. 4, leading to the exclusion of masses smaller than
$m_{\chi}\simeq 0.34\,{\rm meV}$ for pure gravitational production. Note that
for very small couplings, the lower bound on the DM mass stems from the
Ly-$\alpha$ result, whereas for the largest coupling, it is determined by the
saturation of $\Omega_{\rm DM}$.

Figure 5: Allowed regions of scalar DM mass and its coupling to the inflaton,
constrained by the relic abundance, structure formation, and isocurvature
perturbations, assuming perturbative reheating that satisfies Big Bang
Nucleosynthesis constraints ($T_{\rm BBN}<T_{\rm reh}<m_{\phi}$). The green
strip of the boundary of the isocurvature constraint corresponds to varying
number of $e$-folds $47\leq N_{*}\leq 55$. The blue shaded region corresponds
to a region where the power spectrum is sensitive to the DM mass for
$m_{\chi}\ll H_{\rm end}$, while the orange region corresponds to particle
production via broad resonance, resulting in a fragmented inflaton field (see
Garcia:2022vwm  for details).

The last constraint to be satisfied, and the main result of the present work,
corresponds to the bound on the parameter $\beta_{\rm iso}$. This constraint
depends on the mass of the dark matter particle $m_{\chi}$, the coupling ratio
$\sigma/\lambda$, and the ratio $k_{*}/k_{\rm end}$ (see Figs. 2 and 3). Fig.
3 illustrates that for $\sigma/\lambda<10^{-2}$ and $m_{\chi}\ll H_{\rm end}$,
only super-heavy scalar dark matter can avoid the Planck bound on $\beta_{\rm
iso}$ Ling:2021zlj . For larger couplings, the allowed parameter space is
determined by simultaneously satisfying the constraints on the relic
abundance, structure formation, and isocurvature perturbations, and we show it
for T-model in Fig. 5. The gray shaded region is excluded by the Ly-$\alpha$
bound and the requirement that the dark matter density parameter $\Omega_{\rm
DM}h^{2}\simeq 0.12$, subject to the additional constraint $T_{\rm BBN}<T_{\rm
reh}<m_{\phi}$. This assumption is necessary to ensure that the post-
inflationary reheating can be described by perturbation theory, without
conflicting with the predictions of BBN. The upper bound on $T_{\rm reh}$
determines the shape of the lower exclusion region for
$10^{-1}\lesssim\sigma/\lambda<10^{4}$. The lower bound on $T_{\rm reh}$ in
turn controls the shape of the upper exclusion region at large couplings.

The light green region in Fig. 5 corresponds to the masses and couplings
excluded by the upper bound on $\mathcal{P}_{\mathcal{S}}(k_{*})$ for $N_{*}$.
As discussed above, the constraints cannot be satisfied for weak couplings
when $m_{\chi}\ll H_{\rm end}$, particularly for pure gravitational
production, which is therefore completely excluded. We note the weak
dependence of the exclusion region on $m_{\chi}$, which is primarily
controlled by the inflaton-DM coupling strength. The combination of
Ly-$\alpha$ and $\mathcal{P}_{\mathcal{S}}(k_{*})$ limits result in a lightest
allowed DM mass $m_{\chi}\simeq 5\,{\rm meV}$ for $N_{*}=55$ and
$\sigma/\lambda\simeq 0.024$.

It is worth noting that the green strip in the forbidden domain has a sizable
width. This width corresponds to the dependence of the isocurvature spectrum
at the Planck pivot scale on the post-inflationary dynamics, specifically, the
relation between $k_{*}$ and $N_{*}$. For T-model of inflation, the Planck
constraint on the curvature spectral tilt $n_{s}$ and the tensor-to-scalar
ratio $r$ at the $2\sigma$ level imposes the lower bound $N_{*}\gtrsim 47$ for
the number of $e$-folds after horizon crossing of the pivot scale. On the
other hand, $N_{*}=55$ is close to the perturbative limit Ellis:2021kad . The
width of the green strip corresponds then to the saturation of $\beta_{\rm
iso}$ between 55 and 47 $e$-folds, with the latter corresponding to the
rightmost part of this boundary.

##  VI Discussion and Conclusions

The hypothesis of an independent population of the visible and invisible
sectors of the universe, arising from the decay of the inflaton field during
and after inflation, is not only consistent with the lack of a direct
detection signal for DM in the absence of a non-gravitational direct
interaction between the Standard Model and dark matter, but also a plausible
scenario. In our previous study Garcia:2022vwm , we determined the range of
inflaton-DM couplings and DM masses for which the measured DM relic abundance
can be produced from inflaton decay, that accounts for the observed DM relic
abundance from inflaton decay, while also avoiding the constraint associated
with the non-observation of structure formation suppression at small scales.
In this study, we further narrow down the available parameter space by taking
into account the observation of pure adiabatic fluctuations at CMB scales.

Pure gravitational production. In the absence of a direct inflaton-dark matter
coupling, i.e. $\sigma=0$, the isocurvature power spectrum is nearly scale-
invariant and must satisfy the constraint
$\mathcal{P}_{\mathcal{S}}(k_{*})\lesssim 8.3\times 10^{-11}$ at the Planck
pivot scale. The constraint (IV) can be directly translated to the following
bound

$m_{\chi}\,\gtrsim\,0.54\,H_{*}\,,\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ [\sigma=0]$ (29)

where $H_{*}$ is the Hubble parameter evaluated at the CMB crossing scale. Our
results are consistent with previous studies in the pure gravitational regime
for $m_{\chi}<H_{*}$ Chung:2004nh ; Chung:2011xd ; Ling:2021zlj ; Redi:2022zkt
. If we apply this result to the Starobinsky model (T-model), we find the
constraint

$m_{\chi}\,>\,8.5\,(8.4)\times 10^{12}\leavevmode\nobreak\
\text{GeV}\,.\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ [\sigma=0]\,.$
(30)

Light mass limit. Furthermore, we explore the form and amplitude of the
spectrum in full numerical detail for non-vanishing couplings. In the limit
where the bare dark matter mass is much smaller than the Hubble scale,
$m_{\chi}\ll H$, we found that for small couplings, the moderation of the
tachyonic growth of DM due to the large field value of the inflaton is
inadequate to effectively suppress the isocurvature spectrum for long-
wavelength modes. This results in a constraint on the direct coupling given by

$\sigma\gtrsim(0.01-0.02)\left(\frac{H_{*}}{M_{P}}\right)^{2}\,,$ (31)

where the range comes from a weak $m_{\chi}$ and $N_{*}$-dependence. In terms
of the dimensionless ratio $\sigma/\lambda$, for the Starobinsky and T-models,
this bound is given by Eq. (IV). It is worth nothing that these constraints
still leave open a range of small masses/couplings where scalar DM can be as
light as $5\,{\rm meV}$ with couplings of
$\sigma/\lambda\sim\mathcal{O}(10^{-1})$. Smaller masses and/or couplings
might be achievable in scenarios with non-minimal conformal coupling, and we
intend to explore such models in detail in future work.

Generality and scope of our results. We have introduced an infrared cutoff
$k_{\text{IR}}$ that selects only modes corresponding to wavelengths of the
order of maximal causal distances at present. Our numerical results,
supplemented by analytical estimates, have demonstrated that the isocurvature
bound derived in this paper exhibits a logarithmic dependence on the infrared
cutoff, but an exponential dependence on the effective dark matter mass during
inflation. As a result, our bounds are only weakly dependent on the IR cutoff.
Moreover, we have shown that our findings for the Starobinsky and T-model of
inflation can be translated to one another by rescaling the parameter
$\lambda$, associated with the Hubble scale at the end of inflation,
$H_{\rm{end}}$. However, we note that the isocurvature bound depends on the
Hubble parameter at horizon exit, $H_{*}$, as shown by Eqs. (29) and (31).
Therefore, the isocurvature bound derived in this study can be extended to a
wider range of plateau-like inflationary potentials with different values of
$H_{*}$.

Looking towards the future, upcoming experiments such as CMB-S4
Abazajian:2019eic  and LiteBIRD Hazumi:2019lys  will not only be geared
towards the detection of B-modes in the CMB, but will also offer a more
comprehensive analysis of the CMB scalar power spectrum. This will lead to
stronger constraints on the scalar tilt, $n_{s}$, and tensor-to-scalar ratio,
$r$, paving the way for a new era of precision cosmology. Moreover, the
improved analysis is also expected to strengthen the constraints on
isocurvature modes. If no isocurvature modes are detected by these
experiments, the resulting bound will narrow down the parameter space for
gravitationally produced dark matter. Such an outcome will be of paramount
importance for future dark matter searches and will undoubtedly provide
important insights about the nature of dark matter.

###### Acknowledgements.

We would like to thank Mustafa Amin, Veronica Guidetti, Andrew Long, Michele
Redi, and Sai Chaitanya Tadepalli for helpful discussions. MG is supported by
the DGAPA-PAPIIT grant IA103123 at UNAM. MP acknowledges support by the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under
Germany’s Excellence Strategy – EXC 2121 “Quantum Universe” – 390833306. The
work of S.V. was supported in part by DOE grant DE-SC0022148. This work was
made possible by the support of the Institut Pascal at Université Paris-Saclay
during the Paris-Saclay Astroparticle Symposium 2022, with the support of the
P2IO Laboratory of Excellence (program “Investissements 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). Numerical results were partially obtained
from a custom Fortran code utilizing the thread-safe arbitrary precision
package MPFUN-For mpfun , at the Holiday cluster at Instituto de Física, UNAM.

##  Appendix A Quadratic Isocurvature

In this appendix, we show how to evaluate the power spectrum of the DM
isocurvature fluctuation at second order in perturbation theory. As mentioned
in Sec. IV, in the absence of a DM misalignment, it may be evaluated as
Liddle:1999pr ; Chung:2004nh ; Ling:2021zlj

$\displaystyle\mathcal{P}_{\mathcal{S}}(k)\;$
$\displaystyle=\;\frac{1}{\rho_{\chi}^{2}}\frac{k^{3}}{2\pi^{2}}\int\mathop{}\\!\mathrm{d}^{3}\boldsymbol{x}\
\langle\delta\rho_{\chi}(\boldsymbol{x})\delta\rho_{\chi}(0)\rangle
e^{-i\boldsymbol{k}\cdot\boldsymbol{x}}$ (32)
$\displaystyle=\;\frac{k^{3}}{(2\pi)^{5}\rho_{\chi}^{2}a^{8}}\int\mathop{}\\!\mathrm{d}^{3}\boldsymbol{p}\
P_{X}(p,|\boldsymbol{p}-\boldsymbol{k}|)\,,$ (33)

where the integrand in the second line is defined as

$\displaystyle P_{X}(p,q)\;$
$\displaystyle=\;|X^{\prime}_{p}|^{2}|X^{\prime}_{q}|^{2}+a^{4}m_{\chi}^{4}|X_{p}|^{2}|X_{q}|^{2}$
$\displaystyle+a^{2}m_{\chi}^{2}\left[(X_{p}X_{p}^{\prime*})(X_{q}X_{q}^{\prime*})+{\rm
h.c.}\right]\,.$ (34)

The integral (33) needs to be treated carefully. When $m_{\chi}\ll H_{\rm
end}$ and $\sigma/\lambda\ll 1$, $P_{X}(p,q)$ is divergent in the IR for both
arguments. However, this divergence can be naturally regularized by imposing a
lower momentum (IR) cutoff. As discussed in the main text, this cutoff scale
can be associated with the mode that left the horizon at the beginning of
inflation (our choice for $\mathcal{P}_{\mathcal{S}}$), or more
conservatively, with the wavenumber of the present-day horizon scale.
Additionally, it is important to evaluate the power spectrum
$\mathcal{P}_{\mathcal{S}}$ after the decoupling time of DM, which in this
case corresponds to the end of reheating. In scenarios where the amplitude of
isocurvature at CMB scales violates present constraints, the majority of the
contribution to the spectrum comes from modes excited during inflation. Hence,
it is sufficient to evaluate the spectrum during reheating, but well after the
end of inflation.

For finite DM masses and couplings, it is necessary to numerically evaluate
the spectrum and correctly exclude the IR ($p\rightarrow 0$) and the colinear
($|\boldsymbol{p}-\boldsymbol{k}|\rightarrow 0$) divergences from the
integral. We find it convenient to use the following integration strategy for
Eq. (33),

$\displaystyle\mathcal{P}_{\mathcal{S}}(k)\;$
$\displaystyle=\;\frac{k^{2}}{(2\pi)^{4}\rho_{\chi}^{2}a^{8}}\int_{0}^{\infty}\mathop{}\\!\mathrm{d}p\
p\int_{|k-p|}^{k+p}\mathop{}\\!\mathrm{d}q\ q\ P_{X}(p,q)$
$\displaystyle=\;\frac{2k^{2}}{(2\pi)^{4}\rho_{\chi}^{2}a^{8}}\,\mathcal{I}\big{(}P_{X}(p,q)\big{)}\,,$
(35)

with

$\mathcal{I}\big{(}f(p,q)\big{)}\,\equiv\,\left(\int_{k_{\rm
IR}}^{k/2}\mathop{}\\!\mathrm{d}p\int_{k-p}^{k+p}\mathop{}\\!\mathrm{d}q+\int_{k/2}^{k_{\rm
UV}-k}\mathop{}\\!\mathrm{d}p\int_{p}^{k+p}\mathop{}\\!\mathrm{d}q+\int_{k_{\rm
UV}-k}^{k_{\rm UV}}\mathop{}\\!\mathrm{d}p\int_{p}^{k_{\rm
UV}}\mathop{}\\!\mathrm{d}q\right)pq\,f(p,q)\,.$ (36)

We have included a UV cutoff for definiteness, which may be identified with
the decoupling scale, last scale excited at the end of reheating. However,
this scale is not necessary for the convergence of the integral in the UV. The
expression (35) is valid under the assumption that $2k_{\rm IR}<k<k_{\rm
UV}-k_{\rm IR}$. The spectrum for wavenumbers near the boundary values can be
determined by identifying the proper integration limits, as shown in Fig. 6.

$k_{\rm UV}$$k_{\rm UV}$$k_{\rm IR}$$k$$k$$q$$p$ Figure 6: Illustration of the
integration domain in Eq. (36), corresponding to the case with $2k_{\rm
IR}<k<k_{\rm UV}-k_{\rm IR}$.

##  Appendix B Analytical Derivation of the Isocurvature Power Spectrum

In this section, we provide a detailed analytical derivation of the
isocurvature power spectrum. We approximate the dynamics of the universe
during inflation as a de Sitter phase characterized by a constant accelerated
expansion with $a(t)=e^{Ht}$ and a constant Hubble parameter. We assume that
this phase begins at some initial time $a(t_{i})=a_{i}$ and ends at
$a(t_{\rm{end}})=a_{\rm{end}}$. To obtain the isocurvature power spectrum, we
solve the mode equation for $\chi_{k}=X_{k}/a$ expressed as a function of
cosmic time, which is given by

$\ddot{\chi}_{k}+3H\dot{\chi}_{k}+\left(\dfrac{k^{2}}{a^{2}}+m_{\chi}^{2}\right)\chi_{k}\,=\,0\,,$
(37)

where dots denote derivatives with respect to cosmic time and $m_{\chi}$ is
the effective dark matter mass, which is assumed to be constant. Assuming
Bunch-Davies boundary conditions and a pure de Sitter phase, the solution to
the previous equation can be expressed as

$\chi_{k}(t)\,=\,\dfrac{\sqrt{\pi}}{2a^{3/2}\sqrt{H}}e^{i\frac{\pi}{2}(\nu+\frac{1}{2})}H_{\nu}^{(1)}\left(\dfrac{k}{aH}\right)\,,$
(38)

where $H_{\nu}^{(1)}$ is the Hankel function of the first kind, and

$\nu\,\equiv\,\sqrt{\dfrac{9}{4}-\dfrac{m_{\chi}^{2}}{H^{2}}}\,.$ (39)

In the large wavelength limit $k\ll aH$, this solution reduces to

$\chi_{k}(k\ll
aH)\,\simeq\,\dfrac{1}{2\sqrt{\pi}a^{3/2}\sqrt{H}}e^{i\frac{\pi}{2}(\nu-\frac{1}{2})}\left(e^{-i\pi\nu}\Gamma(-\nu)\left(\dfrac{k}{2aH}\right)^{\nu}+\Gamma(\nu)\left(\dfrac{k}{2aH}\right)^{-\nu}\right)\,,$
(40)

in agreement with the results from Ref. Kolb:2023dzp . The power spectrum of
isocurvature perturbations can be estimated at the end of inflation by
neglecting the time derivative of mode functions in Eq. (34). Then Eq. (35)
becomes

$\displaystyle\mathcal{P}_{\mathcal{S}}(k)\;$
$\displaystyle=\;\frac{k^{2}}{(2\pi)^{4}\langle\chi^{2}\rangle^{2}}\int_{0}^{\infty}\mathop{}\\!\mathrm{d}p\
p\int_{|k-p|}^{k+p}\mathop{}\\!\mathrm{d}q\ q\
|\chi_{p}|^{2}|\chi_{q}|^{2}\,,$
$\displaystyle=\;\frac{2k^{2}}{(2\pi)^{4}\langle\chi^{2}\rangle^{2}}\,\mathcal{I}\big{(}|\chi_{p}|^{2}|\chi_{q}|^{2}\big{)}\,,$
(41)

where we approximated $\rho_{\chi}\simeq m_{\chi}^{2}\langle\chi^{2}\rangle$,
which is valid at later times. The variance of the dark matter field $\chi$ is
given by

$\langle\chi^{2}\rangle\,\equiv\,\int_{k_{\text{IR}}}^{k_{\text{UV}}}\,\mathcal{P}_{\chi}\,\dfrac{\mathop{}\\!\mathrm{d}k}{k}\,,$
(42)

with the power spectrum defined as

$\mathcal{P_{\chi}}\,\equiv\,\dfrac{k^{3}}{2\pi^{2}}|\chi_{k}|^{2}\,.$ (43)

Here, we introduce long wavelength (IR) and short wavelength (UV) cutoffs to
regularize the power spectrum, following the approach outlined in
Vilenkin:1983xp , and define $k_{\text{IR}}=a_{i}H$ and
$k_{\text{UV}}=a_{\rm{end}}H$. We use the notation
$N_{\text{tot}}\equiv\log(a_{\rm{end}}/a_{i})$ to represent the total duration
of the de Sitter era in terms of $e$-folds, while $\Delta N_{*}>0$ denotes the
duration of inflation in $e$-folds between the CMB pivot scale $k_{*}$
crossing and the end of inflation. In what follows, we calculate the power
spectra for $\chi$, variances, and the corresponding isocurvature power
spectrum evaluated at $k_{*}$. We consider three different regimes.

Negligible mass limit $m_{\chi}\ll H$. If the DM mass can be neglected
compared to the Hubble scale during the de Sitter phase, then $\nu\simeq 3/2$,
and the mode function can be expressed as

$|\chi_{k}(k\ll
aH)|\,\simeq\,\frac{1}{\sqrt{2}\sqrt{a^{3}H}}\left(\frac{k}{aH}\right)^{-3/2}\,=\frac{H}{\sqrt{2}k^{3/2}}\,,$
(44)

which leads to a well-known scale-invariant result for the power spectrum
$\mathcal{P_{\chi}}\,=\,H^{2}/(4\pi^{2})$. Combining Eq. (44) with Eq. (19)
and neglecting the time derivative of the mode function, we find that in the
long wavelength (IR, $k\ll aH$) regime, the phase space distribution
$f_{\chi}(k,t)$ scales as $f_{\chi}(k,t)\sim|\chi_{k}|^{2}\sim k^{3}$. The
variance of the field at the end of the de Sitter phase is

$\langle\chi^{2}\rangle\,\simeq\,\dfrac{H^{2}}{4\pi^{2}}\log\left(\dfrac{a_{\rm{end}}}{a_{i}}\right)\,=\,\dfrac{H^{2}}{4\pi^{2}}N_{\text{tot}}\,.$
(45)

Evaluated for the CMB crossing scale, the isocurvature power spectrum (41) is

$\mathcal{P}^{m_{\chi}\ll
H}_{\mathcal{S}}(k_{*})\,\simeq\,\dfrac{N_{\text{tot}}-\Delta
N_{*}}{N_{\text{tot}}^{2}}\,,$ (46)

which scales as $1/N_{\text{tot}}$ for $\Delta N_{*}\ll N_{\text{tot}}$.

Small mass limit $m_{\chi}\lesssim H$. If the DM mass is small but not
negligible compared to the Hubble scale, $\nu\simeq 3/2-\beta$ with
$\beta\equiv m_{\chi}^{2}/(3H^{2})$, then the mode function can be
approximated by Mukhanov:1990me ; Liddle:1999pr

$|\chi_{k}(k\ll
aH)|\,\simeq\,\dfrac{1}{\sqrt{2}}H(aH)^{-3/2}\left(\frac{k}{aH}\right)^{-\nu}\,,$
(47)

which corresponds to a blue-tilted power spectrum

$\mathcal{P_{\chi}}\,=\dfrac{H^{2}}{4\pi^{2}}\left(\dfrac{k}{aH}\right)^{2\beta}\,,$
(48)

when $\beta>0$. In this regime, the phase space distribution $f_{\chi}(k,t)$
scales as $f_{\chi}(k,t)\sim|\chi_{k}|^{2}\sim k^{-2\nu}\sim k^{2\beta-3}$.
The variance of the field at the end of the de Sitter phase is given by

$\langle\chi^{2}\rangle\,\simeq\,\left(\frac{1}{aH}\right)^{2\beta}\frac{H^{2}}{8\pi^{2}\beta}\left(k_{\text{UV}}^{2\beta}-k_{\text{IR}}^{2\beta}\right)\,.$
(49)

If $\beta>10^{-2}$ and $k_{\text{IR}}=a_{i}H\simeq 10^{-10}k_{*}$ and
$k_{\text{UV}}=a_{\rm{end}}H\simeq 10^{24}k_{*}$, the variance becomes
insensitive to the IR. In this regime, the isocurvature power spectrum can be
expressed as a sum of hypergeometric functions. Alternatively, one can isolate
the $\beta$ dependence and derive the following approximate expression

$\mathcal{P}^{m_{\chi}\lesssim
H}_{\mathcal{S}}(k_{*})\;\simeq\;\dfrac{\mathcal{P}^{m_{\chi}\ll
H}_{\mathcal{S}}(k_{*})}{\mathscr{F}_{0}}\,\mathscr{F}\left(\dfrac{m_{\chi}^{2}}{3H^{2}}\right)\,,$
(50)

where

$\mathscr{F}(x)\,\equiv\,\dfrac{x^{2}e^{-4\Delta
N_{*}x}}{16\left(1-e^{-2N_{\text{tot}}x}\right)^{2}}\,.$ (51)

and $\mathscr{F}_{0}=\lim_{x\rightarrow
0}\mathscr{F}(x)=1/(64N_{\text{tot}}^{2})$. This function shows that the
isocurvature decays exponentially with the ratio of the DM mass to the Hubble
scale, as shown in Redi:2022zkt . Note that the function $\mathscr{F}(x)$ is
insensitive to $N_{\text{tot}}$ in the limit $N_{\text{tot}}x\gg 1$.
Therefore, the only significant infrared sensitivity of Eq. (50) arises from
the scaling $\mathcal{P}^{m_{\chi}\ll
H}_{\mathcal{S}}(k_{*})/\mathscr{F}_{0}\sim N_{\text{tot}}$ at large
$N_{\text{tot}}$, which is equivalent to a logarithmic scaling in terms of the
infrared cutoff, $\sim\log(k/k_{\text{IR}})$. A comparison between the Eq.
(50) and a numerical evaluation of the isocurvature power spectrum of Eq. (41)
is shown in Fig. 7, indicating good agreement between both results.

Figure 7: Isocurvature power spectrum as a function of dark matter mass. The
red dashed line represents results based on the analytical solution for the
mode function in the negligible mass limit of Eq. (46). The orange dashed line
represents the approximated expression of Eq. (50) in the small mass regime,
and the black line a numerical evaluation of Eq. (41) in the same regime for
comparison.

Large mass limit $m_{\chi}\gtrsim H$. For $m_{\chi}^{2}>9H^{2}/4$,
$\nu\in\mathbb{C}$. We define $\nu\equiv-i\tilde{\nu}$ with
$\tilde{\nu}\in\mathbb{R}^{+}$ and $\tilde{\nu}^{2}\simeq m_{\chi}^{2}/H^{2}$.
In this case, the mode function can be expressed as

$\chi_{k}(k\ll
aH)\,\simeq\,\dfrac{1}{2\sqrt{\pi}a^{3/2}\sqrt{H}}e^{\frac{\pi}{2}(-\tilde{\nu}-\frac{1}{2}i)}\big{[}e^{\pi\tilde{\nu}}\Gamma(-i\tilde{\nu})e^{i\tilde{\nu}\log(x/2)}+\Gamma(i\tilde{\nu})e^{-i\tilde{\nu}\log(x/2)}\big{]}\,,$
(52)

where $x\equiv k/(aH)$, and

$|\chi_{k}(k\ll
aH)|\,\simeq\,\dfrac{e^{-\pi\tilde{\nu}/2}}{2\sqrt{\pi}a^{3/2}\sqrt{H}}\,\mathcal{F}(\tilde{\nu})\,,$
(53)

where the function $\mathcal{F}(\tilde{\nu})$ is defined as

$\mathcal{F}(\tilde{\nu})\,\equiv\,\left(e^{\pi\tilde{\nu}}\Big{[}2\cos\big{(}2\tilde{\nu}\log(x/2)\big{)}\Big{(}\mathcal{R}_{\tilde{\nu}}^{2}-\mathcal{I}_{\tilde{\nu}}^{2}\Big{)}+4\mathcal{R}_{\tilde{\nu}}\mathcal{I}_{\tilde{\nu}}\sin\big{(}2\tilde{\nu}\log(x/2)\big{)}\Big{]}+\left(1+e^{2\pi\tilde{\nu}}\right)\left(\dfrac{\pi}{\tilde{\nu}\sinh(\pi\tilde{\nu})}\right)\right)^{1/2}\,.$
(54)

Here $\mathcal{R}_{\tilde{\nu}}\equiv\text{Re}(\Gamma(i\tilde{\nu}))$ and
$\mathcal{I}_{\tilde{\nu}}\equiv\text{Im}(\Gamma(i\tilde{\nu}))$. The first
term in squared brackets in this expression is an oscillatory term with
frequency $2\tilde{\nu}\log(k/(2aH))$, which can typically be neglected. This
function can be well approximated in the limit of large $\tilde{\nu}$ by

$\mathcal{F}(\tilde{\nu}\gg
1)\,\simeq\,\sqrt{\dfrac{2\pi}{\tilde{\nu}}}e^{\pi\tilde{\nu}/2}\,.$ (55)

which is a good approximation even for $\tilde{\nu}\sim\mathcal{O}(1)$. The
modulus of the mode function can be approximated as

$|\chi_{k}(k\ll
aH)|\,\simeq\,\dfrac{1}{a^{3/2}\sqrt{H}}\dfrac{1}{\sqrt{2\tilde{\nu}}}\,\simeq\dfrac{1}{\sqrt{2m_{\chi}a^{3}}}\,,$
(56)

and the power spectrum can be expressed as

$\mathcal{P_{\chi}}\,\simeq\,\dfrac{1}{4\pi^{2}a^{3}m_{\chi}}k^{3}\,\simeq\dfrac{1}{4\pi^{2}m_{\chi}}H^{3}e^{3(N_{k}-N)}\,\,.$
(57)

Here $N_{k}$ is the number of $e$-folds at which the scale $k$ crosses the
horizon $a(N_{k})H=k$. The spectrum is well behaved in the IR but decreases
expentially with the number of $e$-folds after the horizon crossing
$N-N_{k}>0$, redshifting as non-relativistic matter. The variance of the field
at the end of de Sitter phase is Enqvist:1987au

$\langle\chi^{2}\rangle\,\simeq\,\dfrac{H^{3}}{12\pi^{2}m_{\chi}}\,,$ (58)

where we consider $m_{\chi}\sim H$.111111In this appendix, we use the pure de
Sitter solutions to the mode equations as an approximation to compute the
power spectra expected after a finite quasi-de Sitter phase. It should be
noted that when $\nu$ is imaginary, a naive UV cutoff in a pure de Sitter
universe may not be sufficient to subtract all the UV divergences, and higher-
order adiabatic subtraction terms may be required instead. However, here we
only want to illustrate that the isocurvature power spectrum is suppressed and
a naive UV cutoff is sufficient for our treatment in the limit $m_{\chi}\sim
H$. For a detailed discussion of the optimal order of truncation of the
adiabatic expansion, see Refs. Dabrowski:2014ica ; Corba:2022ugu  a detailed
discussion of the optimal order of truncation of the adiabatic expansion. In
this limit, the isocurvature power spectrum is exponentially suppressed:

$\mathcal{P}^{m_{\chi}\gtrsim
H}_{\mathcal{S}}(k_{*})\;\simeq\;\frac{3e^{-3\Delta N_{*}}}{2}\,,$ (59)

and for $\Delta N_{*}=55$ we find $\mathcal{P}^{m_{\chi}\gtrsim
H}_{\mathcal{S}}(k_{*})\simeq\mathcal{O}(10^{-72})$.

Estimating the Isocurvature Constraint. From Eq. (59), one can readily see
that the isocurvature constraint
$\mathcal{P}_{\mathcal{S}}(k_{*})<\beta_{\text{iso}}\mathcal{P}_{\mathcal{R}}(k_{*})\simeq
10^{-11}$ is always satisfied for $m_{\chi}^{2}>9H^{2}/4$.

In the limit of negligible mass, where $m_{\chi}\ll H$, as described in Eq.
(46), the isocurvature power spectrum is too large by orders of magnitude,
given by $\mathcal{P}^{m_{\chi}\ll
H}_{\mathcal{S}}(k_{*})\,\sim\,\mathcal{O}(10^{-3})$ for
$N_{\text{tot}}\,\sim\,\mathcal{O}(100)$.

In the intermediate regime, where $m_{\chi}\lesssim H$, we can use the
approximate expression (50) to compute the limit

$m_{\chi}\,\gtrsim\,0.57\,H\,,$ (60)

where we used $\Delta N_{*}=55$ and $N_{\rm{tot}}=76.5$ $e$-folds. This result
agrees extremely well with the fully numerical solution Eq. (29). However, for
the Starobinsky and T-models of inflation, considered in this work, the Hubble
rate varies with time. One can therefore evaluate Eq. (60) in terms of the
Hubble scale at the end of inflation,

$\displaystyle m_{\chi}\gtrsim
1.18H_{\rm{end}}\qquad(\text{Starobinsky})\,,\leavevmode\nobreak\ $
$\displaystyle m_{\chi}\gtrsim 1.41H_{\rm{end}}\qquad(\text{T-model})\,.$ (61)

These results agree extremely well with the fully numerical result (IV).
Assuming that the bare DM mass is negligible compared to the contribution
induced by the inflaton $m_{\rm{eff}}^{2}\simeq\sigma\phi^{2}$, the bound (60)
corresponds to

$\displaystyle\dfrac{\sigma}{\lambda}\,>\,0.006\,(0.02)\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ {\textrm{for
Starobinsky\leavevmode\nobreak\ (T-model)}}\,,\leavevmode\nobreak\ $

which is in good agreement with the fully numerical result (IV).

## References

  * (1) XENON Collaboration, E. Aprile et al., “Dark Matter Search Results from a One Ton-Year Exposure of XENON1T,” Phys. Rev. Lett. 121 no. 11, (2018) 111302, arXiv:1805.12562 [astro-ph.CO].
  * (2) LUX Collaboration, D. S. Akerib et al., “Results from a search for dark matter in the complete LUX exposure,” Phys. Rev. Lett. 118 no. 2, (2017) 021303, arXiv:1608.07648 [astro-ph.CO].
  * (3) PandaX-4T Collaboration, Y. Meng et al., “Dark Matter Search Results from the PandaX-4T Commissioning Run,” Phys. Rev. Lett. 127 no. 26, (2021) 261802, arXiv:2107.13438 [hep-ex].
  * (4) LZ Collaboration, J. Aalbers et al., “First Dark Matter Search Results from the LUX-ZEPLIN (LZ) Experiment,” arXiv:2207.03764 [hep-ex].
  * (5) G. Arcadi, M. Dutra, P. Ghosh, M. Lindner, Y. Mambrini, M. Pierre, S. Profumo, and F. S. Queiroz, “The waning of the WIMP? A review of models, searches, and constraints,” Eur. Phys. J. C 78 no. 3, (2018) 203, arXiv:1703.07364 [hep-ph].
  * (6) M. Escudero, A. Berlin, D. Hooper, and M.-X. Lin, “Toward (Finally!) Ruling Out Z and Higgs Mediated Dark Matter Models,” JCAP 12 (2016) 029, arXiv:1609.09079 [hep-ph].
  * (7) Y. Mambrini, Particles in the dark Universe.  No. ISBN 978-3-030-78139-2. Springer, 2021.
  * (8) L. J. Hall, K. Jedamzik, J. March-Russell, and S. M. West, “Freeze-In Production of FIMP Dark Matter,” JHEP 03 (2010) 080, arXiv:0911.1120 [hep-ph].
  * (9) N. Bernal, M. Heikinheimo, T. Tenkanen, K. Tuominen, and V. Vaskonen, “The Dawn of FIMP Dark Matter: A Review of Models and Constraints,” Int. J. Mod. Phys. A 32 no. 27, (2017) 1730023, arXiv:1706.07442 [hep-ph].
  * (10) L. Parker, “Particle creation in expanding universes,” Phys. Rev. Lett. 21 (1968) 562–564.
  * (11) L. Parker, “Quantized fields and particle creation in expanding universes. 1.,” Phys. Rev. 183 (1969) 1057–1068.
  * (12) L. Parker, “Quantized fields and particle creation in expanding universes. 2.,” Phys. Rev. D 3 (1971) 346–356. [Erratum: Phys.Rev.D 3, 2546–2546 (1971)].
  * (13) L. H. Ford, “Gravitational Particle Creation and Inflation,” Phys. Rev. D 35 (1987) 2955.
  * (14) Y. Ema, R. Jinno, K. Mukaida, and K. Nakayama, “Gravitational Effects on Inflaton Decay,” JCAP 05 (2015) 038, arXiv:1502.02475 [hep-ph].
  * (15) Y. Ema, R. Jinno, K. Mukaida, and K. Nakayama, “Gravitational particle production in oscillating backgrounds and its cosmological implications,” Phys. Rev. D 94 no. 6, (2016) 063517, arXiv:1604.08898 [hep-ph].
  * (16) Y. Ema, K. Nakayama, and Y. Tang, “Production of Purely Gravitational Dark Matter,” JHEP 09 (2018) 135, arXiv:1804.07471 [hep-ph].
  * (17) M. Herranen, T. Markkanen, S. Nurmi, and A. Rajantie, “Spacetime curvature and Higgs stability after inflation,” Phys. Rev. Lett. 115 (2015) 241301, arXiv:1506.04065 [hep-ph].
  * (18) T. Markkanen and S. Nurmi, “Dark matter from gravitational particle production at reheating,” JCAP 02 (2017) 008, arXiv:1512.07288 [astro-ph.CO].
  * (19) M. A. G. Garcia, M. Pierre, and S. Verner, “Scalar dark matter production from preheating and structure formation constraints,” Phys. Rev. D 107 no. 4, (2023) 043530, arXiv:2206.08940 [hep-ph].
  * (20) D. J. H. Chung, E. W. Kolb, and A. Riotto, “Superheavy dark matter,” Phys. Rev. D 59 (1998) 023501, arXiv:hep-ph/9802238.
  * (21) D. J. H. Chung, E. W. Kolb, and A. Riotto, “Nonthermal supermassive dark matter,” Phys. Rev. Lett. 81 (1998) 4048–4051, arXiv:hep-ph/9805473.
  * (22) V. Kuzmin and I. Tkachev, “Matter creation via vacuum fluctuations in the early universe and observed ultrahigh-energy cosmic ray events,” Phys. Rev. D 59 (1999) 123006, arXiv:hep-ph/9809547.
  * (23) S. Ling and A. J. Long, “Superheavy scalar dark matter from gravitational particle production in $\alpha$-attractor models of inflation,” Phys. Rev. D 103 no. 10, (2021) 103532, arXiv:2101.11621 [astro-ph.CO].
  * (24) K. Kaneta, S. M. Lee, and K.-y. Oda, “Boltzmann or Bogoliubov? Approaches compared in gravitational particle production,” JCAP 09 (2022) 018, arXiv:2206.10929 [astro-ph.CO].
  * (25) D. J. H. Chung, E. W. Kolb, A. Riotto, and L. Senatore, “Isocurvature constraints on gravitationally produced superheavy dark matter,” Phys. Rev. D 72 (2005) 023511, arXiv:astro-ph/0411468.
  * (26) D. J. H. Chung and H. Yoo, “Isocurvature Perturbations and Non-Gaussianity of Gravitationally Produced Nonthermal Dark Matter,” Phys. Rev. D 87 (2013) 023516, arXiv:1110.5931 [astro-ph.CO].
  * (27) D. J. H. Chung and H. Yoo, “Elementary Theorems Regarding Blue Isocurvature Perturbations,” Phys. Rev. D 91 (2015) 083530, arXiv:1501.05618 [astro-ph.CO].
  * (28) N. Herring, D. Boyanovsky, and A. R. Zentner, “Nonadiabatic cosmological production of ultralight dark matter,” Phys. Rev. D 101 no. 8, (2020) 083516, arXiv:1912.10859 [gr-qc].
  * (29) L. E. Padilla, J. A. Vázquez, T. Matos, and G. Germán, “Scalar Field Dark Matter Spectator During Inflation: The Effect of Self-interaction,” JCAP 05 (2019) 056, arXiv:1901.00947 [astro-ph.CO].
  * (30) M. Redi and A. Tesi, “Dark photon Dark Matter without Stueckelberg mass,” JHEP 10 (2022) 167, arXiv:2204.14274 [hep-ph].
  * (31) Planck Collaboration, Y. Akrami et al., “Planck 2018 results. X. Constraints on inflation,” Astron. Astrophys. 641 (2020) A10, arXiv:1807.06211 [astro-ph.CO].
  * (32) C. Gordon, D. Wands, B. A. Bassett, and R. Maartens, “Adiabatic and entropy perturbations from inflation,” Phys. Rev. D 63 (2000) 023506, arXiv:astro-ph/0009131.
  * (33) S. Renaux-Petel and K. Turzynski, “On reaching the adiabatic limit in multi-field inflation,” JCAP 06 (2015) 010, arXiv:1405.6195 [astro-ph.CO].
  * (34) A. A. Starobinsky, “A New Type of Isotropic Cosmological Models Without Singularity,” Phys. Lett. B 91 (1980) 99–102.
  * (35) J. Ellis, D. V. Nanopoulos, and K. A. Olive, “Starobinsky-like Inflationary Models as Avatars of No-Scale Supergravity,” JCAP 10 (2013) 009, arXiv:1307.3537 [hep-th].
  * (36) R. Kallosh and A. Linde, “Non-minimal Inflationary Attractors,” JCAP 10 (2013) 033, arXiv:1307.7938 [hep-th].
  * (37) M. A. G. Garcia, K. Kaneta, Y. Mambrini, and K. A. Olive, “Inflaton Oscillations and Post-Inflationary Reheating,” JCAP 04 (2021) 012, arXiv:2012.10756 [hep-ph].
  * (38) J. Ellis, M. A. G. Garcia, D. V. Nanopoulos, K. A. Olive, and S. Verner, “BICEP/Keck Constraints on Attractor Models of Inflation and Reheating,” arXiv:2112.04466 [hep-ph].
  * (39) T. Opferkuch, P. Schwaller, and B. A. Stefanek, “Ricci Reheating,” JCAP 07 (2019) 016, arXiv:1905.06823 [gr-qc].
  * (40) Y. Mambrini and K. A. Olive, “Gravitational Production of Dark Matter during Reheating,” Phys. Rev. D 103 no. 11, (2021) 115009, arXiv:2102.06214 [hep-ph].
  * (41) M. R. Haque and D. Maity, “Gravitational dark matter: free streaming and phase space distribution,” arXiv:2112.14668 [hep-ph].
  * (42) S. Clery, Y. Mambrini, K. A. Olive, and S. Verner, “Gravitational portals in the early Universe,” Phys. Rev. D 105 no. 7, (2022) 075005, arXiv:2112.15214 [hep-ph].
  * (43) S. Clery, Y. Mambrini, K. A. Olive, A. Shkerin, and S. Verner, “Gravitational portals with nonminimal couplings,” Phys. Rev. D 105 no. 9, (2022) 095042, arXiv:2203.02004 [hep-ph].
  * (44) B. Barman and N. Bernal, “Gravitational SIMPs,” JCAP 06 (2021) 011, arXiv:2104.10699 [hep-ph].
  * (45) B. Barman, N. Bernal, Y. Xu, and O. Zapata, “Ultraviolet Freeze-in with a Time-dependent Inflaton Decay,” arXiv:2202.12906 [hep-ph].
  * (46) M. R. Haque and D. Maity, “Gravitational Reheating,” arXiv:2201.02348 [hep-ph].
  * (47) M. A. G. Garcia, K. Kaneta, Y. Mambrini, K. A. Olive, and S. Verner, “Freeze-in from preheating,” JCAP 03 no. 03, (2022) 016, arXiv:2109.13280 [hep-ph].
  * (48) L. Kofman, A. D. Linde, and A. A. Starobinsky, “Towards the theory of reheating after inflation,” Phys. Rev. D 56 (1997) 3258–3295, arXiv:hep-ph/9704452.
  * (49) J. Ellis, M. A. G. García, D. V. Nanopoulos, and K. A. Olive, “Two-Field Analysis of No-Scale Supergravity Inflation,” JCAP 01 (2015) 010, arXiv:1409.8197 [hep-ph].
  * (50) L. H. Ford, “Cosmological particle production: a review,” Rept. Prog. Phys. 84 no. 11, (2021) , arXiv:2112.02444 [gr-qc].
  * (51) A. A. Starobinsky and J. Yokoyama, “Equilibrium state of a selfinteracting scalar field in the De Sitter background,” Phys. Rev. D 50 (1994) 6357–6368, arXiv:astro-ph/9407016.
  * (52) G. N. Felder, L. Kofman, and A. D. Linde, “Gravitational particle production and the moduli problem,” JHEP 02 (2000) 027, arXiv:hep-ph/9909508.
  * (53) T. Tenkanen, “Dark matter from scalar field fluctuations,” Phys. Rev. Lett. 123 no. 6, (2019) 061302, arXiv:1905.01214 [astro-ph.CO].
  * (54) A. R. Liddle and A. Mazumdar, “Perturbation amplitude in isocurvature inflation scenarios,” Phys. Rev. D 61 (2000) 123507, arXiv:astro-ph/9912349.
  * (55) Planck Collaboration, N. Aghanim et al., “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641 (2020) A6, arXiv:1807.06209 [astro-ph.CO]. [Erratum: Astron.Astrophys. 652, C4 (2021)].
  * (56) V. K. Narayanan, D. N. Spergel, R. Dave, and C.-P. Ma, “Constraints on the mass of warm dark matter particles and the shape of the linear power spectrum from the Ly$\alpha$ forest,” Astrophys. J. Lett. 543 (2000) L103–L106, arXiv:astro-ph/0005095.
  * (57) M. Viel, J. Lesgourgues, M. G. Haehnelt, S. Matarrese, and A. Riotto, “Constraining warm dark matter candidates including sterile neutrinos and light gravitinos with WMAP and the Lyman-alpha forest,” Phys. Rev. D 71 (2005) 063534, arXiv:astro-ph/0501562.
  * (58) M. Viel, G. D. Becker, J. S. Bolton, and M. G. Haehnelt, “Warm dark matter as a solution to the small scale crisis: New constraints from high redshift Lyman-$\alpha$ forest data,” Phys. Rev. D 88 (2013) 043502, arXiv:1306.2314 [astro-ph.CO].
  * (59) J. Baur, N. Palanque-Delabrouille, C. Yèche, C. Magneville, and M. Viel, “Lyman-alpha Forests cool Warm Dark Matter,” JCAP 08 (2016) 012, arXiv:1512.01981 [astro-ph.CO].
  * (60) V. Irsic et al., “New Constraints on the free-streaming of warm dark matter from intermediate and small scale Lyman-$\alpha$ forest data,” Phys. Rev. D 96 no. 2, (2017) 023522, arXiv:1702.01764 [astro-ph.CO].
  * (61) N. Palanque-Delabrouille, C. Yèche, N. Schöneberg, J. Lesgourgues, M. Walther, S. Chabanier, and E. Armengaud, “Hints, neutrino bounds and WDM constraints from SDSS DR14 Lyman-$\alpha$ and Planck full-survey data,” JCAP 04 (2020) 038, arXiv:1911.09073 [astro-ph.CO].
  * (62) A. Garzilli, O. Ruchayskiy, A. Magalich, and A. Boyarsky, “How warm is too warm? Towards robust Lyman-$\alpha$ forest bounds on warm dark matter,” arXiv:1912.09397 [astro-ph.CO].
  * (63) G. Ballesteros, M. A. G. Garcia, and M. Pierre, “How warm are non-thermal relics? Lyman-$\alpha$ bounds on out-of-equilibrium dark matter,” JCAP 03 (2021) 101, arXiv:2011.13458 [hep-ph].
  * (64) K. Abazajian et al., “CMB-S4 Science Case, Reference Design, and Project Plan,” arXiv:1907.04473 [astro-ph.IM].
  * (65) M. Hazumi et al., “LiteBIRD: A Satellite for the Studies of B-Mode Polarization and Inflation from Cosmic Background Radiation Detection,” J. Low Temp. Phys. 194 no. 5-6, (2019) 443–452.
  * (66) D. H. Bailey, “MPFUN2015: A thread-safe arbitrary precision package,”. http://www.davidhbailey.com/dhbpapers/mpfun2015.pdf.
  * (67) E. W. Kolb, S. Ling, A. J. Long, and R. A. Rosen, “Cosmological gravitational particle production of massive spin-2 particles,” arXiv:2302.04390 [astro-ph.CO].
  * (68) A. Vilenkin, “Quantum Fluctuations in the New Inflationary Universe,” Nucl. Phys. B 226 (1983) 527–546.
  * (69) V. F. Mukhanov, H. A. Feldman, and R. H. Brandenberger, “Theory of cosmological perturbations. Part 1. Classical perturbations. Part 2. Quantum theory of perturbations. Part 3. Extensions,” Phys. Rept. 215 (1992) 203–333.
  * (70) K. Enqvist, K. W. Ng, and K. A. Olive, “Scalar Field Fluctuations in the Early Universe,” Nucl. Phys. B 303 (1988) 713.
  * (71) R. Dabrowski and G. V. Dunne, “Superadiabatic particle number in Schwinger and de Sitter particle production,” Phys. Rev. D 90 no. 2, (2014) 025021, arXiv:1405.0302 [hep-th].
  * (72) S. P. Corbà and L. Sorbo, “On the adiabatic subtraction of cosmological perturbations,” arXiv:2209.14362 [hep-th].