# Inflaton Production of Scalar Dark Matter through Fluctuations and
Scattering

Gongjun Choi  William I. Fine Theoretical Physics Institute, School of Physics
and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA  Marcos A.
G. Garcia  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  Wenqi Ke
William I. Fine Theoretical Physics Institute, School of Physics and
Astronomy, University of Minnesota, Minneapolis, MN 55455, USA  Yann Mambrini
Sorbonne Université, CNRS, Laboratoire de Physique Théorique et Hautes
Energies, LPTHE, F-75005 Paris, France  Keith A. Olive  William I. Fine
Theoretical Physics Institute, School of Physics and Astronomy, University of
Minnesota, Minneapolis, MN 55455, USA  Sarunas Verner  Institute for
Fundamental Theory, Physics Department, University of Florida, Gainesville, FL
32611, USA

###### Abstract

We study the effects on particle production of a Planck-suppressed coupling
between the inflaton and a scalar dark matter candidate, $\chi$. In the
absence of this coupling the dominant source for the relic density of $\chi$
is the long wavelength modes produced from the scalar field fluctuations
during inflation. In this case, there are strong constraints on the mass of
the scalar and the reheating temperature after inflation from the present-day
relic density of $\chi$ (assuming $\chi$ is stable). When a coupling
$\sigma\phi^{2}\chi^{2}$ is introduced, with
$\sigma={\tilde{\sigma}}m_{\phi}^{2}/M_{P}^{2}\sim 10^{-10}{\tilde{\sigma}}$,
where $m_{\phi}$ is the inflaton mass, the allowed parameter space begins to
open up considerably even for ${\tilde{\sigma}}$ as small as $\gtrsim
10^{-7}$. For ${\tilde{\sigma}}\gtrsim\frac{9}{16}$, particle production is
dominated by the scattering of the inflaton condensate, either through single
graviton exchange or the contact interaction between $\phi$ and $\chi$. In
this regime, the range of allowed masses and reheating temperatures is
maximal. For $0.004<{\tilde{\sigma}}<50$, constraints from isocurvature
fluctuations are satisfied, and the production from parametric resonance can
be neglected.

††preprint: UMN–TH–4321/24††preprint: FTPI–MINN–24/13

##  I Introduction

Quantum scalar fluctuations are an inevitable consequence of a de Sitter
phase, including inflation fluc . During inflation, these fluctuations grow
linearly in time, reaching an asymptotic value $\langle\chi^{2}\rangle\sim
H_{I}^{4}/m_{\chi}^{2}$ for a massive scalar field, $\chi$, with a mass (much)
smaller than the Hubble scale during inflation, $m_{\chi}\ll H_{I}$. For
sufficiently light scalar fields, these fluctuations are very large, and the
long wavelength contributions of these fluctuations obey the classical
equations of motion longw , making them indistinguishable from a homogeneous
scalar field background.

Although there are far too many cosmological consequences of these
fluctuations to delineate here, we note their particular importance along
(SUSY) flat directions gkm  and their role in Affleck-Dine Baryogenesis AD ;
Linde:1985gh ; Campbell:1986qg ; Enqvist:2003gh ; Garcia:2013bha , in avoiding
the washout of density perturbations Graziani:1988bp ; Enqvist:2003gh ;
Enqvist:2011pt , and on reheating and thermalization Allahverdi:2005mz ;
Olive:2006uw . All of these effects rely on the fact that the flat directions
(or very light scalar fields) remain so during inflation. In general
supergravity models, this is difficult Dine:1995uk  as Hubble scale masses are
generally induced for all scalar fields. This is not the case for theories
with a Heisenberg symmetry Binetruy:1987xj  which includes no-scale
supergravity noscale , where flat directions are maintained during inflation
GMO .

Of particular interest to us here are the consequences of the fluctuations of
relatively light stable scalar fields, which can eventually play the role of
cold dark matter Turner:1987vd ; Peebles:1999fz ; Enqvist:2014zqa ;
Nurmi:2015ema ; Bertolami:2016ywc ; Alonso-Alvarez:2018tus ; Markkanen:2018gcw
; Tenkanen:2019aij ; Choi:2019mva ; Cosme:2020nac ; Lebedev:2022cic . In
general, the relic abundance of the scalar field dark matter will depend on
the scale of inflation, $H_{I}$, the equation of state during inflaton
oscillations, $w$, the reheating temperature, $T_{\rm RH}$, the mass of the
scalar field, $m_{\chi}$, and its self-interactions, $\lambda_{p}\chi^{p}$.
Qualitatively, scalar field fluctuations drive $\langle\chi^{2}\rangle$ to
large values, which source the field oscillations of $\chi$ after inflation.
Depending on the parameters, $\chi$-oscillations may begin before or after
reheating, and their subsequent evolution will depend on their interactions.
Once these parameters are specified, it is relatively straightforward to
calculate the relic density of $\chi$. We note that this production is in
addition to (and may dominate over) the production of scalar dark matter
directly from the inflaton condensate after inflation ema ; Dimopoulos:2006ms
; Graham:2015rva ; Garny:2015sjg ; Tang:2017hvq ; Bernal:2018qlk ; Ema:2019yrd
; Chianese:2020yjo ; Ahmed:2020fhc ; Kolb:2020fwh ; Redi:2020ffc ;
Ling:2021zlj ; MO ; Bernal:2021kaj ; Barman:2021ugy ; Haque:2021mab ; cmov ;
Clery:2022wib ; Mambrini:2022uol ; Basso:2022tpd ; Haque:2023yra ;
Garcia:2022vwm ; Kaneta:2022gug ; Kolb:2023dzp ; Kaneta:2023kfv ;
Garcia:2023qab ; kkmov .

Most studies assume that the scalar, $\chi$, is secluded and interacts only
through self-interactions and (minimally) through gravity. As a result, the
field is often called a spectator Markkanen:2018gcw . The presence or absence
of self-interactions plays a significant role in determining the relic
density, as the initial amplitude of $\chi$ oscillations is set by the
effective mass of $\chi$. Indeed, in the absence of self-interactions,
$\langle\chi^{2}\rangle\sim H_{I}^{4}/m_{\chi}^{2}\gg H_{I}^{2}$ for
$m_{\chi}\ll H_{I}$. However, the presence of a self-interaction of the form
$\lambda_{4}\chi^{4}$ provides a contribution to the effective mass, and we
expect $\langle\chi^{2}\rangle\sim H_{I}^{2}/\sqrt{\lambda_{4}}\sim H_{I}^{2}$
for $\lambda_{4}\sim 1$ Starobinsky:1994bd , significantly reducing the relic
abundance of $\chi$.

On the other hand, the minimal coupling to gravity ensures that the spectator
cannot be entirely absent from interactions with other fields. The minimal
coupling through single graviton exchange provides a coupling between the
inflaton which (on shell) takes the form $\sigma\phi^{2}\chi^{2}$ with
$\sigma=(m_{\phi}/2M_{P})^{2}$ cmov ; Clery:2022wib  for the case of an
inflaton potential $V(\phi)=\frac{1}{2}m_{\phi}^{2}\phi^{2}$, where
$M_{P}=1/\sqrt{8\pi G}\simeq 2.4\times 10^{18}\,\rm{GeV}$ is the reduced
Planck mass. This term, however, does not contribute to the effective mass of
$\chi$, as it originates from a term in the Lagrangian of the form,
$\phi^{2}(\partial\chi)^{2}$.111There is, in fact, a term $\phi^{2}\chi^{2}$
in the Lagrangian, but it is highly suppressed with a coupling of order
$m_{\chi}^{2}/M_{P}^{2}$ and will not appreciably affect the generation of the
fluctuation in $\chi^{2}$.  This coupling will also be present if $\chi$ is
coupled to curvature ($\mathcal{L}\ni\xi_{\chi}\chi^{2}R$, where $R$ is the
Ricci curvature). In this case, a coupling
$\xi_{\chi}(m_{\phi}^{2}/M_{P}^{2})\phi^{2}\chi^{2}$ appears when the theory
is brought to the Einstein frame Clery:2022wib .

Nevertheless, as both the inflaton and spectator are gauge singlets, there is
no reason not to expect Planck-suppressed operators of the form
$\chi^{2}\phi^{2}$, and these may greatly affect the generation and subsequent
evolution of $\langle\chi^{2}\rangle$. In this work, we study the effect of
such operators on the relic density of spectator dark matter. We consider the
simplest case where inflaton oscillations are driven by a quadratic term about
minimum of the inflaton potential ($w=0$), motivated by a simple inflationary
model Staro . Once the inflaton potential parameters are fixed from CMB
observations, we have three independent free parameters: the spectator mass,
$m_{\chi}$, the spectator coupling to the inflaton, $\sigma$, and the
reheating temperature, $T_{\rm RH}$. The relative values of these parameters
will determine when the effective spectator mass is dominated by its bare
mass, when the spectator oscillations begin, and when reheating is achieved.
The corresponding values of the cosmological scale factor for these events are
denoted respectively by $a_{\chi}$, $a_{\rm osc}$, and $a_{\rm RH}$. Without
either self-interactions or couplings to other fields, as already shown in
Cosme:2020nac ; Lebedev:2022cic , there can be a severe overproduction of
scalar dark matter, as we will also see below. Adding a coupling to the
inflaton greatly relaxes this constraint and opens up the allowed values of
$m_{\chi}$ and $T_{\rm RH}$.

We will show that the parameter space begins to open for Planck-suppressed
couplings222In the theory with a non-minimal coupling to curvature, we can
make the associations ${\tilde{\sigma}}=-\xi_{\chi}$ in the effective mass
Clery:2022wib . This association is only valid when $\xi_{\chi}\chi^{2}\ll
M_{P}^{2}$. $\sigma={\tilde{\sigma}}m_{\phi}^{2}/M_{P}^{2}$, with
${\tilde{\sigma}}\gtrsim 10^{-7}$. Furthermore, for
${\tilde{\sigma}}\gtrsim\frac{9}{16}$, particle production becomes dominated
by inflaton scattering through either single graviton exchange or the contact
interaction between $\chi$ and $\phi$. In this case the allowed range on the
mass of $\chi$ is very large, $\sim 0.01-\sim 10^{13}$ GeV and the allowed
range on $T_{\rm RH}$ also opens, reaching roughly $0.004-10^{15}$ GeV. This
can be accomplished with couplings ${\tilde{\sigma}}$ between 0.004 and 50, so
that constraints from isocurvature perturbations are avoided and non-
perturbative production from scattering can be neglected.

In what follows, we first define the model of inflation, the coupling to the
scalar, and the cosmological evolution of the system in Section II. In Section
III, we define our model parameters and consider the production of long
wavelength modes (important for relatively small couplings), combined with the
gravitational production of $\chi$ through scattering of inflaton field. We
also discuss the constraints from isocurvature perturbations. Our results are
collected in Section IV, where we distinguish between large and small coupling
scenarios. A discussion and summary of our results is given in Section V.

##  II Model and Field Evolution after Inflation

We consider a model with two real scalar fields: the inflaton, $\phi$, and the
spectator, $\chi$. The action is given by

$\displaystyle\mathcal{S}\;=\;\int
d^{4}x\sqrt{-g}\bigg{[}\frac{1}{2}\partial^{\mu}\phi\partial_{\mu}\phi$
$\displaystyle-V(\phi)+\frac{1}{2}\partial_{\mu}\chi\partial^{\mu}\chi$ (1)
$\displaystyle-\frac{1}{2}m_{\chi}^{2}\chi^{2}-\frac{1}{2}\sigma\phi^{2}\chi^{2}\bigg{]}\,,$

where $V(\phi)$ is the inflationary potential, $m_{\chi}$ is the bare mass of
the spectator field, and $\sigma$ is a direct coupling between the inflaton
and spectator field. Although this coupling appears renormalizable, we assume
it is an effective and Planck-suppressed coupling. Therefore, we expect
$\sigma$ to be of order $m_{\phi}^{2}/M_{P}^{2}$ in analogy with the on-shell
coupling from single graviton exchange cmov ; Clery:2022wib . We use the
metric signature $(+,-,-,-)$.

Varying this action with respect to $\phi$ and $\chi$, the equations of motion
are333We neglected the $\sigma\langle\chi^{2}\rangle\phi$ contribution because
it is of order $\sigma\langle\chi^{2}\rangle\ll m_{\phi}^{2}$.

$\ddot{\phi}+3H\dot{\phi}+\frac{dV(\phi)}{d\phi}\;\simeq\;0\,,$ (2)
$\ddot{\chi}+3H\dot{\chi}+m_{\chi,\rm eff}^{2}\chi\;=\;0\,,$ (3)

where $H=\frac{\dot{a}}{a}$ is the Hubble parameter, and the effective mass of
the spectator field $\chi$ is given by $m_{\chi,\rm
eff}^{2}\;=\;m_{\chi}^{2}+\sigma\langle\phi^{2}\rangle$. When accelerated
expansion ends, the inflaton begins a series of oscillations, dominating the
energy density of the Universe until reheating, when the radiation density
(produced, for example, by inflaton decays or scatterings gkmo1 ; gkmo2 )
equals the energy density in the inflaton condensate. Similarly, the spectator
field, $\chi$, will also begin its series of oscillations when $mt\simeq 1$ or
$m_{\chi,\rm eff}(a)=\frac{3}{2}H(a)$ in a matter-dominated (inflaton-
dominated) universe, occurring at $a=a_{\rm osc}$. The redshifting of the
energy density in $\chi$ is highly dependent on the content of the Universe,
and we show how the relative values of $a_{\rm osc}$ and $a_{\rm RH}$ are
crucial in determining the relic density of $\chi$ at the present time.

For numerical analysis, we consider the Starobinsky model of inflation Staro ,

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

which, when expanded about the minimum, behaves as

$V(\phi)\;\simeq\;\frac{1}{2}m_{\phi}^{2}\phi^{2}\,,\qquad\phi\lesssim
M_{P}\,,$ (5)

with $m_{\phi}=2H_{I}$, where $H_{I}$ is the Hubble parameter during
inflation. The normalization of the inflaton potential can be obtained from
the overall amplitude of scalar perturbations, $A_{S^{*}}=2.1\times 10^{-9}$
Planck:2018jri  and

$\frac{m_{\phi}^{2}}{M_{P}^{2}}\;\simeq\;\frac{24\pi^{2}A_{S^{*}}}{N_{*}^{2}}\,.$
(6)

For a nominal choice of $55$ $e$-folds,444We take $55$ $e$-folds for
simplicity. See egnov  for a detailed analysis of the Starobinsky model that
constrains the number of $e$-folds while accounting for reheating. for the
Starobinsky potential, one finds $m_{\phi}\simeq 1.25\times 10^{-5}M_{P}\simeq
3\times 10^{13}~{}{\rm GeV}$ eno6 ; EGNO5 ; egnov ; Ema:2024sit .
Additionally, the scalar tilt, $n_{s}=0.965$, and scalar-to-tensor anisotropy
ratio, $r=0.0035$, are in excellent agreement with Planck results
Planck:2018jri .

The solution to Eq. (2) is well known and leads to harmonic oscillations of
$\phi$. The energy density can be expressed in terms of the amplitude of the
oscillations and $\rho_{\phi}=\frac{1}{2}\dot{\phi}^{2}+V(\phi)=\rho_{\rm
end}(a_{\rm end}/a)^{3}$. Here, the scale factor at the end of inflation is
denoted as $a_{\rm end}$, which is derived from the fact that at $a_{\rm
end}$, $V(\phi_{\rm end})=\dot{\phi}_{\rm end}^{2}$ and $\rho_{\rm
end}=\frac{3}{2}V(\phi_{\rm end})$.555Equivalently, inflation ends when the
accelerated expansion stops and $\ddot{a}(t_{\rm end})=0$. The inflaton field
value at the end of inflation is denoted as $\phi_{\rm end}=\phi(a_{\rm
end})$, and similarly, $H_{\rm end}=H(a_{\rm end})$.

Supposing that $\rho_{\phi}$ dominates the energy budget of the Universe until
reheating is complete, the Hubble parameter scales as
$H(a)\propto\sqrt{\rho_{\phi}}\sim a^{-\frac{3}{2}}$ for $a_{\rm end}<a<a_{\rm
RH}$. It is then convenient to rewrite the equation of motion Eq. (3) for
$\chi$ as

$\chi^{\prime\prime}+\frac{5}{2}\frac{\chi^{\prime}}{a}+4\tilde{\sigma}\frac{\chi}{a^{2}}\;=\;0\,,\qquad\tilde{\sigma}=\sigma\frac{M_{P}^{2}}{m_{\phi}^{2}}\,,$
(7)

where ${}^{\prime}\equiv d/da=(aH)^{-1}d/dt$, and we assumed that the
contribution from the bare mass term $m_{\chi}^{2}$ is negligible. We also
used the relation, $4H_{\rm end}^{2}M_{P}^{2}=m_{\phi}^{2}\phi_{\rm end}^{2}$
at the end of inflation.

If we choose a general ansatz of the form $\chi(a)\propto a^{-A+iB}$, the
solution for real $B$ is given by:

$\displaystyle\chi(a)=\chi_{\rm end}\left(\frac{a_{\rm
end}}{a}\right)^{\frac{3}{4}}\left[\cos\left(B\log\left[\frac{a}{a_{\rm
end}}\right]\right)\right.$ (8)
$\displaystyle+\left.\frac{3}{4B}\sin\left(B\log\left[\frac{a}{a_{\rm
end}}\right]\right)\right]\,,$

with $\chi_{\rm end}=\chi(a_{\rm end})$, $A=\frac{3}{4}$, and
$B=2\sqrt{{\tilde{\sigma}}-\frac{9}{64}}$. The constants of integration have
been set by assuming $\chi({a_{\rm end}})=\chi_{\rm end}$ and
$\chi^{\prime}(a_{\rm end})=0$. Real $B$ then requires ${\tilde{\sigma}}\geq
9/64$. For smaller ${\tilde{\sigma}}$, the solution is no longer oscillatory
and takes the form

$\chi(a)=\chi_{\rm end}\frac{1}{2c}\left(\frac{a_{\rm
end}}{a}\right)^{\frac{c+3}{4}}\left(c-3+\left(\frac{a}{a_{\rm
end}}\right)^{\frac{c}{2}}(c+3)\right)\,,$ (9)

with $c=8\sqrt{\frac{9}{64}-{\tilde{\sigma}}}$. These solutions break down
when the bare mass term for $\chi$ cannot be neglected in Eq. (7). This occurs
when ${\tilde{\sigma}}$ is very small, or at large $a$, when $\phi^{2}$ is
small. Then, Eq. (7) is modified (by replacing the third term with
$4{\tilde{\sigma}}/a^{2}\rightarrow m_{\chi}^{2}a/H_{\rm end}^{2}a_{\rm
end}^{3}$) and the spectator field scales like a matter field, with
$\rho_{\chi}\propto a^{-3}$ and the amplitude of $\chi\propto a^{-3/2}$. We
also discuss this possibility below.

##  III Production Regimes

In this work, we focus on a broad range of $\tilde{\sigma}$ up to
$\tilde{\sigma}\lesssim 50$, where the broad parametric resonance effects
become significant and a full non-perturbative analysis is required
Garcia:2021iag ; Garcia:2022vwm . However, as shown in Garcia:2021iag , one
can treat the dark matter production perturbatively up to this value
(including the narrow resonance regime with
$\tilde{\sigma}\sim\mathcal{O}(1)$). Most importantly, we demonstrate that a
direct coupling between the inflaton and dark matter field opens up a broad
parameter space, allowing for a wide range of dark matter masses and reheating
temperatures, while successfully avoiding isocurvature constraints.666We note
that similar effects can be achieved by introducing a non-minimal coupling
between the dark matter field and gravity Garcia:2023qab .

###  III.1 Long Wavelength Contribution

When $\tilde{\sigma}\lesssim\mathcal{O}(1)$, the production of the spectator
field, $\chi$, is primarily dominated by the long wavelength contribution
(superhorizon modes). We assume that during inflation, the energy density of
the spectator field is largely subdominant, and inflation is driven solely by
the inflaton, $\phi$. At the end of inflation, we further assume that the
expectation value of the spectator field, $\chi$, has reached its asymptotic
value given by fluc

$\langle\chi^{2}\rangle_{\rm end}\;\simeq\frac{3H_{\rm
end}^{4}}{8\pi^{2}m^{2}_{\chi,\rm eff}(a_{\rm end})}\,.$ (10)

This result corresponds to the spectator field variance that is averaged over
long wavelength superhorizon modes, and we show the derivation of this result
using the stochastic approach in Appendix A.

Assuming that in a horizon patch, the field is approximately homogeneous, the
field value after inflation is given by $\chi_{\rm
end}=\sqrt{\langle\chi^{2}\rangle_{\rm end}}$. When $\sigma\phi^{2}(a)\gg
m_{\chi}^{2}$, one can show that the energy density of the spectator field
evolves as

$\rho_{\chi}(a)\;=\;\frac{1}{2}m^{2}_{\chi,\rm
eff}\chi^{2}(a)\;\simeq\;\frac{3H_{\rm
end}^{4}}{16\pi^{2}}\left(\frac{\phi(a)}{\phi_{\rm
end}}\right)^{2}\left(\frac{a_{\rm end}}{a}\right)^{\frac{3}{2}}\,.$ (11)

We note that because the inflaton field scales as $\phi(a)^{2}\propto a^{-3}$
and $\chi\propto a^{-3/2}$, this leads to an overall spectator field energy
density scaling of $\rho_{\chi}(a)\propto a^{-\frac{9}{2}}$. This means that
$\rho_{\chi}$ redshifts faster than $\rho_{\phi}$, and will not dominate the
energy budget of the Universe before reheating is complete.

We note that for large values of the effective mass, the fluctuations are cut
off exponentially with $\langle\chi^{2}\rangle\propto e^{-c\pi(m^{2}_{\chi,\rm
eff}/H_{\rm end}^{2}-\frac{9}{4})}$, where $c\sim\mathcal{O}(1)$ model-
dependent constant Chung:1998bt ; Markkanen:2016aes ; Chung:2018ayg ;
Racco:2024aac .777During inflation, the asymptotic value of
$\langle\chi^{2}\rangle$ is determined by $H_{I}$, but numerical studies
indicate that the subsequent evolution and density of $\chi$ is better fit
using $H_{\rm end}$ paper2 . Therefore, the exponential suppression becomes
significant when

$\frac{m^{2}_{\chi,\rm eff}(a_{\rm end})}{H_{\rm
end}^{2}}\gtrsim\frac{9}{4}\,,~{}~{}m^{2}_{\chi,\rm eff}(a_{\rm
end})\;=\;m^{2}_{\chi}+4\tilde{\sigma}H^{2}_{\rm end}\,,$ (12)

where we used $\phi^{2}_{\rm end}=4H_{\rm end}^{2}M_{P}^{2}/m_{\phi}^{2}$.
This implies that the exponential suppression needs to be accounted for when
$m_{\chi}^{2}\gtrsim\frac{9}{4}H_{\rm end}^{2}$ when $4\tilde{\sigma}H_{\rm
end}^{2}\ll m_{\chi}^{2}$, and $\tilde{\sigma}\gtrsim\frac{9}{16}$ when
$4\tilde{\sigma}H_{\rm end}^{2}\gg m_{\chi}^{2}$. Due to the exponential
suppression, the long wavelength contribution becomes subdominant for large
values of $\tilde{\sigma}\gg\frac{9}{16}$, and the contribution from the
gravitational scattering begins to dominate.

###  III.2 Gravitational Production After Inflation

When the direct coupling between the inflaton and spectator field becomes
large, i.e., $\tilde{\sigma}\gtrsim\frac{9}{16}$, the universal gravitational
interactions between the inflationary and dark sectors start dominating the
dark matter production. To compute this gravitational production of dark
matter during reheating, we expand the spacetime metric around the flat
Minkowski background, with
$g_{\mu\nu}\simeq\eta_{\mu\nu}+\frac{2}{M_{P}}h_{\mu\nu}$, where $h_{\mu\nu}$
is the canonically-normalized linear perturbation. Combining gravitational
interactions with the action in Eq. (1), we find the following interaction
terms (see e.g. hol )

$\mathcal{L}_{I}\;=\;-\frac{1}{M_{P}}h_{\mu\nu}\left(T_{\phi}^{\mu\nu}+T_{\chi}^{\mu\nu}\right)-\frac{1}{2}\sigma\phi^{2}\chi^{2}\,,$
(13)

where $T_{\phi}^{\mu\nu}$ and $T_{\chi}^{\mu\nu}$ are the energy-momentum
tensors of the inflaton and the spectator dark matter field, respectively,
given by

$T_{S=\phi,\,\chi}^{\mu\nu}\;=\;\partial^{\mu}S\partial^{\nu}S-g^{\mu\nu}\left[\frac{1}{2}\partial^{\alpha}S\partial_{\alpha}S-\frac{1}{2}m_{\rm{eff}}^{2}S^{2}\right]\,,$
(14)

where $m_{\rm eff}$ is the effective mass. The field perturbation $h_{\mu\nu}$
characterizes a massless spin-2 (canonically-normalized) graviton. Its
propagator in de Donder gauge that carries the momentum $k$ can be expressed
as Giudice:1998ck

$\Pi^{\mu\nu\rho\sigma}(k)\,=\,\dfrac{\eta^{\rho\nu}\eta^{\sigma\mu}+\eta^{\rho\mu}\eta^{\sigma\nu}-\eta^{\rho\sigma}\eta^{\mu\nu}}{2k^{2}}\,.$
(15)

The matrix element that corresponds to the final state of a pair of dark
matter particles, can be computed by treating the inflaton condensate as a
collection of particles, leading to dark matter production via the $s$-channel
scattering process MO ; cmov ; Clery:2022wib .

Importantly, because the effective mass of dark matter is
$\tilde{\sigma}$-dependent, it will affect the scattering process and lead to
kinematic blocking when $m_{\chi,\rm eff}^{2}\gtrsim m_{\phi}^{2}$, delaying
the scattering processes. However, even though the channel
$\phi\phi\rightarrow h_{\mu\nu}\rightarrow\chi\chi$ is sensitive to the
reheating temperature $T_{\rm RH}$ at the early stages of reheating cmov , as
the value of a direct coupling $\tilde{\sigma}$ increases, the scattering
contribution arising from the $\sigma\phi^{2}\chi^{2}$ coupling begins to
dominate dark matter production, leading to efficient particle production. We
discuss this in detail in the next section.

###  III.3 Model Parameters

We summarize the key parameters for the Starobinsky model of inflation with
$N_{*}=55$ $e$-folds that we use for our numerical results, along with the
field and energy density dependence on the scale factor $a$:

$\displaystyle m_{\phi}=1.25\times 10^{-5}M_{P}\,,$ $\displaystyle\phi_{\rm
end}=0.61\,M_{P}\,,$ (16) $\displaystyle H_{\rm end}=3.1\times
10^{-6}M_{P}\,,$ $\displaystyle|\dot{\phi}_{\rm end}|=\sqrt{2}H_{\rm
end}\,M_{P}\,,$ $\displaystyle\chi\propto a^{-\frac{3}{4}}\,,$
$\displaystyle\phi\propto a^{-\frac{3}{2}}\,,$
$\displaystyle\rho_{\chi}\propto a^{-\frac{9}{2}}\,,$
$\displaystyle\rho_{\phi}\propto a^{-3}\,.$

The scaling for $\chi$ applies when $a_{\rm osc}<a<a_{\chi}$, that is after
oscillations begin and before the effective mass is dominated by the bare
mass. Subsequently, $\rho_{\chi}\propto a^{-3}$ as its mass is constant.

As noted earlier, there are three important scale parameters: $a_{\rm osc}$,
$a_{\chi}$, and $a_{\rm RH}$. Their value relative to $a_{\rm end}$ can be
determined from their definitions: $m_{\chi,\rm eff}(a_{\rm
osc})=\frac{3}{2}H(a_{\rm osc})$, $\sigma\phi^{2}(a_{\chi})=m_{\chi}^{2}$, and
$\rho_{\phi}(a_{\rm RH})=\rho_{\rm end}(a_{\rm end}/a_{\rm RH})^{3}=\rho_{\rm
radiation}=\alpha T_{\rm RH}^{4}$, where $\alpha=g_{\rm RH}\pi^{2}/30$ and
$g_{\rm RH}=427/4$ is the number of Standard Model degrees of freedom.888This
assumes $T_{\rm RH}>m_{t}$ (the top quark mass). For lower $T_{\rm RH}$,
$g_{\rm RH}$ must be appropriately reduced. The solutions are given by

$\displaystyle\left(\frac{a_{\rm osc}}{a_{\rm end}}\right)^{3}$
$\displaystyle={\rm Max}\left(\frac{9H_{\rm
end}^{2}}{4m_{\chi}^{2}}\left(1-\frac{16}{9}{\tilde{\sigma}}\right),1\right)\,,$
(17) $\displaystyle\left(\frac{a_{\chi}}{a_{\rm end}}\right)^{3}$
$\displaystyle={\rm Max}\left(\frac{\sigma\phi_{\rm
end}^{2}}{m_{\chi}^{2}},1\right)={\rm Max}\left(\frac{4{\tilde{\sigma}}H_{\rm
end}^{2}}{m_{\chi}^{2}},1\right)\,,$ (18) $\displaystyle\left(\frac{a_{\rm
RH}}{a_{\rm end}}\right)^{3}$ $\displaystyle=\frac{3H_{\rm
end}^{2}M_{P}^{2}}{\alpha T_{\rm RH}^{4}}\,.$ (19)

Note that for $\tilde{\sigma}\gtrsim\frac{9}{16}$, the oscillations begin at
the end of inflation, $a_{\rm osc}=a_{\rm end}$. For smaller $\tilde{\sigma}$,
$(a_{\rm osc}/a_{\rm end})^{3}\simeq(1.1\times 10^{13}~{}{\rm
GeV}/m_{\chi})^{2}$. Similarly, $(a_{\chi}/a_{\rm
end})^{3}\simeq{\tilde{\sigma}}(1.5\times 10^{13}~{}{\rm GeV}/m_{\chi})^{2}$.

The order of these scale factors depends on $\sigma$, $m_{\chi}$, and $T_{\rm
RH}$, assuming that $\phi_{\rm end}$ and $H_{\rm end}$ are set by the
inflationary model. Thus, we have

$\displaystyle\qquad{\tilde{\sigma}}<\frac{9}{32}\,,$ $\displaystyle
a_{\chi}<a_{\rm osc}\,,$ (20)
$\displaystyle\qquad{\tilde{\sigma}}<\frac{3}{4}\frac{m_{\chi}^{2}M_{P}^{2}}{\alpha
T_{\rm RH}^{4}}\,,$ $\displaystyle a_{\chi}<a_{\rm RH}\,,$ (21)
$\displaystyle\qquad{\tilde{\sigma}}<\frac{9}{16}\left(1-\frac{4M_{P}^{2}m_{\chi}^{2}}{3\alpha
T_{\rm RH}^{4}}\right)\,,$ $\displaystyle a_{\rm RH}<a_{\rm osc}\,.$ (22)

Each of these cases implies a different evolution for $\rho_{\chi}$. We will
discuss them one by one in the following.

###  III.4 Isocurvature Constraints

We also discuss the isocurvature constraints. The current isocurvature power
spectrum constraints from Planck are given by

$\beta_{\rm
iso}\;\equiv\;\frac{\mathcal{P}_{\mathcal{S}}(k_{*})}{\mathcal{P}_{\mathcal{R}}(k_{*})+\mathcal{P}_{\mathcal{S}}(k_{*})}<0.038\,,$
(23)

at $95\%$ CL with the pivot scale $k_{*}=0.05\,\rm{Mpc}^{-1}$ Planck:2018jri ,
where $\mathcal{P}_{\mathcal{R}}$ represents the curvature spectra and
$\mathcal{P}_{\mathcal{S}}$ denotes the isocurvature spectra. During
inflation, if dark matter is light compared to the inflationary scale, it
leads to the production of isocurvature modes in tension with the current
Planck constraints Chung:2004nh ; Garcia:2023awt .

The curvature spectrum of a massive scalar field can be expressed as

$\mathcal{P}_{\mathcal{R}}(k)\;=\;\frac{H_{I}^{2}}{4\pi^{2}}\left(\frac{k}{aH_{I}}\right)^{3-2\sqrt{\frac{9}{4}-\frac{m^{2}_{\chi,\rm
eff}}{H_{I}^{2}}}}\,,$ (24)

when $m^{2}_{\chi,\rm eff}/H_{I}^{2}\lesssim 9/4$, and the isocurvature power
spectrum is given by

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

where $\rho_{\chi}$ is the dark matter density and $\delta\rho_{\chi}$ is the
dark matter fluctuation. Because the dark matter energy density scales as
$\rho_{\chi}(x)\propto\chi^{2}(x)$, and assuming Gaussianity, we can use the
approximation

$\frac{\langle\delta\rho_{\chi}(x)\delta\rho_{\chi}(0)\rangle}{\langle\rho_{\chi}\rangle^{2}}=\frac{\left\langle\chi^{2}(x)\chi^{2}(0)\right\rangle-\left\langle\chi^{2}\right\rangle^{2}}{\left\langle\chi^{2}\right\rangle^{2}}=2\frac{\langle\chi(x)\chi(0)\rangle^{2}}{\left\langle\chi^{2}\right\rangle^{2}}\,.$
(26)

This allows us to recast the isocurvature power spectrum in the following form

$\displaystyle\mathcal{P}_{\mathcal{S}}(k)$
$\displaystyle=\frac{k^{3}}{2\pi^{2}}\frac{2}{\left\langle\chi^{2}\right\rangle^{2}}\int
d^{3}xe^{-i\boldsymbol{k}\cdot\boldsymbol{x}}\langle\chi(x)\chi(0)\rangle^{2}$
(27)
$\displaystyle=\frac{k^{3}}{2\pi^{2}}\frac{2}{\left\langle\chi^{2}\right\rangle^{2}}\int\frac{d^{3}q}{(2\pi)^{3}}P(q)P(k-q)\,.$

If we directly compute this integral, we obtain Redi:2022zkt

$\mathcal{P}_{\mathcal{S}}(k)\;\simeq\;\frac{8m_{\chi,\rm eff}^{2}}{3H_{\rm
end}^{2}}\left(\frac{k}{a_{\rm end}H_{I}}\right)^{\frac{4m_{\chi,\rm
eff}^{2}}{3H_{I}^{2}}}\,.$ (28)

Using this expression, we observe that the modes that exit the horizon early
during inflation will be strongly suppressed. Therefore, for the $k_{*}$ mode
that exits after $N_{*}=55$ $e$-folds of inflation, we can use the
approximation

$\mathcal{P}_{\mathcal{S}}(k_{*})\;\simeq\;\frac{8m_{\chi,\rm
eff}^{2}}{3H_{I}}\exp\left(\frac{-4N_{*}m^{2}_{\chi,\rm
eff}}{3H_{I}^{2}}\right)\,.$ (29)

If we now use the $N_{*}=55$ $e$-folds, we find the constraint

$m_{\chi,\rm eff}(t_{*})\gtrsim 0.5H_{I}\simeq H_{\rm end}\,,$ (30)

where we used the fact that for the Starobinsky model of inflation,
$H_{I}\simeq 2H_{\rm end}$. This constraint implies that

$\displaystyle m_{\chi}$ $\displaystyle\gtrsim H_{\rm
end}\,,\qquad\sigma\phi_{*}^{2}\ll m^{2}_{\chi,\rm eff}\,,$ (31)
$\displaystyle\tilde{\sigma}$ $\displaystyle\gtrsim
0.003\,,\qquad\sigma\phi_{*}^{2}\gg m^{2}_{\chi,\rm eff}\,,$ (32)

where in the second line, we used the ratio $\phi_{*}^{2}/\phi_{\rm
end}^{2}\simeq 75$ for the Starobinsky model of inflation.

A fully numerical computation of the isocurvature constraint for the
Starobinsky model of inflation, that was performed in Ref. Garcia:2023awt ,
shows that the isocurvature constraints are given by

$\displaystyle m_{\chi}$ $\displaystyle\gtrsim 1.1H_{\rm end}\,,$
$\displaystyle\sigma\phi_{*}^{2}\ll m^{2}_{\chi,\rm eff}\,,$ (33)
$\displaystyle\tilde{\sigma}$ $\displaystyle\gtrsim 0.004\,,$
$\displaystyle\sigma\phi_{*}^{2}\gg m^{2}_{\chi,\rm eff}\,,$ (34)

which are in excellent agreement with the analytical results. We use the fully
numerical results when discussing the dark matter constraints. However, the
analytical procedure carried out in this section is general and can be easily
applied to various models of inflation.

##  IV Scenarios

###  IV.1 Case I: $\sigma\phi_{\rm end}^{2}\ll m_{\chi}^{2}$

When $\sigma$ is very small (or absent), the bare mass term comes to dominate
early (and will then always dominate at later times), and $a_{\chi}<a_{\rm
osc},a_{\rm RH}$. The spectator field starts oscillating when
$\frac{3}{2}H(a_{\rm osc})\simeq m_{\chi}$. Because the Hubble parameter
scales as $H(a)\propto a^{-\frac{3}{2}}$, the oscillations start before
reheating if $\frac{3}{4}\frac{\rho_{\rm RH}}{M_{P}^{2}}<m_{\chi}^{2}$, or
$3\alpha T_{\rm RH}^{4}<4m_{\chi}^{2}M_{P}^{2}$, when $a_{\rm osc}/a_{\rm
end}=(3H_{\rm end}/2m_{\chi})^{2/3}$. We refer to this case as IA, where
$a_{\rm\chi}<a_{\rm osc}<a_{\rm RH}$.

The initial value when the field starts oscillating is given by $\chi(a_{\rm
osc})=\sqrt{\langle\chi^{2}\rangle_{\rm end}}$. We want to compute the energy
density of the spectator field at the time of reheating. By combining Eqs.
(10), (17) and (19), we find

$\displaystyle\rho_{\chi}(a_{\rm
RH})\;=\;\frac{1}{2}m_{\chi}^{2}\langle\chi^{2}\rangle_{\rm
end}\left(\frac{a_{\rm osc}}{a_{\rm end}}\right)^{3}\left(\frac{a_{\rm
end}}{a_{\rm RH}}\right)^{3}$ (35)
$\displaystyle\;=\;\frac{9}{64\pi^{2}}H_{\rm end}^{4}\frac{\alpha T_{\rm
RH}^{4}}{M_{P}^{2}({m_{\chi}^{2}+\sigma\phi_{\rm
end}^{2}})}\left(1-\frac{16}{9}\tilde{\sigma}\right)$
$\displaystyle\;=\;\frac{9\alpha H_{\rm end}^{4}T_{\rm
RH}^{4}}{64\pi^{2}m_{\chi}^{2}M_{P}^{2}}\,,\qquad\sigma\phi_{\rm end}\ll
m_{\chi}^{2};\,{\tilde{\sigma}\ll 1}\,.$

The numerical evolution of $\rho_{\chi}$ is shown in Fig. 1 when multiplied by
$a^{3}$. The choice of parameters is given in the caption and have been chosen
to yield a present density of $\Omega_{\chi}h^{2}=0.12$. The dots indicate the
values of $a_{\chi}=1$, $a_{\rm osc}=5$. Since the reheating temperature in
this case is low, $a_{\rm RH}=\mathcal{O}(10^{18})$, far off the scale of the
plot (values of $a$ are all relative to $a_{\rm end}$). The star shows the
result from Eq. (35) and is in excellent agreement with the numerical result.

Figure 1:  Energy density evolution of the spectator field $\rho_{\chi}a^{3}$.
Here we assumed $\sigma\phi_{\rm end}^{2}\ll m_{\chi}^{2}$, $m_{\chi}=10^{12}$
GeV, and $T_{\rm RH}=77$ GeV. The star shows the relic density obtained from
Eq. (35)

To obtain the present-day relic abundance of dark matter, the energy density
in Eq. (35) must be redshifted to the present by

$\left(\frac{a_{\rm RH}}{a_{0}}\right)^{3}=\left(\frac{T_{0}}{T_{\rm
RH}}\right)^{3}\frac{g_{0}}{g_{\rm RH}}\,,$ (36)

where $T_{0}$ is the present temperature of the CMB, and $g_{0}=(43/4)(4/11)$.
The fraction of critical density at the present time is then given by mybook

$\displaystyle\frac{\Omega_{\chi}h^{2}}{0.12}\;=\;\frac{\rho_{\chi}(a_{0})}{\rho_{c}}~{}\frac{h^{2}}{0.12}\;\simeq
10^{47}\frac{\rho_{\chi}(a_{0})}{\rm GeV^{4}}$ (37)
$\displaystyle\simeq\;4.9\times 10^{7}~{}{\rm
GeV}^{-1}~{}\frac{\rho_{\chi}(a_{\rm RH})}{T_{\rm RH}^{3}}\,,$

where we used

$\rho_{c}=8.1\times 10^{-47}\,h^{2}\,{\rm GeV^{4}}=1.054\times
10^{-5}\,h^{2}\,\rm{GeV}\,\rm{cm}^{-3}\,.$ (38)

Combining these results, we find

$\displaystyle\frac{\Omega_{\chi}h^{2}}{0.12}\;\simeq\;2.5\times
10^{7}\,\mathrm{GeV}^{-1}\frac{H_{\rm end}^{4}T_{\rm
RH}}{M_{P}^{2}({m_{\chi}^{2}+\sigma\phi_{\rm
end}^{2}})}\left(1-\frac{16}{9}\tilde{\sigma}\right)$
$\displaystyle\;\simeq\;\frac{T_{\rm RH}}{77~{}{\rm
GeV}}\left(\frac{10^{12}~{}{\rm
GeV}}{m_{\chi}}\right)^{2}\,,\qquad\sigma\phi_{\rm end}\ll
m_{\chi}^{2};\,{\tilde{\sigma}\ll 1}\,.$ (39)

Figure 2: $T_{\rm RH}$ vs $\tilde{\sigma}$ satisfying
$\Omega_{\chi}h^{2}=0.12$ for $m_{\chi}=10^{2}$(red), $10^{4}$ (orange),
$10^{7}$ (green), $10^{10}$ (blue), $10^{12}$ (purple) GeV. The gray dashed
line on the left panel shows the BBN bound on $T_{\rm RH}$. The cyan shaded
region in the left panel is excluded by the isocurvature constraint in Eq.
(34). The right panel magnifies the range $10^{-2}\leq\tilde{\sigma}\leq 50$
of the left panel.

The solutions to $\Omega_{\chi}h^{2}=0.12$ are shown by the horizontal
portions of the curves in Fig. 2. In this regime, the upper limit in the
reheating temperature (to ensure that $\Omega_{\chi}h^{2}\leq 0.12$) from Eq.
(39) is $T_{\rm RH}\lesssim 7.7\times 10^{-{23}}m_{\chi}^{2}$ (all quantities
in GeV). Thus, even a spectator mass equal to the inflaton mass would require
a reheating temperature $\lesssim 69$ TeV, using $H_{\rm end}$ given in Eq.
(16). This is consistent with the results found in Cosme:2020nac . Note also
the constraint $T_{\rm RH}\gtrsim 4$ MeV to ensure successful nucleosynthesis
tr4  imposes a lower bound on $m_{\chi}\gtrsim 7\times 10^{9}$ GeV.

In addition to the contribution to the relic density from the oscillations of
$\chi$, there are two additional sources of $\chi$ production. One is the
single graviton exchange process involving particle production from the
inflaton condensate. The other is the four-point interaction. However, for
small $\tilde{\sigma}$, the relic density from long wavelength far exceeds
that obtained from particle production through single graviton exchange.

The rate per unit volume for gravitational production cmov ; Clery:2022wib ;
Mambrini:2022uol  including the four-point interaction is given by

$\displaystyle R_{g}=$ $\displaystyle\frac{2\times\rho_{\phi}^{2}}{256\pi
M_{P}^{4}}\left(1-4\tilde{\sigma}+\frac{m_{\chi,{\rm
eff}}^{2}}{2m^{2}_{\phi}}\right)^{2}\sqrt{1-\frac{m_{\chi,{\rm
eff}}^{2}}{m_{\phi}^{2}}}\,,$ (40)

where the explicit factor of 2 indicates that two $\chi$ particles are
produced per interaction. To obtain the relic density, this rate should be
integrated in the Boltzmann equation,

$\frac{dY_{\chi}}{da}=\frac{\sqrt{3}M_{P}}{\sqrt{\alpha}T_{\rm
RH}^{2}}a^{2}\left(\frac{a}{a_{\rm RH}}\right)^{\frac{3}{2}}R_{g}(a)\,,$ (41)

where the energy density of $\chi$ is $\rho_{\chi}=m_{\chi,{\rm
eff}}Y_{\chi}/a^{3}$. When $m_{\chi,{\rm eff}}^{2}\ll m_{\phi}^{2}$ (and
${\tilde{\sigma}}\ll 1$), there is a simple analytic solution cmov , and the
relic density can be approximated by

$\frac{\Omega_{\chi}h^{2}}{0.12}\simeq 8.6\times 10^{6}\frac{H_{\rm end}T_{\rm
RH}m_{\chi}}{{\rm GeV}M_{P}^{2}}\,.$ (42)

This becomes comparable to the density given in Eq. (39) only when
$m_{\chi}\simeq 1.1\times 10^{13}$ GeV. Note that the relic density from
single graviton exchange is proportional to $m_{\chi}$, whereas the density
from large-scale fluctuations is inversely proportional to the mass squared.
Combining both sources of production, we see that there is an additional an
upper limit to $T_{\rm RH}\lesssim 4.3$ TeV. Thus, in the absence of even a
small (Planck-suppressed) coupling of the spectator to the inflaton, the
reheating temperature is restricted to a range $0.004<T_{\rm RH}/{\rm
GeV}<4300$, and the mass of the spectator is then restricted to the range
$7\times 10^{9}<m_{\chi}/{\rm GeV}<3\times 10^{13}$ (the upper limit
corresponding to the inflaton mass and hence kinematic limit). As we will see,
the gravitational strength coupling, when included, greatly opens up the
allowed range for the reheating temperature and spectator mass.

The spectator field may begin oscillating after reheating ($a_{\rm osc}>a_{\rm
RH}$ corresponding to $3\alpha T_{\rm RH}^{4}>m_{\chi}^{2}M_{P}^{2}$). We
refer to this as case IB. In this case, we can evolve $\rho_{\chi}$ directly
from $a_{\rm osc}$ to present time using

$\displaystyle\rho_{\chi}(a_{0})\;=\;\rho_{\chi}(a_{\rm
osc})\left(\frac{a_{\rm osc}}{a_{0}}\right)^{3}$ (43)
$\displaystyle\;=\;\frac{3}{16\pi^{2}}H_{\rm
end}^{4}\frac{m_{\chi}^{2}}{m_{\chi}^{2}+\sigma\phi_{\rm
end}^{2}}\left(\frac{g_{0}}{g_{\rm osc}}\right)\left(\frac{T_{0}}{T_{\rm
osc}}\right)^{3}$ $\displaystyle\;\simeq\;\frac{3}{16\pi^{2}}H_{\rm
end}^{4}\frac{(3\alpha)^{\frac{3}{4}}T_{0}^{3}}{2^{\frac{3}{2}}M_{P}^{\frac{3}{2}}m_{\chi}^{\frac{3}{2}}}\left(\frac{g_{0}}{g_{\rm
osc}}\right)\,,$

when $a_{\rm\chi}<a_{\rm RH}<a_{\rm osc}$ and is independent of the reheating
temperature. To obtain the last line of Eq. (43), we used $m_{\chi,{\rm
eff}}^{2}(a_{\rm osc})\simeq m_{\chi}^{2}=\frac{9}{4}H^{2}=\frac{3}{4}\alpha
T_{\rm osc}^{4}/M_{P}^{2}$.

The present fraction of critical density is then given by Eq. (37)

$\displaystyle\frac{\Omega_{\chi}h^{2}}{0.12}\;=\;1.1\times 10^{7}{\rm
GeV}^{-1}\frac{H_{\rm end}^{4}}{M_{P}^{\frac{3}{2}}m_{\chi}^{\frac{3}{2}}}$
(44) $\displaystyle\simeq\left(\frac{4.3\times 10^{20}~{}{\rm
GeV}}{m_{\chi}}\right)^{\frac{3}{2}}\,,~{}~{}~{}~{}\sigma\phi_{\rm end}^{2}\ll
m_{\chi}^{2}\,.$

It is then clear that the Universe is largely over-closed in this case, as the
condition $\Omega h^{2}\lesssim 0.12$ requires $m_{\chi}>M_{P}$. We do not
consider this case any further.

###  IV.2 Case II: $\sigma\phi_{\rm end}^{2}\gg m_{\chi}^{2}$

As $\sigma$ is increased, the dependence on $m_{\chi}$ in Eq. (35) disappears
(as long as $m_{\chi}^{2}\ll\sigma\phi_{\rm end}^{2}$) and $T_{\rm RH}$
increases. As long as ${\tilde{\sigma}}<9/32$, we maintain, $a_{\rm
end}<a_{\chi}<a_{\rm osc}$ and the relic density (39) becomes

$\frac{\Omega_{\chi}h^{2}}{0.12}\;\simeq\;2.5\times
10^{7}\,\mathrm{GeV}^{-1}\frac{H_{\rm end}^{4}T_{\rm
RH}}{M_{P}^{2}\sigma\phi_{\rm end}^{2}}\simeq\tilde{\sigma}^{-1}\frac{T_{\rm
RH}}{25~{}{\rm TeV}}\,.$ (45)

Thus for a fixed relic density, we require an increase in $T_{\rm RH}\sim
25{\tilde{\sigma}}$ GeV as seen in Fig. 2. The curves depart from being
horizontal at low ${\tilde{\sigma}}$ when ${\tilde{\sigma}}\simeq
m_{\chi}^{2}/4H_{\rm end}^{2}$. They then appear to merge (due to the log
scale of the plot) as $T_{\rm RH}$ increases with $\tilde{\sigma}$ almost
independently of $m_{\chi}$. This behavior is very similar to Eq. (39), with
$\tilde{\sigma}H_{\rm end}^{2}$ playing the role of $m_{\chi}^{2}$.

For ${\tilde{\sigma}}>9/32$, the spectator field begins to oscillate while
dominated by the effective mass $\propto\phi$, until $a=a_{\chi}$, where the
oscillations are dominated by the bare mass. Thus, $a_{\rm osc}<a_{\chi}$ (see
Eq. (20)). For lower values of $\tilde{\sigma}$, $a_{\chi}<a_{\rm osc}$, and
case I is applicable. Note that, as $\sigma\phi^{2}$ redshifts as $a^{-3}$, as
does $H^{2}$, once the oscillations of $\chi$ begin ($m_{\chi,{\rm
eff}}>\frac{3}{2}H$) they never stop . Indeed, the effective mass term
$\propto\phi$ decreases with time, and the bare mass dominates the
oscillations at $a=a_{\chi}$. It is easy to understand that for larger
$\tilde{\sigma}$ the phase of oscillations begin earlier because
$\sigma\phi^{2}>m_{\chi}$. During this phase of oscillations, $\rho_{\chi}$
undergoes a pre-dilution between $a_{\rm end}$ and $a_{\chi}$, opening the
$T_{\rm RH}$ parameter space relative to that allowed when only the bare mass
term is present.

Similarly to case IA, we will consider first the late reheating scenario,
i.e., $a_{\chi}<a_{\rm RH}$, or $\sigma\phi^{2}(a_{\rm RH})\lesssim
m_{\chi}^{2}$, with the upper limit on $\tilde{\sigma}$ given by Eq. (21). We
refer to this as case IIA, which allows for much larger $T_{\rm RH}$ than case
IA.

The energy density must now be evolved from $a_{\rm end}$ as follows:

$\displaystyle\rho_{\chi}(a_{0})\;=$
$\displaystyle\frac{3\beta}{16\pi^{2}}H_{\rm end}^{4}\left(\frac{a_{\rm
end}}{a_{\rm osc}}\right)^{3}\left(\frac{a_{\rm osc}}{a_{\chi}}\right)^{3+2A}$
(46) $\displaystyle\times\left(\frac{a_{\chi}}{a_{\rm
RH}}\right)^{3}\left(\frac{a_{\rm RH}}{a_{0}}\right)^{3}$ $\displaystyle=$
$\displaystyle\frac{3}{16\pi^{2}}H_{\rm end}^{4}\left(\frac{a_{\rm
osc}}{a_{\rm end}}\right)^{2A}\left(\frac{a_{\rm end}}{a_{\chi}}\right)^{2A}$
$\displaystyle\;\times\;\left(\frac{a_{\rm end}}{a_{\rm
RH}}\right)^{3}\left(\frac{a_{\rm RH}}{a_{0}}\right)^{3}\,,$

where $A$ is defined from Eq. (8) For $\tilde{\sigma}=\frac{9}{32}$, using the
relations in Eqs. (17)-(19), this solution matches that given in Eq. (37). For
the limited range $9/32<{\tilde{\sigma}}<9/16$, $a_{\rm end}<a_{\rm osc}$. For
larger ${\tilde{\sigma}}>9/16$, $a_{\rm osc}=a_{\rm end}$. However, in this
case, the exponential suppression in the initial value of
$\langle\chi^{2}\rangle$ alluded to earlier,

$\beta={\rm
Min}\left(\exp\left[{-1.25\pi\left(4{\tilde{\sigma}}+\frac{m_{\chi}^{2}}{H_{\rm
end}^{2}}-\frac{9}{4}\right)}\right],1\right)\,,$ (47)

must be included.

In this case, the resulting density fraction is

$\Omega_{\chi}h^{2}\;=\;6.7\times 10^{5}~{}\beta~{}\frac{H_{\rm
end}m_{\chi}T_{\rm RH}}{\mathrm{GeV}M_{P}^{2}\sqrt{\tilde{\sigma}}}\,.$ (48)

or

$\frac{\Omega_{\chi}h^{2}}{0.12}\simeq\frac{\beta}{\sqrt{\tilde{\sigma}}}\left(\frac{T_{\rm
RH}}{10^{10}~{}{\rm GeV}}\right)\left(\frac{m_{\chi}}{1.4\times 10^{7}~{}{\rm
GeV}}\right)\,.$ (49)

The numerical solution to the equation of motion for $\chi$ in case IIA with
${\tilde{\sigma}}=1.7$ is shown in Fig. 3. For this value of $\tilde{\sigma}$,
$a_{\rm osc}=a_{\rm end}$, $a_{\chi}\simeq 6.1\times 10^{5}a_{\rm end}$, and
$a_{\text{RH}}\simeq 1.3\times 10^{6}a_{\rm end}$. In contrast to the case
shown in Fig. 1, oscillations begin when inflation ends and the density
decreases as $a^{-\frac{9}{2}}$. For $a>a_{\chi}$, the density drops as
$a^{-3}$. Once again, we see excellent agreement between the numerical and
analytical results.

Figure 3:  The evolution of the energy density of the spectator field
$\rho_{\chi}a^{3}$ for case IIA (black), case IIB with $m_{\chi,{\rm
eff}}>2H(a_{\rm RH})$ (red), and case IIB with $m_{\chi,{\rm eff}}<2H(a_{\rm
RH})$ (blue). $a_{2}$ denotes $a_{\text{osc-2}}$. Parameters chosen are:
$\tilde{\sigma}=1.7$, $m_{\chi}=5\times 10^{4}$ GeV, and $T_{\rm
RH}=5.95\times 10^{10}$ GeV for the black curve; $\tilde{\sigma}=1.8$,
$m_{\chi}=3\times 10^{4}$ GeV, and $T_{\rm RH}=8.7\times 10^{10}$ GeV for the
red curve; and $\tilde{\sigma}=1.5$, $m_{\chi}=1.5\times 10^{4}$ GeV, and
$T_{\rm RH}=2.64\times 10^{11}$ GeV (blue curve). In each case, the star
denotes the present day relic density from Eqs. (48), (50), and (51),
respectively. In all three cases, $a_{\rm osc}=a_{\rm end}$. The parameters in
each case are chosen so that $\Omega_{\chi}h^{2}=0.12$ when both production
mechanisms are considered.

As one can discern from Eq. (49), the allowed range for $m_{\chi}$ and $T_{\rm
RH}$ is considerably greater than that found for case IA, when the coupling
$\sigma$ could be ignored. In this case, we also see that the relic density is
proportional to $m_{\chi}T_{\rm RH}$, whereas at small ${\tilde{\sigma}}$, the
relic density of $\chi$ was dominated by the large-scale fluctuations, as
${\tilde{\sigma}}$ is increased, the contribution from particle production
(either from single graviton exchange or directly from the contact term) comes
to dominate. Recall that for small $\tilde{\sigma}$, the two contributions
become comparable only at large bare mass $m_{\chi}\simeq 10^{13}$ GeV. When
$\tilde{\sigma}>9/16$, the density produced from single graviton exchange and
the contact interaction begin to dominate for all masses. Indeed, for
$\tilde{\sigma}>9/16$ there are two sources of suppression of the large
wavelength modes. One is the exponential suppression when the effective mass
is large and this occurs when ${\tilde{\sigma}}=9/16-m_{\chi}^{2}/4H_{\rm
end}^{2}$. However, as noted earlier, when ${\tilde{\sigma}}=9/16$, $a_{\rm
osc}=a_{\rm end}$, and we lose a strong enhancement of $(a_{\rm osc}/a_{\rm
end})^{3/2}\gg 1$ at smaller $\tilde{\sigma}$.

These effects are displayed in both panels of Fig. 2. At
$\tilde{\sigma}=9/16$, we see a meteoric rise in the required reheating
temperature to attain $\Omega_{\chi}h^{2}=0.12$. The degree to which the
temperature rises depends on the bare mass and is roughly inversely
proportional to $m_{\chi}$. At large $\tilde{\sigma}$, the relic density is
predominantly due to particle production and obtained from a numerical
integration of Eq. (41). Because $m_{\chi,{\rm
eff}}^{2}=m_{\chi}^{2}+\sigma\phi^{2}(a)$ is a function of the scale factor
$a$, a useful analytical form for $\rho_{\chi}$ is not available. Furthermore,
as ${\tilde{\sigma}}$ is increased, the effective mass of $\chi$ increases,
and eventually the particle production processes are kinematically forbidden
when inflation ends. However, because the effective mass depends on $a$, as
$a$ increases, particle production is again kinematically allowed. This effect
first occurs when $\tilde{\sigma}\simeq 2.5$ and is visible in the plot as a
small change in slope at $\log{\tilde{\sigma}}\simeq 0.4$.

We now see the full effect of adding the interaction term for $\chi$. In its
absence, we are restricted to very large scalar masses ($10^{10}-10^{13}$) GeV
and very low reheating temperatures ($0.004-4300)$ GeV as seen in the
horizontal portions of the curves in the left panel of Fig. 2. For non-
negligible ${\tilde{\sigma}}\leq 9/16$, the reheating temperature may be as
large as the allowed maximum
$=(3/\alpha)^{\frac{1}{4}}(H_{I}M_{P})^{\frac{1}{2}}=2.3\times 10^{15}$ GeV,
where the maximal value is attained for low masses $m_{\chi}\simeq 50$ GeV. At
${\tilde{\sigma}}=9/16$, the contribution of the large-scale fluctuations is
highly suppressed and the relic density is produced by particle production
from the inflaton condensate. The inclusion of a coupling between the inflaton
and the spectator now allows for a wide range of masses and reheating
temperatures all satisfying $\Omega_{\chi}h^{2}=0.12$.

In Fig. 4, we show the allowed parameter space for ${\tilde{\sigma}}\geq 9/16$
in the $(T_{\rm RH},m_{\chi})$ plane. The curves correspond to
$\Omega_{\chi}h^{2}=0.12$ for fixed values of $\tilde{\sigma}$ as indicated.
As noted above, for ${\tilde{\sigma}}\geq 9/16$, the relic density is
saturated by the contribution from inflaton scattering. Here, we see that the
temperature is maximal at low masses and that this maxiumum occurs at lower
$m_{\chi}$ as $\tilde{\sigma}$ is increased. The dotted line in Fig. 4 shows
the result from cmov  due to inflaton scattering through single graviton
exchange (i.e., with $\tilde{\sigma}=0$) neglecting the (dominant)
contribution from large scale fluctuations. In this case there is direct
relation between $m_{\chi}$ and $T_{\rm RH}$ as given in Eq. (42). For masses
lower than $\sim 100$ GeV, there are no solutions that yield the correct relic
density (smaller masses and/or reheating temperatures give
$\Omega_{\chi}h^{2}<0.12$). As $\tilde{\sigma}$ is increased, the lower limit
on the mass decreases. For example, for ${\tilde{\sigma}}=1$, $m_{\chi}\gtrsim
35$ GeV, the maximal temperature cannot be obtained as kinematic blocking due
to the large effective mass prevents particle production, and so the reheating
temperature must be lower so that $a_{\rm RH}$ is large enough to allow a
finite range of integration in Eq. (41). In this case, $m_{\chi}\gtrsim 10$
GeV and the maximum reheating temperature is $T_{\rm RH}\sim 1.1\times
10^{15}$. For ${\tilde{\sigma}}=50$, $m_{\chi}>0.03$ GeV and $T_{\rm RH}\simeq
8\times 10^{14}$ GeV is required to obtain $\Omega_{\chi}h^{2}=0.12$.

Figure 4: The reheating temperature, $T_{\rm RH}$ vs the bare mass,
$m_{\chi}$, satisfying $\Omega_{\chi}h^{2}=0.12$ for
$\tilde{\sigma}=\frac{9}{16}$(black), 1 (purple), $5$(blue), $10$(red),
$50$(orange), $0$(gray dashed). The line for $\tilde{\sigma}=0$ neglects the
contribution from large-scale fluctuations and corresponds to the result in
cmov . For $m_{\chi}<m_{\rm\chi,min}$, (corresponding to the endpoint of the
lines) there is no $T_{\rm RH}$ satisfying $\Omega_{\chi}h^{2}=0.12$.

Included in the solutions shown in Figs. 2 and 4, are solutions corresponding
to case IIA as well as other evolutionary patterns. For completeness, we
describe these possibilities. When ${\tilde{\sigma}}>\frac{9}{32}$ and
${\tilde{\sigma}}>3m_{\chi}^{2}M_{P}^{2}/4\alpha T_{\rm RH}^{4}$, $a_{\rm
osc}<a_{\rm RH}<a_{\chi}$ and we refer to this as case IIB. As before, for
${\tilde{\sigma}}>\frac{9}{16}$, $a_{\rm osc}=a_{\rm end}$ (for smaller
${\tilde{\sigma}}$ there is an additional suppression that can be obtained
from Eq. (17)). However, when $a_{\chi}>a_{\rm RH}$, the effective mass after
reheating quickly becomes equal to the bare mass as the inflaton decays. The
spectator density can then be computed using
$m_{\chi}^{2}\langle\chi^{2}\rangle$ and allowing $\chi$ to evolve as
$a^{-2A}$ and providing a factor $(a_{\rm end}/a_{\rm RH})^{3/2}$. After
reheating (or $a>a_{\chi}$ which is only slightly larger than $a_{\rm RH}$)
the density falls as $a^{-3}$ providing the factor in Eq. (36). Putting all
the factors together, we arrive at

$\frac{\Omega_{\chi}h^{2}}{0.12}\;=\;\frac{3.2\times 10^{6}H_{\rm
end}^{3}m_{\chi}^{2}\beta}{\mathrm{GeV}M_{P}T_{\rm
RH}(m_{\chi}^{2}+\sigma\phi_{\rm end}^{2})}\,.$ (50)

The evolution of $\rho_{\chi}$ for this case with $a_{\rm osc}<a_{\rm
RH}<a_{\chi}$ is shown by the red curve in Fig. 3. Note that due to inflaton
decay, the value of $a_{\chi}$ is slightly larger than that obtained from Eq.
(18), though this is not used in the numerical integration. The analytic
result from Eq. (50) gives very good agreement with the numerical result.

In this case, we must further distinguish an additional possibility. When
$m_{\chi,{\rm eff}}>2H(a_{\rm RH})$ oscillations of $\chi$ continue
uninterrupted.999As the Universe is radiation-dominated, the condition for
oscillations is $m=2H$ rather than $m=3H/2$ as used earlier in the matter-
dominated era. This led to the solution given in Eq. (50). If, however,
$m_{\chi,{\rm eff}}$ is small and $m_{\chi,{\rm eff}}<2H(a_{\rm RH})$ after
reheating, there will be a delay in the oscillations of $\chi$, until $a_{\rm
osc-2}$ when $m_{\chi,{\rm eff}}=2H(a_{\rm osc-2})$. Relative to the previous
result in Eq. (50), there is a slight enhancement by a factor $(a_{\rm
osc-2}/a_{\rm RH})^{3}=(T_{\rm RH}/T_{\rm osc-2})^{3}$, where $T_{\rm osc-2}$
is given by a similar expression used in Eq. (43), $4\alpha T_{\rm
osc-2}^{4}=3m_{\chi}^{2}M_{P}^{2}$. The result is

$\frac{\Omega_{\chi}h^{2}}{0.12}\;=\;\frac{5.7\times 10^{7}H_{\rm
end}^{3}\sqrt{m_{\chi}}T_{\rm
RH}^{2}\beta}{\mathrm{GeV}M_{P}^{5/2}(m_{\chi}^{2}+\sigma\phi_{\rm
end}^{2})}\,.$ (51)

The evolution of $\rho_{\chi}$ for this case is shown by the blue curve in
Fig. 3. Note that for $a_{\rm osc-2}>a>a_{\chi}$, the evolution of $\chi$
stops and $\rho_{\chi}$ is constant, leading to the observed increase in
$a^{3}\rho_{\chi}$ seen in the figure. Once again, for this parameter set, the
agreement between the numerical evaluation and the analytical result is
excellent.

##  V Discussion

In this work, we have considered the production of a massive (stable) scalar
field during inflation and the subsequent reheating process. It is well known
fluc  that for scalar fields with masses significantly below the Hubble scale
during inflation, large-scale fluctuations are produced with
$\langle\chi^{2}\rangle\gg H_{I}^{2}$ as given in Eq. (10). When inflation
ends, these fluctuations can be viewed as long-wavelength oscillations, and
their contribution to the energy density takes the form of matter, which can
account for all or part of the dark matter. Indeed, for light fields, this
form of production is so efficient that there are strong constraints on the
mass and reheating temperature to avoid the overproduction of dark matter
Turner:1987vd ; Peebles:1999fz ; Enqvist:2014zqa ; Nurmi:2015ema ;
Bertolami:2016ywc ; Alonso-Alvarez:2018tus ; Markkanen:2018gcw ;
Tenkanen:2019aij ; Cosme:2020nac ; Lebedev:2022cic , see Eq. (39). This can be
seen graphically in Fig. 2 in the limit $\tilde{\sigma}\rightarrow 0$, which
displays curves of constant $\Omega_{\chi}h^{2}=0.12$ for fixed $m_{\chi}$.

In addition, after inflation, there is the inevitable production of scalars
through their minimal coupling to gravity MO ; cmov . This source of
production can be understood as a single graviton exchange coupling the
inflaton condensate to the scalar. In the absence of scalar self-interactions
or direct coupling to the inflaton the energy density resulting from the
large-scale fluctuations dominate the total.

In this work, we have considered the possibility of a Planck-suppressed
coupling between the scalar and the inflaton. Indeed, such couplings are
generated by graviton exchange, but these are derivative couplings that do not
affect the effective mass of $\chi$. Nevertheless, without a symmetry
preventing such a coupling, we would expect Planck-suppressed coupling,
$\sigma\phi^{2}\chi^{2}$, to be present. Such a coupling will also be
generated if $\chi$ is non-minimally coupled to curvature Clery:2022wib . We
have normalized the coupling to the inflaton mass so that
$\sigma={\tilde{\sigma}}m_{\phi}^{2}/M_{P}^{2}$.

For extremely small $\tilde{\sigma}$, we have seen that the relic density is
indeed determined almost solely by the large-scale fluctuations and therefore
there are very strong constraints on the mass and reheating temperature due to
the overproduction of $\chi$ during the phase of inflation. These are
summarized in Eq. (39). However, even for very small coupling,
$\tilde{\sigma}\gtrsim 10^{-7}$, the allowed reheating temperature range is
increased as sensitivity to the scalar mass is decreased, as can be seen in
Fig. 2. This is due to an earlier phase of oscillations initiated by the
$\sigma\phi^{2}$ term and enhanced dilution.

For larger ${\tilde{\sigma}}\geq\frac{9}{16}$, direct particle production from
the condensate becomes important and soon dominates. In this regime, there is
a wide range of scalar masses (from 50 MeV to $10^{13}$ GeV) and a wide range
in reheating temperatures (from $\sim 1$ to $\sim 10^{15}$ GeV) for
$0.004<\tilde{\sigma}<50$, where the lower bound comes from the constraint
obtained from isocurvature perturbations and the upper bound from ensuring
that non-perturbative processes such as parametric resonance is unimportant.
These are our main results and are displayed in the right panel of Fig. 2 and
Fig. 4.

We have provided a combination of analytical and numerical results to
calculate the relic density. A more complete numerical analysis will be
presented elsewhere paper2 . These results can naturally be generalized by
including scalar self-interactions $\lambda_{p}\chi^{p}$, as well as
considering background cosmologies in which the inflaton potential is not
quadratic about its minimum. The case for $p=4$ was considered in
Markkanen:2018gcw ; Cosme:2020nac .

###### Acknowledgements.

This project has received support from the European Union’s Horizon 2020
research and innovation program under the Marie Sklodowska-Curie grant
agreement No 860881-HIDDeN and the CNRS-IRP project UCMN. The authors would
like to acknowledge Simon Clery, Mathieu Gross, and Jong-Hyun Yoon for useful
discussions. The work of G.C. and K.A.O. was supported in part by DOE grant
DE-SC0011842 at the University of Minnesota. The work of M.A.G.G. was
supported by the DGAPA-PAPIIT grant IA103123 at UNAM, the CONAHCYT “Ciencia de
Frontera” grant CF-2023-I-17, and the Programa de Investigación del Instituto
de Física 2023 (PIIF23). S.V. was supported in part by DOE grant DE-SC0022148
at the University of Florida.

## References

  * (1) T.S. Bunch and P.C.W. Davies, Proc. R. Soc. London A360 (1978) 117; A. Vilenkin, Phys. Lett. 115B, 91 (1982); A. Vilenkin and L. H. Ford, Phys. Rev. D 26, 1231 (1982); A. D. Linde, Phys. Lett. 116B, 335 (1982); A. A. Starobinsky, Phys. Lett. 117B, 175 (1982); A. Vilenkin, Nucl. Phys. B 226, 527 (1983); A. D. Linde, Phys. Lett. 160B, 243 (1985); K. Enqvist, K. W. Ng and K. A. Olive, Nucl. Phys. B 303, 713 (1988).
  * (2) A. D. Linde, Phys. Lett. B 131 (1983), 330-334
  * (3) T. Gherghetta, C. F. Kolda and S. P. Martin, Nucl. Phys. B 468 (1996), 37-58 [arXiv:hep-ph/9510370 [hep-ph]].
  * (4) I. Affleck and M. Dine, Nucl. Phys. B 249 (1985), 361-380
  * (5) A. D. Linde, Phys. Lett. B 160 (1985), 243-248
  * (6) B. A. Campbell, J. R. Ellis, D. V. Nanopoulos and K. A. Olive, Mod. Phys. Lett. A 1 (1986), 389; J. R. Ellis, D. V. Nanopoulos and K. A. Olive, Phys. Lett. B 184 (1987), 37-44; J. R. Ellis, K. Enqvist, D. V. Nanopoulos and K. A. Olive, Phys. Lett. B 191 (1987), 343-348
  * (7) K. Enqvist and A. Mazumdar, Phys. Rept. 380 (2003), 99-234 [arXiv:hep-ph/0209244 [hep-ph]].
  * (8) M. A. G. Garcia and K. A. Olive, JCAP 09 (2013), 007 [arXiv:1306.6119 [hep-ph]].
  * (9) F. Graziani and K. A. Olive, Phys. Lett. B 216 (1989), 31-36
  * (10) K. Enqvist, D. G. Figueroa and G. Rigopoulos, JCAP 01 (2012), 053 [arXiv:1109.3024 [astro-ph.CO]].
  * (11) R. Allahverdi and A. Mazumdar, JCAP 10 (2006), 008 [arXiv:hep-ph/0512227 [hep-ph]]; R. Allahverdi, K. Enqvist, J. Garcia-Bellido, A. Jokinen and A. Mazumdar, JCAP 06 (2007), 019 [arXiv:hep-ph/0610134 [hep-ph]]; R. Allahverdi and A. Mazumdar, Phys. Rev. D 78 (2008), 043511 [arXiv:0802.4430 [hep-ph]].
  * (12) K. A. Olive and M. Peloso, Phys. Rev. D 74 (2006), 103514 [arXiv:hep-ph/0608096 [hep-ph]]; A. E. Gumrukcuoglu, K. A. Olive, M. Peloso and M. Sexton, Phys. Rev. D 78 (2008), 063512 [arXiv:0805.0273 [hep-ph]].
  * (13) M. Dine, L. Randall and S. D. Thomas, Phys. Rev. Lett. 75 (1995), 398-401 [arXiv:hep-ph/9503303 [hep-ph]]; M. Dine, L. Randall and S. D. Thomas, Nucl. Phys. B 458 (1996), 291-326 [arXiv:hep-ph/9507453 [hep-ph]].
  * (14) P. Binetruy and M. K. Gaillard, Phys. Lett. B 195 (1987), 382-388
  * (15) E. Cremmer, S. Ferrara, C. Kounnas and D. V. Nanopoulos, Phys. Lett. B 133 (1983) 61; A. B. Lahanas and D. V. Nanopoulos, Phys. Rept. 145 (1987) 1.
  * (16) M. K. Gaillard, H. Murayama and K. A. Olive, Phys. Lett. B 355 (1995) 71 [hep-ph/9504307]; B. A. Campbell, M. K. Gaillard, H. Murayama and K. A. Olive, Nucl. Phys. B 538 (1999), 351-367 [arXiv:hep-ph/9805300 [hep-ph]].
  * (17) M. S. Turner and L. M. Widrow, Phys. Rev. D 37 (1988), 3428
  * (18) P. J. E. Peebles and A. Vilenkin, Phys. Rev. D 60 (1999), 103506 [arXiv:astro-ph/9904396 [astro-ph]].
  * (19) K. Enqvist, S. Nurmi, T. Tenkanen and K. Tuominen, JCAP 08 (2014), 035 [arXiv:1407.0659 [astro-ph.CO]].
  * (20) S. Nurmi, T. Tenkanen and K. Tuominen, JCAP 11, 001 (2015) [arXiv:1506.04048 [astro-ph.CO]].
  * (21) O. Bertolami, C. Cosme and J. G. Rosa, Phys. Lett. B 759 (2016), 1-8 [arXiv:1603.06242 [hep-ph]].
  * (22) G. Alonso-Álvarez and J. Jaeckel, JCAP 10 (2018), 022 [arXiv:1807.09785 [hep-ph]].
  * (23) T. Markkanen, A. Rajantie and T. Tenkanen, Phys. Rev. D 98 (2018) no.12, 123532 [arXiv:1811.02586 [astro-ph.CO]].
  * (24) T. Tenkanen, Phys. Rev. Lett. 123 (2019) no.6, 061302 [arXiv:1905.01214 [astro-ph.CO]].
  * (25) G. Choi, H. J. He and E. D. Schiappacasse, JCAP 10, 043 (2019) [arXiv:1906.02094 [astro-ph.CO]].
  * (26) C. Cosme and T. Tenkanen, Phys. Rev. D 102 (2020) no.12, 123534 [arXiv:2009.01149 [astro-ph.CO]].
  * (27) O. Lebedev, JCAP 02, 032 (2023) [arXiv:2210.02293 [hep-ph]].
  * (28) Y. Ema, R. Jinno, K. Mukaida and K. Nakayama, JCAP 05, 038 (2015) [arXiv:1502.02475 [hep-ph]]; Y. Ema, R. Jinno, K. Mukaida and K. Nakayama, Phys. Rev. D 94, no.6, 063517 (2016) [arXiv:1604.08898 [hep-ph]]; Y. Ema, K. Nakayama and Y. Tang, JHEP 09, 135 (2018) [arXiv:1804.07471 [hep-ph]].
  * (29) K. Dimopoulos, Phys. Rev. D 74, 083502 (2006) [arXiv:hep-ph/0607229 [hep-ph]].
  * (30) P. W. Graham, J. Mardon and S. Rajendran, Phys. Rev. D 93, no.10, 103520 (2016) [arXiv:1504.02102 [hep-ph]].
  * (31) M. Garny, M. Sandora and M. S. Sloth, Phys. Rev. Lett. 116, no.10, 101302 (2016) [arXiv:1511.03278 [hep-ph]]; M. Garny, A. Palessandro, M. Sandora and M. S. Sloth, JCAP 02, 027 (2018) [arXiv:1709.09688 [hep-ph]].
  * (32) Y. Tang and Y. L. Wu, Phys. Lett. B 774, 676-681 (2017) [arXiv:1708.05138 [hep-ph]].
  * (33) N. Bernal, M. Dutra, Y. Mambrini, K. Olive, M. Peloso and M. Pierre, Phys. Rev. D 97 (2018) no.11, 115020 [arXiv:1803.01866 [hep-ph]].
  * (34) Y. Ema, K. Nakayama and Y. Tang, JHEP 07, 060 (2019) [arXiv:1903.10973 [hep-ph]].
  * (35) M. Chianese, B. Fu and S. F. King, JCAP 06, 019 (2020) [arXiv:2003.07366 [hep-ph]]; M. Chianese, B. Fu and S. F. King, JCAP 01, 034 (2021) [arXiv:2009.01847 [hep-ph]].
  * (36) A. Ahmed, B. Grzadkowski and A. Socha, JHEP 08, 059 (2020) [arXiv:2005.01766 [hep-ph]].
  * (37) E. W. Kolb and A. J. Long, JHEP 03, 283 (2021) [arXiv:2009.03828 [astro-ph.CO]].
  * (38) M. Redi, A. Tesi and H. Tillim, JHEP 05, 010 (2021) [arXiv:2011.10565 [hep-ph]].
  * (39) S. Ling and A. J. Long, Phys. Rev. D 103, no.10, 103532 (2021) [arXiv:2101.11621 [astro-ph.CO]].
  * (40) Y. Mambrini and K. A. Olive, Phys. Rev. D 103 (2021) no.11, 115009 [arXiv:2102.06214 [hep-ph]].
  * (41) N. Bernal and C. S. Fong, JCAP 06 (2021), 028 [arXiv:2103.06896 [hep-ph]]; X. Sun, [arXiv:2112.04304 [hep-ph]].
  * (42) B. Barman and N. Bernal, JCAP 06, 011 (2021) [arXiv:2104.10699 [hep-ph]].
  * (43) M. R. Haque and D. Maity, Phys. Rev. D 106, no.2, 023506 (2022) [arXiv:2112.14668 [hep-ph]].
  * (44) S. Clery, Y. Mambrini, K. A. Olive and S. Verner, Phys. Rev. D 105, no.7, 075005 (2022) [arXiv:2112.15214 [hep-ph]].
  * (45) S. Clery, Y. Mambrini, K. A. Olive, A. Shkerin and S. Verner, Phys. Rev. D 105 (2022) no.9, 095042 [arXiv:2203.02004 [hep-ph]].
  * (46) Y. Mambrini, K. A. Olive and J. Zheng, JCAP 10 (2022), 055 [arXiv:2208.05859 [hep-ph]].
  * (47) E. Basso, D. J. H. Chung, E. W. Kolb and A. J. Long, JHEP 12, 108 (2022) [arXiv:2209.01713 [gr-qc]].
  * (48) M. R. Haque, D. Maity and R. Mondal, JHEP 09, 012 (2023) [arXiv:2301.01641 [hep-ph]].
  * (49) M. A. G. Garcia, M. Pierre and S. Verner, Phys. Rev. D 107 (2023) no.4, 043530 [arXiv:2206.08940 [hep-ph]].
  * (50) K. Kaneta, S. M. Lee and K. y. Oda, JCAP 09, 018 (2022) [arXiv:2206.10929 [astro-ph.CO]].
  * (51) E. W. Kolb, S. Ling, A. J. Long and R. A. Rosen, JHEP 05, 181 (2023) [arXiv:2302.04390 [astro-ph.CO]].
  * (52) K. Kaneta and K. y. Oda, [arXiv:2304.12578 [hep-ph]].
  * (53) M. A. G. Garcia, M. Pierre and S. Verner, [arXiv:2305.14446 [hep-ph]].
  * (54) K. Kaneta, W. Ke, Y. Mambrini, K. A. Olive and S. Verner, Phys. Rev. D 108 (2023) no.11, 115027 [arXiv:2309.15146 [hep-ph]].
  * (55) A. A. Starobinsky and J. Yokoyama, Phys. Rev. D 50, 6357-6368 (1994) [arXiv:astro-ph/9407016 [astro-ph]].
  * (56) A. A. Starobinsky, Phys. Lett. B 91, 99 (1980).
  * (57) M. A. G. Garcia, K. Kaneta, Y. Mambrini and K. A. Olive, Phys. Rev. D 101 (2020) no.12, 123507 [arXiv:2004.08404 [hep-ph]].
  * (58) M. A. G. Garcia, K. Kaneta, Y. Mambrini and K. A. Olive, JCAP 04 (2021), 012 [arXiv:2012.10756 [hep-ph]].
  * (59) Y. Akrami et al. [Planck], Astron. Astrophys. 641 (2020), A10 doi:10.1051/0004-6361/201833887 [arXiv:1807.06211 [astro-ph.CO]].
  * (60) J. Ellis, D. V. Nanopoulos and K. A. Olive, Phys. Rev. Lett. 111 (2013) 111301 [arXiv:1305.1247 [hep-th]].
  * (61) J. Ellis, M. A. G. Garcia, D. V. Nanopoulos and K. A. Olive, JCAP 07, 050 (2015) [arXiv:1505.06986 [hep-ph]].
  * (62) J. Ellis, M. A. G. Garcia, D. V. Nanopoulos, K. A. Olive and S. Verner, Phys. Rev. D 105, no.4, 043504 (2022) [arXiv:2112.04466 [hep-ph]].
  * (63) Y. Ema, M. A. G. Garcia, W. Ke, K. A. Olive and S. Verner, [arXiv:2404.14545 [hep-ph]].
  * (64) M. A. G. Garcia, K. Kaneta, Y. Mambrini, K. A. Olive and S. Verner, JCAP 03, no.03, 016 (2022) [arXiv:2109.13280 [hep-ph]].
  * (65) D. J. H. Chung, Phys. Rev. D 67, 083514 (2003) [arXiv:hep-ph/9809489 [hep-ph]].
  * (66) T. Markkanen and A. Rajantie, JHEP 01, 133 (2017) [arXiv:1607.00334 [gr-qc]].
  * (67) D. J. H. Chung, E. W. Kolb and A. J. Long, JHEP 01, 189 (2019) [arXiv:1812.00211 [hep-ph]].
  * (68) D. Racco, S. Verner and W. Xue, [arXiv:2405.13883 [hep-ph]].
  * (69) G. Choi, M. A. G. Garcia, W. Ke, Y. Mambrini, K. A. Olive, and S. Verner (in preparation).
  * (70) B. R. Holstein, Am. J. Phys. 74, 1002-1011 (2006) [arXiv:gr-qc/0607045 [gr-qc]].
  * (71) G. F. Giudice, R. Rattazzi and J. D. Wells, Nucl. Phys. B 544, 3-38 (1999) [arXiv:hep-ph/9811291 [hep-ph]].
  * (72) D. J. H. Chung, E. W. Kolb, A. Riotto and L. Senatore, Phys. Rev. D 72, 023511 (2005) [arXiv:astro-ph/0411468 [astro-ph]].
  * (73) M. A. G. Garcia, M. Pierre and S. Verner, Phys. Rev. D 107, no.12, 123508 (2023) [arXiv:2303.07359 [hep-ph]].
  * (74) M. Redi and A. Tesi, JHEP 10, 167 (2022) [arXiv:2204.14274 [hep-ph]].
  * (75) Y. Mambrini, Particles in the dark Universe, Springer Ed., ISBN 978-3-030-78139-2 (2021).
  * (76) M. Kawasaki, K. Kohri and N. Sugiyama, Phys. Rev. D 62, 023506 (2000) [arXiv:astro-ph/0002127 [astro-ph]]. P. F. de Salas, M. Lattanzi, G. Mangano, G. Miele, S. Pastor and O. Pisanti, Phys. Rev. D 92, no.12, 123534 (2015) [arXiv:1511.00672 [astro-ph.CO]]. S. Hannestad, Phys. Rev. D 70, 043506 (2004) [arXiv:astro-ph/0403291 [astro-ph]]; T. Hasegawa, N. Hiroshima, K. Kohri, R. S. L. Hansen, T. Tram and S. Hannestad, JCAP 12, 012 (2019) [arXiv:1908.10189 [hep-ph]].
  * (77) G. Choi, W. Ke and K. A. Olive, [arXiv:2402.04310 [hep-ph]].
  * (78) A. A. Starobinsky, Lect. Notes Phys. 246, 107-126 (1986) doi:10.1007/3-540-16452-9_6
  * (79) A. Y. Kamenshchik, A. A. Starobinsky and T. Vardanyan, Eur. Phys. J. C 82, no.4, 345 (2022) doi:10.1140/epjc/s10052-022-10295-z [arXiv:2109.05625 [gr-qc]].
  * (80) P. R. Anderson and E. Mottola, Phys. Rev. D 89, 104038 (2014) [arXiv:1310.0030 [gr-qc]].
  * (81) P. R. Anderson and E. Mottola, Phys. Rev. D 89, 104039 (2014) [arXiv:1310.1963 [gr-qc]].
  * (82) L. Li, T. Nakama, C. M. Sou, Y. Wang and S. Zhou, JHEP 07, 067 (2019) [arXiv:1903.08842 [astro-ph.CO]].
  * (83) S. P. Corbà and L. Sorbo, JCAP 07, 005 (2023) [arXiv:2209.14362 [hep-th]].
  * (84) T. Kobayashi and N. Afshordi, JHEP 10, 166 (2014) [arXiv:1408.4141 [hep-th]].
  * (85) T. Hayashinaka, T. Fujita and J. Yokoyama, JCAP 07, 010 (2016) [arXiv:1603.04165 [hep-th]].
  * (86) M. Banyeres, G. Domènech and J. Garriga, JCAP 10, 023 (2018) [arXiv:1809.08977 [hep-th]].
  * (87) G. W. Gibbons and S. W. Hawking, Phys. Rev. D 15, 2738-2751 (1977)

##  Appendix A Fokker-Planck Approach

To compute the expectation value of $\chi$, we use the stochastic approach
that treats the long wavelength superhorizon contribution of the field
$\chi(t,\mathbf{x})$ in a de Sitter background as a classical stochastic
variable, $\varphi$, characterized by the probability distribution
$\rho(t,\varphi)$ Starobinsky:1986fx ; Starobinsky:1994bd ;
Kamenshchik:2021tjh . This probability distribution must satisfy the Fokker-
Planck equation

$\frac{\partial\rho}{\partial
t}\;=\;\frac{H_{I}^{3}}{8\pi^{2}}\frac{\partial^{2}\rho}{\partial\varphi^{2}}+\frac{1}{3H_{I}}\frac{\partial}{\partial\phi}\left(\frac{\partial
V}{\partial\varphi}\rho(t,\varphi)\right)\,.$ (52)

This expression implies that the long wavelength contribution of an effective
potential $V(s(t,\mathbf{x}))$ is equal to the expectation value of the
effective potential of the classical stochastic variable $V(\varphi)$.
Therefore, using the effective mass, $m_{\chi,\rm eff}$, the spectator field
potential in terms of the classical stochastic variable can be expressed as

$V(\varphi)\;=\;\frac{1}{2}m_{\chi,\rm eff}^{2}\varphi^{2}\,.$ (53)

At late times, the solution to the Fokker-Planck equation reaches static
equilibrium and becomes

$\rho_{\mathrm{eq}}(\varphi)=N^{-1}\exp\left(-\frac{8\pi^{2}}{3H_{I}^{4}}V(\varphi)\right)\,,$
(54)

where $N$ is fixed according to the normalization condition

$\int_{-\infty}^{\infty}\rho_{\mathrm{eq}}(\varphi)d\varphi=1\,.$ (55)

For the effective potential (53), we find that the normalization constant $N$
can be expressed as

$N\;=\;\int_{-\infty}^{\infty}\exp\left[-\frac{8\pi^{2}}{3H_{I}^{4}}\left(\frac{m_{\chi,\rm
eff}^{2}\varphi^{2}}{2}\right)\right]d\varphi\;=\;\sqrt{\frac{3H_{I}^{4}}{4\pi
m_{\chi,\rm eff}^{2}}}\,.$ (56)

Next, we compute the expectation value of $\langle\chi^{2}\rangle$, given by

$\displaystyle\langle\chi^{2}\rangle\;=\;\langle\varphi^{2}\rangle\;=\;\int_{-\infty}^{\infty}\varphi^{2}\rho_{\rm
eq}(\varphi)d\varphi\;=\;\frac{1}{N}\frac{3\sqrt{3}H_{I}^{6}}{16\pi^{\frac{5}{2}}m_{\chi,\rm
eff}^{3}}\;=\;\frac{3H_{I}^{4}}{8\pi^{2}m_{\chi,\rm eff}^{2}}\,.$ (57)

We note that when considering a pure de Sitter background, the asymptotic
value is determined by the Hubble parameter, $H_{I}$. However, numerical
studies show that for a slow-roll model of inflation, the density of $\chi$ is
a better fit using $H_{\rm end}$ paper2 . Thus, at the end of inflation, we
use the following expectation value for the spectator field:

$\langle\chi^{2}\rangle_{\rm end}\;\simeq\frac{3H_{\rm
end}^{4}}{8\pi^{2}m^{2}_{\chi,\rm eff}(a_{\rm end})}\,.$ (58)