[description= Comoving momentum: ]commov$q\,T_{\star}$ [description= Metric
perturbations in synchronous gauge: ]handeta$h$ and $\eta$ [description= NCDM
“temperature”: ]ftemp$T_{\star}$ [description= Phase space distribution:
]fased$f$ [description= Equation of state: ]eqstate$w$ [description= Phase
space distribution perturbation: ]psif$\Psi$ [description= Transfer function:
]trans$\mathcal{T}$ [description= WDM mass: ]wdmass$m_{\rm WDM}$ [description=
NCDM mass: ]ncdmm$m_{\rm DM}$ [description=Momentum (modulus)]p$p$
[description=Momentum]4$\mathbf{p}$ [description=Energy (of a
particle)]5$p_{0}$ [description= Relative density fluctuation: ]deltar$\delta$
[description= Pressure: ]Pp$P$ [description= Energy density: ]rhop$\rho$
[description= Velocity divergence: ]thetap$\theta$ [description= Anisotropic
stress: ]sigmap$\sigma$ [description= Scale factor (FLRW): ]ap$a$
[description= Hubble parameter (function): ]Hp$H$ [description= Conformal
Hubble parameter (function): ]mathcalH$\mathcal{H}$ [description= Comoving
“energy”: ]comene$\epsilon$ [description=Conformal time’]13$\tau$
[description=Cosmic time]14$t$ [description= Phase space distribution
parametrization: ]fasedp$\alpha,\beta,\gamma$ [description= Reduced Planck
mass: ]mplanck$M_{P}$ [description= Collision term:
]collisionterm$\mathcal{C}$ [description= Number density of the species
$\chi$: ]nchi$n_{\chi}$ [description= Normalization scale of the $2\rightarrow
2$ cross section: ]Lambda$\Lambda$ [description= Rest-frame sound speed:
]cs$\hat{c}_{s}$ [description= Adiabatic sound speed: ]ca$c_{a}$ [description=
Free-streaming scale: ]kFS$k_{\text{FS}}$ [description= Free-streaming
horizon: ]kH$k_{\text{H}}$ [description= Reheating temperature: ]Treh$T_{\rm
reh}$ [description= Highest temperature during reheating: ]Tmax$T_{\rm max}$
[description= Internal degrees of freedom of the species $\chi$:
]gchi$g_{\chi}$ [description= Reheating temperature parametrization: ]b$b$
[description= Effective number of relativistic species: ]Neff$N_{\rm eff}$
[description= Number of effective entropy degrees of freedom:
]gstars$g_{\rm*s}$ [description= Inflaton scalar field: ]inflaton$\phi$
[description= Inflaton decay width: ]Gammaphi$\Gamma_{\phi}$

IFT-UAM/CSIC-20-135

How warm are non-thermal relics?
Lyman-$\alpha$ bounds on out-of-equilibrium dark matter

Guillermo Ballesteros, Marcos A. G. Garcia, Mathias Pierre

Instituto de Física Teórica UAM/CSIC,
Calle Nicolás Cabrera 13-15, Cantoblanco E-28049 Madrid, Spain
Departamento de Física Teórica, Universidad Autónoma de Madrid (UAM)
Campus de Cantoblanco, 28049 Madrid, Spain

Abstract

> We investigate the power spectrum of Non-Cold Dark Matter (NCDM) produced in
> a state out of thermal equilibrium. We consider dark matter production from
> the decay of scalar condensates (inflaton, moduli), the decay of thermalized
> and non-thermalized particles, and from thermal and non-thermal freeze-in.
> For each case, we compute the NCDM phase space distribution and the linear
> matter power spectrum, which features a cutoff analogous to that for Warm
> Dark Matter (WDM). This scale is solely determined by the equation of state
> of NCDM. We propose a mapping procedure that translates the WDM
> Lyman-$\alpha$ mass bound to NCDM scenarios. This procedure does not require
> expensive ad hoc numerical computations of the non-linear matter power
> spectrum. By applying it, we obtain bounds on several NCDM possibilities,
> ranging from $m_{\rm DM}\gtrsim{\rm EeV}$ for DM production from inflaton
> decay with a low reheating temperature, to sub-keV values for non-thermal
> freeze-in. We discuss the phenomenological implications of these results for
> specific examples which include strongly-stabilized and non-stabilized
> supersymmetric moduli, gravitino production from inflaton decay,
> $Z^{\prime}$ and spin-2 mediated freeze-in, and non-supersymmetric spin-3/2
> DM.

###### Contents

  1. 1 Introduction and results
    1. 1.1 Motivation
    2. 1.2 Ly-$\alpha$ constraints on out-of-equilibrium dark matter
    3. 1.3 Dark matter production mechanisms
  2. 2 Non-cold dark matter cosmology
    1. 2.1 Linear cosmological perturbation theory
    2. 2.2 Large scale structure
    3. 2.3 Analytical rescaling and generalized phase space distribution
  3. 3 Decay of a classical condensate
    1. 3.1 Perturbative inflaton decay
      1. 3.1.1 DM phase space distribution
      2. 3.1.2 Power spectrum and Ly-$\alpha$ constraints
      3. 3.1.3 Relic density and phenomenology
    2. 3.2 Moduli decays
      1. 3.2.1 DM phase space distribution
      2. 3.2.2 Power spectrum and Ly-$\alpha$ constraints
      3. 3.2.3 Relic density and phenomenology
  4. 4 Freeze-in via decay
    1. 4.1 Thermal decay
      1. 4.1.1 DM phase space distribution
      2. 4.1.2 Power spectrum and Ly-$\alpha$ constraints
      3. 4.1.3 Relic density and phenomenology
    2. 4.2 Non-thermal decay
      1. 4.2.1 DM phase space distribution
      2. 4.2.2 Power spectrum and Ly-$\alpha$ constraints
      3. 4.2.3 Relic density and phenomenology
  5. 5 Ultraviolet freeze-in via scatterings
    1. 5.1 Thermal freeze-in
      1. 5.1.1 DM phase space distribution
      2. 5.1.2 Power spectrum and Ly-$\alpha$ constraints
      3. 5.1.3 Relic density and phenomenology
    2. 5.2 Non-thermal freeze-in
      1. 5.2.1 DM phase space distribution
      2. 5.2.2 Power spectrum and Ly-$\alpha$ constraints
      3. 5.2.3 Relic density and phenomenology
  6. 6 Light, but not too light, dark matter
  7. 7 Conclusions
  8. A The Boltzmann equation in an expanding universe
    1. A.1 Generalities
    2. A.2 Freeze-in via scatterings
    3. A.3 Non-thermal freeze-in

##  1 Introduction and results

###  1.1 Motivation

After a few decades of remarkable improvement, dark matter (DM) direct
detection experiments have reached a sensitivity on the nucleon-DM scattering
cross section around $10^{-46}\leavevmode\nobreak\ \text{cm}^{2}$ for DM
masses of the order of the electroweak scale [1]. The absence of any confirmed
experimental signal (also in indirect detection and colliders) strongly
constrains the viable parameter space of Weakly Interacting Massive Particle
(WIMP) models of DM based on the vanilla freeze-out mechanism. This calls for
a reassessment of the attractiveness of this framework in the simplest models
[2, 3, 4]. In this context, exploring theoretically and experimentally other
scenarios [5] that can achieve the correct DM abundance is necessary. A well-
known example of such a scenario is the freeze-in mechanism [6]. Other
examples are a dark sector that thermalizes only with itself [7, 8] and a DM
depletion process lead by cannibalization [9, 10, 11]. These proposals assume
feeble couplings between the SM and the DM,111A large energy scale, even
larger than the reheating temperature, can be invoked to justify suppressed
SM-DM interactions. Just to give some examples, this scale can be identified
with the Planck mass in gravitino DM [12, 13, 14, 15, 16, 17, 18], with the
mass of heavy gauge fields in Grand Unified Theories [19, 20, 21] and with a
new physics threshold in scenarios inspired by modified gravity [22, 23, 24,
25]. helping them to satisfy current bouds. Consequently, this tends to reduce
the chances of testing them by traditional means [5]. However, several
phenomenological studies [7, 26, 27, 28, 29, 30, 31] have highlighted various
possibilities for observing such DM candidates.

Scenarios in which the DM is produced by non-standard mechanisms may feature
an important DM self-interaction cross section or a free-streaming scale,
affecting the large scale structure of the universe. These properties could
allow to alleviate purported tensions in the $\Lambda$CDM model at galactic
and sub-galactic scales [32, 33, 34]. Also, non-thermal DM models have been
proposed to address the tensions between early and late time determinations of
the Hubble constant [35] and of the clustering of matter [36, 37].

Indeed, in absence of thermodynamic equilibrium between the DM and the SM, the
DM phase space distribution can differ significantly from the standard freeze-
out case. This opens a possibility for discriminating between different DM
models and production mechanisms. The DM component in the standard
$\Lambda$CDM model of cosmology is assumed to be entirely pressureless. A non-
vanishing DM kinetic energy would then result in a cutoff in the matter power
spectrum on small wavelength Fourier modes (as compared to $\Lambda$CDM
prediction).

An interesting possibility for testing these Non-Cold Dark Matter (NCDM)
models –which do not conform to the standard freeze-out mechanism– is to
measure the Ly-$\alpha$ forest of absorption lines of light emitted by distant
quasars around redshift $z=2-4$, which is produced due to the neutral hydrogen
present in the intergalactic medium. This provides enough information on the
matter power spectrum at sufficiently small scales for probing the
aforementioned cutoff.

###  1.2 Ly-$\alpha$ constraints on out-of-equilibrium dark matter

The well known Ly-$\alpha$ bound on the DM mass for Warm Dark Matter (WDM)
[38, 39, 33, 40, 41, 42, 43],

$\;\gtrsim\;(1.9-5.3)\leavevmode\nobreak\ \text{keV at 95\% C.L.}\,,$ (1.1)

can be mapped into constraints on various out-of-equilibrium NCDM production
mechanisms. To do this, we compute (for the first time) the phase space
distributions in several of these models by integrating the Boltzmann
transport equation, numerically and/or analytically, depending on the
production process. For the large majority of the scenarios that we consider,
the resulting phase space distributions can be remarkably well described by a
generalized distribution of the form

$f(q)\,\propto\,q^{\alpha}\,\exp{\left(-\beta\,q^{\gamma}\right)}\,,$ (1.2)

where $q$ denotes the DM comoving momentum and are model-dependent constants.
We then use CLASS [44, 45] to compute, for each of the NCDM models we
consider, the linear power spectrum $\mathcal{P}_{\text{NCDM}}(k)$, or, more
precisely, the linear transfer function, defined in terms of the ratio to the
$\Lambda$CDM spectrum as follows,

$(k)\,\equiv\,\left(\dfrac{\mathcal{P}_{\text{NCDM}}(k)}{\mathcal{P}_{\Lambda\text{CDM}}(k)}\right)^{1/2}\,.$
(1.3)

We assume that the DM is entirely composed of a single NCDM species (and is
produced by only one mechanism in each scenario). By varying the DM mass, we
match the transfer function to the one of a fermionic WDM scenario (for which
the bound (1.1) applies). We perform this matching numerically and, also, with
an approximate semi-analytical procedure, demonstrating their equivalence. We
find that the matching can be done with great accuracy for all the models we
consider (and for all the relevant ranges of their parameters). In addition,
by approximating the NCDM species as a perfect fluid, we show that the cutoff
in the transfer function can be entirely characterized in terms of the
equation of state parameter of NCDM, $w$, which allows to translate the WDM
mass Ly-$\alpha$ limit (1.1) to the NCDM case. This can be done for each of
the NCDM models, without having to run specifically tailored N-body
simulations or doing a dedicated non-linear analysis of the NCDM perturbations
at small scales. A general analytical expression relates these bounds to each
other via only the knowledge of the first and second moments of the phase
space distribution,

$\,=\,m_{\text{WDM}}\left(\dfrac{T_{\star}}{T_{\text{WDM},0}}\right)\sqrt{\dfrac{\langle
q^{2}\rangle}{\langle
q^{2}\rangle_{\text{WDM}}}}\quad\quad\text{(Lyman-$\alpha$ bound)}\,,$ (1.4)

where $\propto\langle q\rangle$ is the present “temperature” of NCDM,
understood as the energy scale that normalizes the typical momentum of the
distribution. In all cases, we find a remarkable agreement between this
approximation and the numerical computation of the linear power spectrum.
Using this procedure, we achieve (for most of our scenarios and in the range
of scales of interest) a $\lesssim 3\%$ error in the matching to the transfer
function of WDM, see Figure 3. In the (very few) least precise examples that
we consider the matching worsens to at most $10\%$.

Figure 1: Pipeline applied in this paper to derive bounds on a given NCDM
model from the Lyman-$\alpha$ WDM mass limit. The boxes represent the
computational steps or inputs used to derive a bound on the NCDM mass.
Starting from the NCDM collision term and Lyman-$\alpha$ bound on WDM, two
possible paths allow to derive a bound on NCDM. Our matching procedure, going
through the equation of state matching, allows us to obtain the NCDM bound
without computing numerically the NCDM transfer function.  Figure 2:  Out-of-
equilibrium dark matter as a probe for early universe dynamics. The top figure
depicts schematically the phase space distributions for the production
scenarios considered in the present work, mapped over the history of the early
universe. The relic distributions for the thermal freeze-in production with
$n\geq 6$ and non-thermal freeze-in with $n>2$ are set at the thermalization
time-scale, after inflation but well before the end of reheating. The relic
abundance and the distribution of dark matter produced from inflaton decay or
low-$n$ freeze-in are set at the matter-radiation transition at the end of
reheating. Dark matter can be also produced through decays of thermalized and
non-thermalized fields during early radiation domination. Most distributions
have the approximate form $f(q)\propto q^{\alpha}e^{-\beta q^{\gamma}}$. The
bottom left figure shows that, in all cases, the transfer function of the
linear power spectrum determines the lower bound on the dark matter mass,
since $\mathcal{T}(k)$ can be matched to the corresponding WDM bound. The
bottom right figure shows the range of masses that can be bound by Ly-$\alpha$
observations, depending on the dark matter production mechanism and the
reheating temperature (here $T_{\rm reh}=10^{10}\,{\rm GeV}$).

This bound-mapping procedure is summarized in Figs. 1 and 2. The top panel of
the latter shows schematically the different shapes of the distribution
functions corresponding to six distinct NCDM production processes, active
during or after reheating. We now proceed to enumerate these processes, and
summarize our results for each of them.

###  1.3 Dark matter production mechanisms

We consider various non-equilibrium DM production mechanisms. They are all
assumed to proceed perturbatively via the scattering or decay of particles,
with time-scales ranging from the very end of inflation, during the earliest
stages of reheating, to the radiation dominated universe occurring after the
end of reheating. We list these mechanisms below; see also Fig. 2.

Inflaton decay (Section 3.1). It is often assumed that the DM may have been
produced from the decay of the inflaton field, $\phi$. Even in the absence of
tree-level inflaton-DM couplings, DM-SM interactions can generate a non-
vanishing inflaton $\rightarrow$ DM decay channel at higher order in
perturbation theory [46]. Assuming that this decay proceeds perturbatively
through a two-body process, we find that the DM phase space distribution is of
the form $f(q)\propto q^{-3/2}e^{-0.74q^{2}}$, where the power-like behaviour
at low $q$ arises from redshifting during the matter dominated reheating
epoch, and the Gaussian tail comes from the depletion of the inflaton
condensate at the end of reheating. The resulting Ly-$\alpha$ constraint on
the DM mass is proportional to the ratio of the inflaton mass to the reheating
temperature, being $m_{\rm DM}\gtrsim 3.8\,{\rm MeV}$ for $T_{\rm
reh}=10^{10}\,{\rm GeV}$ and $m_{\phi}\simeq 3\times 10^{13}\,{\rm GeV}$, the
later fixed by the measurement of the amplitude of the curvature power
spectrum [47, 48]. Higher reheating temperatures can reduce the limit down to
the keV range, whereas lower reheating temperatures can increase it well
beyond the TeV range.

Moduli decay (Section 3.2). In many SM extensions, in particular in
supergravity and string constructions, there is a plethora of scalar fields
with very weak couplings to the SM (typically of gravitational strength) and
masses that are typically of the order of the weak scale. These fields are
known as moduli, and can have far reaching cosmological consequences if they
are excited away from their vacuum values in the early Universe. We consider
DM production from moduli decays in two scenarios: when the modulus dominates
the energy of the Universe and decays at late times, and when the modulus is
always subdominant to the inflaton/radiation background due to some
stabilization mechanism. In the first case, the shape of the DM phase space
distribution is identical to that for DM produced from inflaton decay, and the
lower bound on $m_{\rm DM}$ is proportional to the ratio of the modulus mass
($m_{Z}$) to its reheating temperature, with $m_{\rm DM}\gtrsim 13\,{\rm GeV}$
for $m_{Z}=10\,{\rm TeV}$ and $T_{\rm reh}=1\,{\rm MeV}$. For the decay of
stabilized moduli, we find non-thermal DM distributions of the form
$f(q)\propto q^{-3/2}e^{-q^{3/2}}$ or $f(q)\propto q^{-1}e^{-q^{2}}$,
depending on whether the modulus decays during or after reheating,
respectively. In these cases, the limit on $m_{\rm DM}$ depends on the ratio
of the modulus mass to the background temperature evaluated at the moment of
its decay, and on the ratio of the inflaton and modulus decay widths.

Thermal and non-thermal decays (Section 4). DM could have been produced also
from the decay of free particles. In such a case the DM phase space
distribution and its present abundance depend strongly on the initial momentum
distribution of the decaying particles. We consider here two possibilities:
the decay of a thermalized particle species during radiation domination
(Section 4.1.1), and the decay of a particle with a non-equilibrium
distribution, assumed to be produced from the decay of the inflaton (Section
4.2). In both cases we assume that the decaying particle is much lighter than
the inflaton, yet much heavier than DM. For the thermal decay case, we find
that the DM inherits a quasi-thermal distribution, $f(q)\sim q^{-1/2}e^{-q}$,
and the bound on its mass is given by $m_{\rm DM}\gtrsim 7\,{\rm keV}$.

For the non-thermal decay, we find that the shape of the distribution is
highly dependent on the momentum of the parent particle when it decays. If
this initial state decays while it is relativistic, the DM inherits the
Gaussian tail of the parent unstable particle, $f(q)\sim
q^{-5/2}e^{-0.74q^{2}}$. The Ly-$\alpha$ constraint is identical to that for
the direct decay of the inflaton to DM, reduced by a factor of $\sim 0.3$. If
instead the decaying particle is non-relativistic, the DM phase space
distribution is highly non-thermal, skewed towards large momenta, and not
suitable for a fit of the form (1.2). The Ly-$\alpha$ constraint depends on
the mass and width of the decaying particles; more specifically proportional
to the ratio of the mass to the temperature $T_{\rm dec}\propto\Gamma^{1/2}$
at which the decay occurs.

Thermal freeze-in via scatterings (Section 5.1). We consider the possibility
of a DM population generated via the freeze-in mechanism by annihilations of
thermalized SM particles. We assume that the typical DM-SM scattering
amplitude can be parametrized by

$|\mathcal{M}|^{2}\;=\;16\pi\frac{s^{\frac{n}{2}+1}}{\Lambda^{n+2}}\,,$ (1.5)

where $n$ is an integer, $\sqrt{s}$, the square root of the Mandelstam
variable, is the center-of-mass energy (in the high-energy limit), and
$\Lambda$ is some high-energy scale.222Small differences on the dependence on
Maldestam variables in the high energy limit can be absorbed into the value of
$\Lambda$.

For $0\leq n<6$, the DM is produced at the end of the reheating process, at
the reheating temperature $T_{\rm reh}$. For $n=0,2$ the resulting DM momentum
distribution is quasi-thermal with $\beta\sim 1$ and $\gamma=1$. Instead, for
$n=4$ it has a nearly Gaussian tail. For these three scenarios, the matching
of the power spectrum to WDM is excellent and the bound translates to $m_{\rm
DM}\gtrsim 6-9\,{\rm keV}$, with the precise value depending on $n$ and the
quantum statistics of the thermalized scatterers.

When $n\geq 6$, most of the DM is produced on the earliest stages of
reheating, at the maximum temperature $T_{\rm max}$. When this is the case,
the fitting expression (1.2) fails. In particular, for $n=6$ $f(q)$
interpolates between a $q^{3}$ behaviour at $q\ll 1$, and an exponential tail
at $q\gg 1$, through a region where $f(q)\sim q^{-3}$. This relatively
complicated form of the distribution translates to an imperfect match with the
WDM power spectrum, which nevertheless leads to a bound of the form $m_{\rm
DM}^{2}\gtrsim 81\,{\rm keV^{2}}/\ln(T_{\rm max}/T_{\rm reh})$.

Non-thermal freeze-in via scatterings (Section 5.2). The delay between the end
of inflation and the onset of thermal equilibrium in the primordial plasma can
leave an imprint on the DM phase space distribution if the parent scatterers
are produced directly from inflaton decays. Inflaton decay products are
typically very energetic, with momenta of the order of the inflaton mass. Only
after a process of soft radiation emission and energy transfer through
scatterings these decay products reach thermal equilibrium. Thermalization
occurs after the beginning of reheating, but well before it ends. As it turns
out, if $n>2$ in (1.5), most of the DM could have been produced non-thermally
by the very first SM particles present in the universe [49].

As a proof of concept, we consider here annihilations with $n=4$. Notably,
under the freeze-in assumption, the transport equation can be solved in a
closed, albeit complicated, form. We find that the approximation $f(q)\sim
q^{-3/2}e^{-2.5q^{2.6}}$ adequately describes the DM phase space distribution.
The power spectrum matching with WDM can be performed accurately, leading to a
minimum DM mass of the form $m_{\rm DM}\propto m_{\phi}^{23/15}T_{\rm
reh}^{-7/15}$, where $m_{\phi}$ is the inflaton mass, and $T_{\rm reh}$ the
reheating temperature. For $T_{\rm reh}=10^{10}\,{\rm GeV}$, $m_{\rm
DM}\gtrsim 0.4\,{\rm keV}$.

Our paper is organized as follows. In Section 2 we review the treatment of
NCDM relics in cosmological linear perturbation theory, and discuss the
properties of the transfer function (1.3) and the re-scaling into NCDM of the
Lyman-$\alpha$ bounds coming from WDM. In Sections 3 to 5 we study the
production mechanisms we just listed, their Ly-$\alpha$ bounds and the
corresponding phenomenological implications. In Section 6 we discuss the
implications for the effective number of relativistic species. We present our
conclusions in Section 7.

Appendix A contains a brief review of the Boltzmann equation in the early
universe, as well a detailed calculation of the generic form of the collision
term for DM production via freeze-in (Appendix A.2), and the integration of
this collision term for the $n=4$ non-thermal freeze-in scenario (Appendix
A.3). A glossary of the main symbols used in this paper is provided in
Appendix LABEL:symbols. We use a natural system of units in which
$k_{B}=\hbar=c=1$.

##  2 Non-cold dark matter cosmology

###  2.1 Linear cosmological perturbation theory

In the standard $\Lambda$CDM model of cosmology, the DM is assumed to be cold
(CDM), i.e. presureless. Therefore, its equation of state parameter $w$ –
defined by the relation $\bar{P}=\,\bar{\rho}$ where $\bar{\rho}$ and
$\bar{P}$ are its (time-dependent) background energy density and pressure – is
exactly vanishing. However, DM particles produced in the early universe, of
thermal or non-thermal origin, would actually possess some momentum
distribution with a non-vanishing averaged momentum $\langle p\rangle\neq 0$,
which could manifests in a deviation from $w=0$ and, possibly, also through
other moments of it. We will now discuss a set of approximations under which
$w$ can be the sole function encoding the deviations from CDM, both at the
level of the background dynamics and for linear perturbations analyses. In
Sections 3–5 we will study concrete examples of such DM creation processes.

The phase space distribution, $f$, of a general cosmological species is a
function of position, momentum and (conformal) time, $\tau$, that
characterizes its energy-momentum tensor. It is convenient to split it into a
time-dependent homogeneous background part, $\bar{f}(|\boldsymbol{p}|,\tau)$,
plus a fluctuation quantified by a function $\ll 1$, such that
$f(\boldsymbol{x},\boldsymbol{p},\tau)=\bar{f}(|\boldsymbol{p}|,\tau)[1+\Psi(\boldsymbol{x},\boldsymbol{p},\tau)]$;
see e.g. [50]. The background energy density and pressure functions,
$\bar{\rho}$ and $\bar{P}$, of a NCDM relic are then

$\bar{\rho}=4\pi\left(\dfrac{T_{\star}}{a}\right)^{4}\int
q^{2}\epsilon\bar{f}(q)\mathop{}\\!\mathrm{d}q\,,\qquad\bar{P}=\dfrac{4\pi}{3}\left(\dfrac{T_{\star}}{a}\right)^{4}\int
q^{2}\dfrac{q^{2}}{\epsilon}\bar{f}(q)\mathop{}\\!\mathrm{d}q\,,$ (2.1)

where is the scale factor of the Universe and $T_{\star}$ is a convenient
energy scale that characterizes the DM density at the present time.

Following the conventions of [45]333Our $T_{\star}$ is a time-independent
quantity denoted by $T_{\text{NCDM},0}$ in [45]. we define

$q=\dfrac{p\,a}{{T_{\star}}}\quad\text{with}\,\quad=\sqrt{q^{2}+\left(\dfrac{m_{\text{DM}}\,a}{T_{\star}}\right)^{2}}\,,$
(2.2)

where the product is the comoving momentum and $p=|\boldsymbol{p}|$ is the
(absolute value) of the momentum of individual NCDM particles. In Fourier
space, the perturbation $\Psi$ of the NCDM phase space distribution can be
expanded in Legendre polynomials $P_{\ell}$ as follows:

$\Psi(\boldsymbol{k},\hat{\boldsymbol{n}},q,\tau)=\sum_{\ell=0}^{\infty}(-i)^{\ell}(2\ell+1)\Psi_{\ell}(\boldsymbol{k},q,\tau)P_{\ell}(\hat{\boldsymbol{k}}\cdot\hat{\boldsymbol{n}})\,,$
(2.3)

where $k$ is the comoving wavenumber of the perturbations in Fourier space,
$\boldsymbol{k}=k\,\hat{\boldsymbol{k}}$ and
$p=\hat{\boldsymbol{n}}\cdot\vec{\boldsymbol{p}}$. The quantities defining the
perturbed energy-momentum tensor are

$\displaystyle\begin{aligned}
\delta\,=\,&4\pi\left(\dfrac{T_{\star}}{a}\right)^{4}\int
q^{2}\epsilon\bar{f}(q)\Psi_{0}\,\mathop{}\\!\mathrm{d}q,\quad&&\text{energy
density fluctuation}\\\
\delta\,=\,&\dfrac{4\pi}{3}\left(\dfrac{T_{\star}}{a}\right)^{4}\int
q^{2}\dfrac{q^{2}}{\epsilon}\bar{f}(q)\Psi_{0}\,\mathop{}\\!\mathrm{d}q,\quad&&\text{pressure
(density) fluctuation}\\\ (\bar{\rho}+\bar{P})\,=\,&4\pi
k\left(\dfrac{T_{\star}}{a}\right)^{4}\int
q^{3}\bar{f}(q)\Psi_{1}\,\mathop{}\\!\mathrm{d}q,\quad&&\text{velocity
divergence}\\\ (\bar{\rho}+\bar{P})\,=\,&\dfrac{8\pi
k}{3}\left(\dfrac{T_{\star}}{a}\right)^{4}\int
q^{2}\dfrac{q^{2}}{\epsilon}\bar{f}(q)\Psi_{2}\,\mathop{}\\!\mathrm{d}q,\quad&&\text{anisotropic
stress.}\end{aligned}$ (2.4)

For decoupled NCDM, the phase space distribution satisfies the collisionless
Boltzmann equation

$\frac{\partial
f}{\partial\tau}+\frac{\mathop{}\\!\mathrm{d}x^{i}}{\mathop{}\\!\mathrm{d}\tau}\frac{\partial
f}{\partial
x^{i}}+\frac{\mathop{}\\!\mathrm{d}q}{\mathop{}\\!\mathrm{d}\tau}\frac{\partial
f}{\partial
q}+\frac{\mathop{}\\!\mathrm{d}n_{i}}{\mathop{}\\!\mathrm{d}\tau}\frac{\partial
f}{\partial n_{i}}\;=\;0\,,$ (2.5)

with $i=1,2,3$ and $\boldsymbol{n}$ being a unitary 3-vector pointing in the
direction of the momentum, as defined above. In the synchronous gauge, this
equation leads to the following system for the quantities $\Psi_{\ell}$,

$\displaystyle\begin{aligned}
\dot{\Psi}_{0}\,=\,&-\dfrac{qk}{\epsilon}\Psi_{1}+\dfrac{1}{6}\dot{h}\dfrac{\mathop{}\\!\mathrm{d}\ln\bar{f}}{\mathop{}\\!\mathrm{d}\ln
q}\,,\\\
\dot{\Psi}_{1}\,=\,&\dfrac{qk}{3\epsilon}\Big{(}\Psi_{0}-2\Psi_{2}\Big{)}\,,\\\
\dot{\Psi}_{2}\,=\,&\dfrac{qk}{5\epsilon}\Big{(}2\Psi_{1}-3\Psi_{3}\Big{)}-\Big{(}\dfrac{1}{15}\dot{h}+\dfrac{2}{5}\dot{\eta}\Big{)}\dfrac{\mathop{}\\!\mathrm{d}\ln\bar{f}}{\mathop{}\\!\mathrm{d}\ln
q}\,,\\\
\dot{\Psi}_{\ell}\,=\,&\dfrac{qk}{(2\ell+1)\epsilon}\Big{(}\ell\Psi_{\ell-1}-(\ell+1)\Psi_{\ell+1}\Big{)}\,,\quad[\ell\geq
3]\end{aligned}$ (2.6)

where are the trace and traceless part of the metric perturbation [50].

For a non-relativistic species, higher multipoles are typically suppressed by
(positive) powers of $q/\epsilon\sim p/m_{\text{DM}}$, making any
$\Psi_{\ell}$ with $\ell\geq 2$ much smaller than $\Psi_{0}$ and $\Psi_{1}$.
In this case, the Boltzmann hierarchy can be truncated imposing
$\Psi_{\ell}=0$ for $\ell>1$, as discussed in [51], whose analysis shows the
validity of this truncation. As argued also in [52], in this (non-
relativistic) case $\Psi_{0}$ depends only mildly on the variable $q$; and the
integrals in (2.4) are dominated by the low $q\ll\epsilon$ regime,444Notice
that for heavy-tailed distributions, the later approximation is no longer
valid and one cannot apply the analytical arguments presented in this section.
For instance in the case where DM particles could have been produced from
Primordial-Black-Hole evaporation [53, 54], where the distribution function
behaves as $\bar{f}(q)\sim 1/q^{5}$ at large $q$. In that case the integral
appearing in the background-pressure expression is always dominated by
$q\gtrsim m_{\text{DM}}a/T_{\star}$, and hence it cannot be restricted to
$q\ll\epsilon$. so that we can identify $\delta
P/\delta\rho\,\simeq\,\bar{P}/\bar{\rho}=w$.

In this limit, the first two equations of (2.6) can be integrated over $q$,
allowing us to describe the NCDM species with a coupled system of (continuity
and Euler) equations:

$\displaystyle\dot{\delta}\,$
$\displaystyle=\,-(1+w)\Big{(}\theta+\dfrac{\dot{h}}{2}\Big{)}-3\mathcal{H}\left(\hat{c}_{s}^{2}-w\right)\delta+9\mathcal{H}^{2}(1+w)\left(\hat{c}_{s}^{2}-c_{a}^{2}\right)\dfrac{\theta}{k^{2}}\,,$
(2.7) $\displaystyle\dot{\theta}\,$
$\displaystyle=\,-\mathcal{H}\left(1-3\hat{c}_{s}^{2}\right)\theta+\dfrac{\hat{c}_{s}^{2}}{1+w}k^{2}\delta\,,$
(2.8)

where $\,\equiv\,\delta\rho/\bar{\rho}$, and $=a$, where
$H\equiv\dot{a}=\mathop{}\\!\mathrm{d}a/\mathop{}\\!\mathrm{d}\tau$. To first
order in $w\ll 1$, the adiabatic sound speed is
${}^{2}\equiv\dot{\bar{P}}/\dot{\bar{\rho}}\simeq 5w/3$. In addition, as shown
in [45], for sufficiently non-relativistic species, the (rest frame) sound
speed555See e.g. [55] for its definition. can be reasonably well approximated
by the adiabatic sounds speed ${}^{2}\simeq c_{a}^{2}$.

Notice that by taking $w=0$ one recovers the usual CDM perturbation equation
$\dot{\delta}=-1/2\dot{h}$. In the NCDM domination era, from the perturbed
Einstein equations, the trace of the metric fluctuation $h$ satisfies the
equation

$\ddot{h}+\mathcal{H}\dot{h}+3(1+3w)\mathcal{H}^{2}\delta\,=\,0\,,$ (2.9)

allowing the system (2.7)–(2.8) to be reduced to

$\ddot{\delta}+\mathcal{H}\dot{\delta}-\frac{3}{2}\mathcal{H}^{2}\left(1-w\dfrac{10}{9}\frac{k^{2}}{\mathcal{H}^{2}}\right)\delta\,=\,0\,.$
(2.10)

In the limit where $w=0$ exactly, overdensities grow “democratically”, i.e.
independently of $k$ (as in $\Lambda$CDM). However, for non-vanishing $w$, at
a given time, there is a suppressed growth for modes larger than the free-
streaming wavenumber $k>k_{\text{FS}}(a)$ with

${}^{2}(a)\,=\,\dfrac{9}{10}\dfrac{\mathcal{H}^{2}}{w}\,=\,\dfrac{3}{2}\dfrac{\mathcal{H}^{2}}{c_{g}^{2}}\,.$
(2.11)

Thus, a cutoff in the power spectrum can be observed at a given time for modes
larger than the free-streaming horizon wavenumber $k_{H}(a)$, which can be
expressed in term of $k_{\text{FS}}$ as [36]

$(a)\,\equiv\,\left[\int_{0}^{a}\dfrac{1}{k_{\text{FS}}(\tilde{a})}\dfrac{\mathop{}\\!\mathrm{d}\tilde{a}}{\tilde{a}}\right]^{-1}\,.$
(2.12)

From these equations, we see that $w$ is the only quantity (together with the
current DM density) that controls in first approximation the suppression of
the power spectrum at large $k$. In the non-relativistic limit, $w$ can be
expressed in terms of the normalized second moment of the distribution
function

$w\,\simeq\dfrac{\delta
P}{\delta\rho}\,=\,\dfrac{T_{\star}^{2}}{3m_{\text{DM}}^{2}}\frac{\langle
q^{2}\rangle}{a^{2}}\,,$ (2.13)

with the $n$-th moment being

$\langle q^{n}\rangle\,\equiv\,\dfrac{\int
q^{n+2}\bar{f}(q)\mathop{}\\!\mathrm{d}q}{\int
q^{2}\bar{f}(q)\mathop{}\\!\mathrm{d}q}\,.$ (2.14)

As a result, given a phase space distribution for the DM, determination of its
second moment is sufficient to estimate the cutoff of the matter power
spectrum.

###  2.2 Large scale structure

For a given NCDM cosmology, the cutoff can be described in terms of the
transfer function $\mathcal{T}(k)$ defined as

$\mathcal{T}(k)=\left(\dfrac{\mathcal{P}(k)}{\mathcal{P}_{\Lambda\text{CDM}}(k)}\right)^{1/2}\,,$
(2.15)

which compares (at the present time) the power spectrum for a given NCDM
cosmology to the typical $\Lambda$CDM case. As we will now discuss, a small
scale cutoff in the matter power spectrum may be one of the few possibilities
at our disposal for distinguishing NCDM cosmologies from the paradigmatic
$\Lambda$CDM model and thus probe the degree of DM “warmness”. Light emitted
by distant quasars and subsequently interacting with the neutral Hydrogen of
the intergalactic medium around redshifts $z\sim 2-6$ generates a pattern of
absorption lines around $\sim 1000\,$Å: the Ly-$\alpha$ forest. This allows to
probe the power spectrum on scales $k\sim(0.1-10)\,h\,\text{Mpc}^{-1}$ at the
present time, by estimating the amount of matter through a determination of
the Ly-$\alpha$ optical depth, thus providing one of the most stringent ways
of testing NCDM models.

Constraints from the Ly-$\alpha$ flux power spectrum on the DM properties are
usually given as a lower bound on the WDM mass parameter, $m_{\text{WDM}}$,
used as a reference. Given $m_{\text{WDM}}$, the WDM phase space is
characterized by a single quantity: $T_{\text{WDM}}$. In spite of our
notation, this quantity, which decreases with time, is not a temperature,
stricto sensu, since we assume the WDM species not to be in thermal
equilibrium at recombination and later times. Such a DM candidate is assumed
to have achieved a state of thermal equilibrium at some earlier time in the
evolution of the Universe and would have subsequently decoupled later on, as
it happens e.g. for neutrinos in the SM. Indeed, a good benchmark scenario for
WDM, which we will assume from now on, is a fermionic DM candidate with two
degrees of freedom having a Fermi-Dirac distribution. In this case the WDM
relic density can be related to its mass and “temperature” $T_{\text{WDM}}$ by

$\Omega_{\text{WDM}}h^{2}\simeq\left(\dfrac{m_{\text{WDM}}}{94\leavevmode\nobreak\
\text{eV}}\right)\left(\dfrac{T_{\text{WDM}}}{T_{\nu}}\right)^{3}\,,$ (2.16)

where $T_{\nu}=(4/11)^{1/3}\,T$ is the neutrino temperature as expected in the
SM after $e^{+}e^{-}$ annihilations, assuming instantaneous decoupling,
expressed as a function of the photon temperature $T$. As usual, $h$ denotes
the reduced Hubble constant, defined by the relation $H_{0}\equiv
100\,h\,\text{km}\,\text{s}^{-1}\,\text{Mpc}^{-1}$.

Assuming that the WDM saturates the DM density determined by Planck [47], a
numerical evaluation of the free-streaming horizon in Eq. (2.11) gives that
the cutoff in the linear matter power spectrum occurs at

$k_{H}(a=1)\,\simeq\,3.5\,h\,\text{Mpc}^{-1}\,.$ (2.17)

for $m_{\text{WDM}}=1$ keV. As shown in [56, 57, 39], an analytical fit for
the transfer function in the WDM case is given by

$\mathcal{T}(k)=\left(1+(\alpha k)^{2\nu}\right)^{-5/\nu}\,,$ (2.18)

with the parameters

$\nu=1.12\quad\text{and}\quad\alpha=0.24\left(\dfrac{}{1\leavevmode\nobreak\
\text{keV}}\dfrac{T_{\text{WDM}}}{T_{\nu}}\right)^{-0.83}\left(\dfrac{\Omega_{\text{WDM}}h^{2}}{0.12}\right)^{-0.16}\text{Mpc}\,.$
(2.19)

Importantly, these fitting parameters are independent of the standard
cosmological parameters (other than the DM abundance, $\Omega_{\text{WDM}}$).
The non-observation of a cutoff in actual data for the matter power spectrum
can be translated into a constraint on the WDM mass. A recent analysis [42]
gives a bound $m_{\text{WDM}}>5.3$ keV at 95% C.L., while the reference [43]
derived a less stringent bound $m_{\text{WDM}}>1.9$ keV at 95% C.L., by
claiming a more conservative treatment of thermal history for the
intergalactic medium. In the following we will take $m_{\text{WDM}}>3$ keV as
a reference but allow our results to be translated for a different value.

The most up-to-date lower bounds on the WDM mass from Ly-$\alpha$ data (1.1)
have been obtained using the medium resolution X-shooter spectrographic
observations of the intermediate redshift ($z:3-4.2$) XQ-100 sample of quasars
[58, 59] and the higher-resolution, higher-redshift ($z:4.2-5.4$) data from
the HIRES/MIKE spectrographs [60, 61]. These data can be used in combination
with probes of the matter power spectrum at smaller comoving scales
($k<(\rm{km/s})^{-1}$) via Lyman-$\alpha$ data, in particular from the Baryon
Oscillation Spectroscopic Survey (BOSS) of the Sloan Digital Sky Survey (SDSS-
III) [62, 63]. For future prospects (including higher redshifts), potentially
allowing an enhanced sensitivity to the cutoff of the spectrum see [64].

In principle, it may seem reasonable to assume that in order to compare the
expected matter power spectrum for a (more general, non Fermi-Dirac) NCDM
cosmology to the WDM case, it should be essential to take into account the
non-linear behaviour of the DM density field on the small cosmological
distances ($1-100$ Mpc) probed by Ly-$\alpha$ data. Performing such a
comparison requires costly N-body simulations for each possible NCDM case of
interest. Nevertheless, the authors of Ref. [65, 66] have performed a large
set of N-body simulations of models featuring an ample variety of transfer
functions, confronting the resulting power spectra to Ly-$\alpha$ data, and
concluding that all the models that are ruled out can also be rejected by
doing a simpler, linear analysis.

As we will show in the next sections, the shape of the linear power spectra of
the various NCDM models we consider turns out to be very similar to the one
for WDM, in spite of having, in some cases, notable differences at the level
of the phase space. Therefore, we can translate directly the WDM Ly-$\alpha$
bounds by computing, numerically, the linear transfer functions for our NCDM
models using a Boltzmann code, such as CLASS [44, 45], and comparing the
result with the linear transfer function in the WDM case.

The shape of the transfer function at the scales relevant for the change
induced in the matter power spectrum by WDM free-streaming can also be probed
by comparing the number of satellite galaxies of the Milky Way with N-body
simulations [67, 68]. This method gives a bound on the WDM mass that is
complementary and comparable to those obtained from Ly-$\alpha$ data. The
initial conditions for these N-body simulations were set in [67, 68] to mimic
the (linear) transfer function (2.18). Assuming that the formation of
satellite galaxies only depends on the nature of the DM through the linear
transfer function, we can also map these WDM mass bounds into constraints on
NCDM models that feature different distribution functions, just as we do with
Ly-$\alpha$ bounds.

###  2.3 Analytical rescaling and generalized phase space distribution

Let us now consider a NCDM model for which, by assumption, $w\ll 1$ is the
only quantity needed to characterize the cutoff in the linear transfer
function. Then, we can estimate the bound on $w$ from Ly-$\alpha$ by finding
the value of $m_{\text{DM}}$ such that

$\displaystyle w(m_{\text{DM}})=w_{\text{WDM}}(m_{\text{WDM}})\,.$ (2.20)

The bound is obtained by assuming that the cutoff scale of the linear matter
power spectrum for WDM can be translated to that of NCDM equating the
equations of state. A correspondence between two NCDM scenarios (a sterile
neutrino and a particle that decouples while being relativistic) was proposed
for the first time (to our knowledge) in Ref. [69]. By equating the power
spectra, the authors found a relation between these two scenarios, which
possess distribution functions with the same analytical expression but with
different parameters. A similar matching procedure using the mean square of
the DM velocity was proposed in [70] and extended in [71, 72] for several
freeze-in models. As we show below, our (generalized) matching relation can be
applied to a wide variety of NCDM scenarios, even for those in which thermal
equilibrium is not established before DM decoupling. From Eq. (2.13) we can
write $w_{\text{WDM}}$ as

$w_{\text{WDM}}(a)\,\simeq\,6\times
10^{-15}\,a^{-2}\,\left(\dfrac{\text{keV}}{m_{\text{WDM}}}\right)^{8/3}\,,$
(2.21)

implying that the bound on $m_{\text{WDM}}\sim$ keV from Ly-$\alpha$ [39, 42,
43] translates into

$\displaystyle w_{\text{WDM}}(a=1)\lesssim 10^{-15}\,,$ (2.22)

showing that DM is indeed very cold. It is worth emphasizing that the
constraint from Eq. (2.22) corresponds to a constraint at recombination of
$w_{\text{WDM}}(a\sim 10^{-3})\lesssim 10^{-9}$ whereas analyses based on CMB
data constrain this value only at the level of $w_{\text{WDM}}(a\sim
10^{-3})\lesssim 10^{-4}$ [52, 73].

For instance, a typical WIMP with a mass $m_{\rm DM}=100\leavevmode\nobreak\
\text{GeV}$ that decoupled at a freeze-out temperature $T_{F}\simeq m_{\rm
DM}/20$ inherits a Maxwell-Boltzmann distribution after decoupling of the form

$f(p,t)\,=\,\dfrac{g_{\rm
DM}}{(2\pi)^{3}}\exp\left[\dfrac{-p^{2}a(t)^{2}}{2m_{\rm
DM}a_{F}^{2}T_{F}}\right]\,.$ (2.23)

In the non-relativistic limit, the energy and pressure densities can be
evaluated analytically and $w$ can be expressed as

$w(a)\,\simeq\,\dfrac{a_{F}^{2}}{a^{2}}\dfrac{T_{F}}{m_{\rm
DM}}\,\simeq\,10^{-29}\left(\dfrac{1}{a^{2}}\right)\left(\dfrac{20\,T_{F}}{m_{\rm
DM}}\right)\left(\dfrac{100\leavevmode\nobreak\ \text{GeV}}{m_{\rm
DM}}\right)^{2}\left(\dfrac{100}{g_{*}^{F}}\right)^{2/3}\,,$ (2.24)

where $g_{*}^{F}$ denotes the effective number of degrees of freedom. This
value for $w$ is several orders of magnitude lower than the typical value
constrained by Ly-$\alpha$. For our NCDM case, Eq. (2.20) leads to

$m_{\text{DM}}\,=\,m_{\text{WDM}}\left(\dfrac{T_{\star}}{T_{\text{WDM}}}\right)\sqrt{\dfrac{\langle
q^{2}\rangle}{\langle
q^{2}\rangle_{\text{WDM}}}}\simeq\,7.56\leavevmode\nobreak\
\text{keV}\,\left(\dfrac{m_{\text{WDM}}}{3\leavevmode\nobreak\
\text{keV}}\right)^{4/3}\left(\dfrac{T_{\star}}{T}\right)\sqrt{\langle
q^{2}\rangle}\,.$ (2.25)

Alternatively, the bound can be expressed in terms of the mean momentum at the
present time, $\langle p\rangle_{0}=\langle q\rangle\,T_{\star}$ where
$\langle q\rangle$ is defined in Eq. (2.14), giving

$m_{\text{DM}}\,\simeq\,7.56\leavevmode\nobreak\
\text{keV}\,\left(\dfrac{m_{\text{WDM}}}{3\leavevmode\nobreak\
\text{keV}}\right)^{4/3}\left(\dfrac{\langle
p\rangle_{0}}{T_{0}}\right)\dfrac{\sqrt{\langle q^{2}\rangle}}{\langle
q\rangle}\,.$ (2.26)

As we will show, most of the NCDM cases discussed in this paper can be well
described with a generalized phase space distribution of the form

$f(q)\,\propto\,q^{\alpha}\,\exp{\left(-\beta\,q^{\gamma}\right)}\,,$ (2.27)

with constant $\alpha>-3$ and $\beta,\gamma>0$ as required for the DM number
density to be finite. For this distribution the normalized $n$-th moment
(2.14) is

$\langle
q^{n}\rangle\,=\,\beta^{\frac{2-n}{\gamma}}\,\dfrac{\Gamma\left(\frac{1+n+\alpha}{\gamma}\right)}{\Gamma\left(\frac{3+\alpha}{\gamma}\right)}\,.$
(2.28)

The rescaling of the mass reproducing the same cutoff as the WDM case then
gives

$m_{\text{DM}}\,\simeq\,7.56\leavevmode\nobreak\
\text{keV}\,\left(\dfrac{m_{\text{WDM}}}{3\leavevmode\nobreak\
\text{keV}}\right)^{4/3}\left(\dfrac{\langle
p\rangle_{0}}{T_{0}}\right)\,\sqrt{\dfrac{\Gamma\left(\frac{3+\alpha}{\gamma}\right)\,\Gamma\left(\frac{5+\alpha}{\gamma}\right)}{\Gamma^{2}\left(\frac{4+\alpha}{\gamma}\right)}}\,,$
(2.29)

which does not depend explicitly on $\beta$.666However, the mean momentum
$\langle p\rangle_{0}$ depends actually on this quantity. As we will show
later for the various examples we consider, this allows to translate any bound
from Ly-$\alpha$ on $m_{\text{WDM}}$ to a bound for a given NCDM model on
$m_{\text{DM}}$ provided that the DM phase space distribution can be well
described by (2.27), with good precision and without requiring a numerical
computation of the power spectrum.777Notably, in the case we consider for
which (2.27) does not apply ($n=6$ thermal freeze-in, for example), we can
still find an analytical form for the constraint using the more general form
of Eq. (2.25).

Provided that the NCDM equation-of-state parameter evolves as $w(a)\simeq
w_{0}\,a^{-2}$ at redshift $z<10^{6}$, for which the wavenumbers $k$ relevant
for Lyman-$\alpha$ data enters the horizon, the NCDM power spectrum should
exhibit the same features as the WDM power spectrum at first order in $w$. A
different $a$-dependence of $w(a)$ would affect our matching procedure. For
instance, in cannibalistic dark matter scenarios, the equation-of-state
parameter behaves as $w(a)\propto 1/\log a$ when number-changing processes are
active and $w(a)\propto a^{-2}$ in the non-relativistic regime at later times
$z<10^{3}$[36]. In this case it has been shown that the NCDM power spectrum
can be matched to a WDM one with a good precision[36], by introducing a
similar matching procedure up to a correction factor of order one[74].

Figure 3: Relative difference between the transfer functions of the various
NCDM models considered in this work and the WDM case (assuming
$m_{\text{WDM}}=3$ keV). The scale $k_{1/2}^{\text{WDM}}$ defined in Section
2.3 is represented as a black vertical dashed line. Our matching procedure
performs at the level of $\sim 3\%$ or better at this scale for most models.
The notation for the modulus decay cases is introduced in Section 3.2.1 (see
Fig. 8). For thermal particle decays, FD and BE stand for the Fermi-Dirac and
Bose-Einstein distributions for the decaying particle, respectively (see Fig.
11). For non-thermal decays, R and NR denote a relativistic or a non-
relativistic decaying particle, respectively (see Eq. 4.14). For thermal and
non-thermal freeze-in, the parameter $n$ has been introduced in Eq. (1.5).

The numerical precision achieved using our rescaling procedure is represented
in Fig. 3, which shows the relative difference between the various NCDM
transfer functions considered in this work and the one for the WDM case. The
transfer functions are computed numerically with CLASS and using our rescaling
procedure to for a given WDM mass, which we assume to be $m_{\text{WDM}}=3$
keV in Fig. 3. The NCDM transfer functions match accurately, mostly with a
precision below the percent level, the WDM transfer function for
$k<20\,h\,\text{Mpc}^{-1}$. The precision decreases for larger modes $k$.

In order to estimate the difference on the cutoff scale expected between NCDM
and WDM using our procedure, with a quantity more relevant for observational
data, we represent in Fig. 3 the scale $k_{1/2}^{\text{WDM}}$ defined such
that
$\mathcal{P}_{\text{WDM}}(k_{1/2}^{\text{WDM}})=(1/2)\mathcal{P}_{\Lambda\text{CDM}}(k_{1/2}^{\text{WDM}})$
for a given WDM mass. Fig. 3 shows that at $k_{1/2}^{\text{WDM}}$, our
rescaling procedure allows to achieve a $3\%$ difference or better on the NCDM
transfer function, relative to the WDM case, for most of our scenarios. The
least precise cases achieve a $\sim 10\%$ difference at
$k_{1/2}^{\text{WDM}}$. These correspond to DM production from the decay of a
non-thermal relativistic particle and DM production via thermal freeze-in with
$n=6$.

An accurate estimate of the transfer function is numerically more challenging
due to the specific shape of the phase space distributions in these cases. For
this reason, we believe that the larger relative difference displayed in Fig.
3 for these cases can be partially attributed to the requirement of having a
reasonable computation time, at the price of a limited precision.

##  3 Decay of a classical condensate

We consider as the first application of our formalism the production of DM
from the perturbative decay of a classical, spatially homogeneous, oscillating
condensate. As a first example, we study the decay of the inflaton field into
DM during reheating, assuming all other DM interactions can be neglected. We
then consider the decay of a modulus field, a scalar field present in the
early Universe with a non-vanishing vacuum misalignment: a displacement from
its post-inflationary global minimum, which leads to a subsequent epoch of
oscillations about this minimum. We explore the scenario in which the
oscillations of the modulus dominate the energy density of the Universe at
late times and, also, the case in which they are subdominant to the inflaton
or radiation background. In all cases we find the non-thermal DM phase space
distributions, and the corresponding Ly-$\alpha$ bounds on the DM mass.

###  3.1 Perturbative inflaton decay

Let us assume that the production of a DM particle $\chi$ proceeds through the
two-body decay of the inflaton field during reheating, i.e. through a process
of the form $\phi\rightarrow\chi+\psi$. The rest frame decay rate for this
process is given by $\Gamma_{\phi\rightarrow\chi\psi}={\rm Br}_{\chi}$, where
${\rm Br}_{\chi}$ denotes the branching ratio to $\chi$, and is the total
decay rate of the inflaton. We assume the coupling of $\chi$ with $\phi$ is
sufficiently weak to disregard the re-population of $\phi$ from inverse
decays.888This is an assumption which is justified a posterori by requiring
the generated DM density to match the observed relic abundance.  Moreover, as
in all other cases, we assume that the couplings of $\chi$ to the visible
sector or to itself are not strong enough to bring it to kinetic and/or
chemical equilibrium, and may therefore be disregarded. It must be emphasized
that, for simplicity, in each of the cases discussed in this paper we assume
that $100\%$ of the DM relic abundance is produced by a single mechanism. In
addition, we assume that the DM particles do not have significant interactions
between them.999As highlighted recently in [36, 75] self-interactions could
affect the power spectrum, in particular for cases with light DM masses.
Moreover, we neglect thermal effects that would give rise to a subdominant
contribution for UV freeze-in but have been shown to alter the produced DM
phase space distribution in IR-dominated freeze-in specific scenarios [76].

In order to apply the procedure described in Section 2 for mapping the WDM
bound on $m_{\text{WDM}}$ into a bound on the mass of $\chi$, we must first
determine the form of the phase space distribution $f_{\chi}$ generated from
decays of the inflaton field, by solving the Boltzmann transport equation

$\frac{\partial f_{\chi}}{\partial t}-H|\boldsymbol{p}|\frac{\partial
f_{\chi}}{\partial|\boldsymbol{p}|}\,=\,\mathcal{C}[f_{\chi}(|\boldsymbol{p}|,t)]\,,$
(3.1)

where $[f_{\chi}]$ denotes the collision term, determined by the inflaton-DM
interaction. In Appendix A we provide the general form of this collision term,
as well as the general solution of (3.1) in the absence of inverse processes
and in the free-streaming limit.

####  3.1.1 DM phase space distribution

Under the assumptions discussed above, the decay of the inflaton to $\chi$
will be perturbative. If this is true for all its decay channels, then $\phi$
is, on average, spatially homogeneous, and the phase space distribution may be
written as $f_{\phi}(k,t)=(2\pi)^{3}n_{\phi}(t)\delta^{(3)}(\boldsymbol{k})$,
with $n_{\phi}$ the instantaneous inflaton number density. Disregarding
inverse decays, the collision term for the transport equation that determines
the distribution function for $\chi$ takes the form

$\displaystyle\mathcal{C}[f_{\chi}(p,t)]\;$
$\displaystyle=\;\frac{1}{2p_{0}}\int\frac{\mathop{}\\!\mathrm{d}^{3}{\boldsymbol{k}}}{(2\pi)^{3}2k_{0}}\frac{g_{\psi}\mathop{}\\!\mathrm{d}^{3}{\boldsymbol{p}}_{\psi}}{(2\pi)^{3}2p_{\psi
0}}(2\pi)^{4}\delta^{(4)}(k-p-p_{\psi})$ $\displaystyle\hskip
145.0pt\times|\mathcal{M}|^{2}_{\phi\rightarrow\chi\psi}f_{\phi}(k)\left(1\pm
f_{\chi}(p)\pm f_{\psi}(p_{\psi})\right)$ (3.2) $\displaystyle=\;\frac{\pi
n_{\phi}}{4m_{\phi}p_{0}}\int_{\rm
RF}\frac{g_{\psi}\mathop{}\\!\mathrm{d}^{3}{\boldsymbol{p}}_{\psi}}{p_{\psi
0}}\delta(m_{\phi}-p_{0}-p_{\psi
0})\delta^{(3)}(\boldsymbol{p}+\boldsymbol{p}_{\psi})|\mathcal{M}|^{2}_{\phi\rightarrow\chi\psi}\left(1\pm
f_{\chi}(p)\pm f_{\psi}(p_{\psi})\right)$
$\displaystyle=\;\frac{2\pi^{2}}{g_{\chi}\varepsilon^{2}_{\psi}}n_{\phi}\Gamma_{\phi\rightarrow\chi\psi}(1\pm
f_{\chi}(p_{0})\pm
f_{\psi}(\varepsilon_{\psi}))\delta(p_{0}-\varepsilon_{\psi})\,,$ (3.3)

where notations and conventions are detailed in the appendix. Here
$\varepsilon_{\psi}=(m_{\phi}^{2}+m_{\psi}^{2}-m_{\text{DM}}^{2})/2m_{\phi}$
denotes the energy of the daughter particle. The collision term can be further
simplified in the limit when $m_{\text{DM}},m_{\psi}\ll m_{\phi}$, so that
$p_{0}\simeq|\boldsymbol{p}|=p$, and if the quantum statistics of the decay
products can be neglected.101010This is ensured provided that the effective
coupling
$y_{\chi}\equiv(8\pi\Gamma_{\phi\rightarrow\chi\psi}/m_{\phi})^{1/2}\ll
10^{-5}$. If this is the case we can simply write

$\mathcal{C}[f_{\chi}(p,t)]\;=\;\frac{8\pi^{2}}{g_{\chi}m_{\phi}^{2}}n_{\phi}\Gamma_{\phi\rightarrow\chi\psi}\delta(p-m_{\phi}/2)\,,$
(3.4)

Substitution of this collision term into the transport equation (3.1) yields
an equation that has an exact solution in terms of the Hubble parameter $H$
and the inflaton occupation number [49, 77],

$\displaystyle f_{\chi}(p,t)\;$
$\displaystyle=\;\frac{16\pi^{2}\Gamma_{\phi\rightarrow\chi\psi}n_{\phi}(\hat{t})}{g_{\chi}m_{\phi}^{3}H(\hat{t})}\theta(t-\hat{t})\,,$
(3.5)

where $\hat{t}$ is the solution to the equation

$\frac{a(t)}{a(\hat{t})}=\frac{m_{\phi}}{2p}\,.$ (3.6)

In order to obtain a closed form for $f_{\chi}$ we need to solve for the
inflaton number density and the expansion rate. This can be achieved by
integrating the Friedmann-Boltzmann system of equations

$\displaystyle\dot{\rho}_{\phi}+3H\rho_{\phi}+\Gamma_{\phi}\rho_{\phi}\;$
$\displaystyle=\;0\,,$ (3.7)
$\displaystyle\dot{\rho}_{r}+4H\rho_{r}-\Gamma_{\phi}\rho_{\phi}\;$
$\displaystyle=\;0\,,$ (3.8) $\displaystyle\rho_{\phi}+\rho_{r}\;=\;3H^{2}$
$\displaystyle M_{P}^{2}\,,$ (3.9)

where the reduced Planck is $=1/\sqrt{8\pi\,G}$ (being $G$ Newton’s
gravitational constant) and where we denote by $\rho_{\phi}$ and $\rho_{r}$
the energy densities of the inflaton condensate and that of its relativistic
decay products, respectively. Note that (3.8) is nothing but the integrated
version of the transport equation (3.1) for an ultrarelativistic species with
${\rm Br}_{r}=1$. Straightforward integration gives [78]

$n_{\phi}(t)\;=\;\frac{\rho_{\phi}(t)}{m_{\phi}}\;=\;\frac{\rho_{\rm
end}}{m_{\phi}}\left(\frac{a(t)}{a_{\rm
end}}\right)^{-3}e^{-\Gamma_{\phi}(t-t_{\rm end})}\,,$ (3.10)

where the sub-index “end” denotes quantities at the end of inflation. For
$t_{\rm end}\ll t\ll\Gamma_{\phi}^{-1}$ the exponential in the previous
expression can be disregarded: the Universe is dominated by the matter-like
oscillations of $\phi$. Therefore, we may also approximate $a\propto t^{2/3}$,
and $\hat{t}\simeq(2p/m_{\phi})^{3/2}t$. Substitution into (3.5) yields the
following expression for the phase space distribution of $\chi$ well before
the end of reheating at $t_{\rm reh}\simeq\Gamma_{\phi}^{-1}$,

($t\ll t_{\rm reh}$) $\displaystyle\begin{aligned}
f_{\chi}(p,t)\;&=\;\frac{24\pi^{2}{\rm
Br}_{\chi}\Gamma_{\phi}}{g_{\chi}m_{\phi}^{3}}\left(\frac{m_{\phi}}{2p}\right)^{3/2}t\,n_{\phi}(t)\,\theta(m_{\phi}/2-p)\\\
&\simeq\;\frac{24\pi^{2}n_{\chi}(t)}{g_{\chi}m_{\phi}^{3}}\left(\frac{m_{\phi}}{2p}\right)^{3/2}\,\theta(m_{\phi}/2-p)\,.\end{aligned}$
(3.11)

Here we have approximated the number density of decay products as

$n_{\chi}(t)\;\simeq\;{\rm Br}_{\chi}\frac{\rho_{\rm
end}}{m_{\phi}}\left(1-e^{-\Gamma_{\phi}(t-t_{\rm
end})}\right)\left(\frac{a(t)}{a_{\rm end}}\right)^{-3}\,.$ (3.12)

obtained by counting the quanta produced from inflaton decay. Note the
consistency of (3.11) with the defining relation (A.3) between $f_{\chi}$ and
$n_{\chi}$.

The distribution (3.11) will come handy for our study of non-thermal freeze-in
in Section 5.2. For our present purposes, though, this distribution is
incomplete, as it lacks the high momentum tail that will be generated when the
inflaton energy density begins to get exhausted. We must therefore extend
(3.11) beyond the end of reheating. As a first approximation, we evaluate
(3.5) at $t_{\rm reh}=\Gamma_{\phi}^{-1}$, the moment of time at which the
energy density in $\phi$ is approximately equal to that in radiation,
$\rho_{\phi}\simeq\rho_{r}$, and where $H_{\rm reh}\simeq 2\Gamma_{\phi}/3$.
With the reheating temperature given by

$T_{\rm reh}\;=\;\left(\frac{30\rho_{\rm rad}}{\pi^{2}g_{*s}^{\rm
reh}}\right)^{1/4}\,,$ (3.13)

where denotes the effective number of relativistic degrees of freedom for
entropy, we can substitute into (3.5) to obtain

$f_{\chi}(p,t_{\rm reh})\;\simeq\;\frac{4\pi^{4}{\rm Br}_{\chi}g_{*s}^{\rm
reh}}{5g_{\chi}}\left(\frac{T_{\rm
reh}}{m_{\phi}}\right)^{4}\left(\frac{m_{\phi}}{2p}\right)^{3/2}e^{1-(2p/m_{\phi})^{3/2}}\theta(m_{\phi}/2-p)\,.$
(3.14)

Naively, this distribution would evolve at later times simply in accordance to
(A.6). However, the production of entropy from inflaton decay does not
suddenly stop at $t_{\rm reh}$, but continues for some time into the radiation
domination era. The continuous transition $w=0\rightarrow 1/3$ makes the
analytical estimation of $f_{\chi}$ beyond $t_{\rm reh}$ complicated, although
not impossible (see e.g. [79, 48]). Nevertheless, Eqs. (3.5) and (3.10) make
an estimate of the shape of the tail of the distribution straightforward.
During radiation domination $a\propto t^{1/2}$, implying that for momenta
which satisfy the relation $t_{\rm reh}\ll\hat{t}=(2p/m_{\phi})^{2}t$, the
time-dependence of $n_{\phi}$ yields

($t_{\rm reh}\ll(2p/m_{\phi})^{2}t$) $\displaystyle
f_{\chi}(p,t)\;\propto\;\exp\left[\left(\frac{2p}{m_{\phi}}\right)^{2}\frac{t}{t_{\rm
reh}}\right]\theta(m_{\phi}/2-p)\,,$ (3.15)

i.e. a Gaussian tail.

A better approximation for $f_{\chi}(p,t)$ beyond the end of reheating can be
constructed by solving numerically the Friedmann-Boltzmann system (3.7)-(3.9)
together with (3.1) with collision term (3.4). This solution is shown as the
continuous black curve in Fig. 4, in the form of the rescaled distribution
$\bar{f}_{{\rm R}}$, defined through the relation

Figure 4: The rescaled distribution function $\bar{f}_{{\rm R}}$, defined in
(3.16), as a function of the rescaled momentum $q$, for DM produced from
inflaton decay. Solid, black: the numerically computed result. Dashed-dotted,
blue: the analytical result (3.14) without the Heaviside function. Dashed,
orange: the phenomenological fit (3.18). The part of the distribution for
which $q<1$ is populated during $t<t_{\rm reh}$. The part of the distribution
for which $q>1$ is populated during $t>t_{\rm reh}$.

$f_{\chi}(p,t)\,\mathop{}\\!\mathrm{d}^{3}\boldsymbol{p}\;=\;\frac{4\pi^{4}{\rm
Br}_{\chi}g_{*s}^{\rm reh}}{5g_{\chi}}\left(\frac{T_{\rm
reh}}{m_{\phi}}\right)^{4}\,\left(\frac{a_{0}}{a(t)}\right)^{3}T_{\star}^{3}\,\bar{f}_{\rm
R}(q)\,\mathop{}\\!\mathrm{d}^{3}\boldsymbol{q}\,.$ (3.16)

Here $q$ is defined as in (2.2), and in this scenario

$T_{\star}\;=\;\frac{m_{\phi}}{2}\frac{a_{\rm
reh}}{a_{0}}\;=\;\left(\frac{g_{*s}^{0}}{g_{*s}^{\rm
reh}}\right)^{1/3}\frac{m_{\phi}}{2T_{\rm reh}}\,T_{0}\,.$ (3.17)

The numerical solution was computed at $t=50t_{\rm reh}$, well beyond the
matter-radiation equality that signals the end of reheating. At this time the
universe is dominated by radiation, and the production of entropy from
inflaton decay has ceased. The particle population that was produced during
$t<t_{\rm reh}$ occupies the distribution at $q<1$, while the population
created during $t>t_{\rm reh}$ corresponds to the $q>1$ tail. Shown in Fig. 4
is also the analytical solution (3.14), ignoring the Heaviside cutoff at
$q=1$. As expected, this expression accurately describes the distribution at
small momenta, $\bar{f}_{\rm R}\propto q^{-3/2}$, but the tail is not matched.
Given that we expect the large momentum regime to be described by (3.15), we
also show in the figure, as an orange dashed curve, a fitting function that
mimics the low- and high-energy behavior of the distribution,

$\bar{f}_{{\rm R}}(q)\;\simeq\;2.28\,q^{-3/2}e^{-0.74q^{2}}\,.$ (3.18)

This approximation is of the form (2.27), and provides an excellent fit to the
exact form of $\bar{f}_{\rm R}$. Note the seeming mismatch between the ratio
$(t/t_{\rm reh})^{1/2}$ and the ratio $a(t)/a_{\rm reh}$ through which $q$ is
defined, quantified by the factor 0.74 in the exponent. This is due to the
relatively complicated dependence of the scale factor on time in the matter-
radiation transition at the end of reheating, affecting the high-energy tail
of the distribution.

####  3.1.2 Power spectrum and Ly-$\alpha$ constraints

With the phase space distribution for DM produced from direct inflaton decay,
we can now make use of Eq. (2.25) to map the WDM Ly-$\alpha$ constraints on
the DM mass for this scenario. Straightforward calculation gives the following
rescaling of the bound on the DM mass,

$\displaystyle
m_{\text{DM}}\;\gtrsim\;\left(\dfrac{m_{\text{WDM}}}{3\leavevmode\nobreak\
\text{keV}}\right)^{4/3}$ $\displaystyle\left(\dfrac{106.75}{g_{*s}^{\rm
reh}}\right)^{1/3}$ $\displaystyle\times\left(\dfrac{m_{\phi}}{3\times
10^{13}\leavevmode\nobreak\
\text{GeV}}\right)\left(\dfrac{10^{10}\leavevmode\nobreak\
\text{GeV}}{T_{\text{reh}}}\right)\begin{cases}3.78\leavevmode\nobreak\
\text{MeV}\,,\leavevmode\nobreak\ &{\rm Numerical}\,,\\\
4.11\leavevmode\nobreak\ \text{MeV}\,,\leavevmode\nobreak\ &{\rm
Analytical}\,,\\\ 3.79\leavevmode\nobreak\ \text{MeV}\,,\leavevmode\nobreak\
&{\rm Fit}\,.\end{cases}$ (3.19)

The numerical, analytical and fit approximations correspond to the numerically
computed distribution shown in Fig. 4, to (3.14), and to (3.18), respectively.
For low reheating temperatures the bound on the NCDM mass becomes
significantly larger than that for WDM. This can be understood by fixing the
inflaton mass and decreasing progressively the reheating temperature. The bulk
of DM is produced around the reheating temperature with typical momentum
$p\sim m_{\phi}/2$, regardless of the radiation temperature. Reducing the
reheating temperature therefore prevents the momentum of the DM particle from
redshifting too much, resulting in a hotter spectrum at the present time than
that expected for large reheating temperatures.

Figure 5: Linear transfer function for the scenario in which DM is produced by
inflaton decay, assuming the numerical, analytical or fitted phase space
distributions described in Section 3.1, by taking the mass estimated in (3.19)
and identical reheating temperatures, $T_{\rm reh}=10^{10}\,{\rm GeV}$. The
transfer function for the WDM case is shown for comparison with a dashed black
line. We depicted as well the transfer function for numerical, analytical or
fitted phase space distributions with identical masses and reheating
temperature, $T_{\rm reh}=10^{12}\,{\rm GeV}$.

Fig. 5 shows the form of the transfer function for the matter power spectrum,
as computed with CLASS [44, 45]. Depicted are the results for the numerical,
analytical and fit approximations to $f_{\chi}$. The rightmost set of curves
shows the form of each $\mathcal{T}(k)$ for a reheating temperature $T_{\rm
reh}=10^{10}\,{\rm GeV}$, and masses given by Eq. (3.19). The overlap of all
three curves with each other, and with the reference WDM transfer function,
demonstrates the validity of our method for this DM production mechanism. At
$k_{1/2}^{\rm WDM}$ the relative difference between WDM and the numerical
result is in particular smaller than $10^{-3}$, c.f. Fig. 3. The leftmost
cluster of curves shows the form of $\mathcal{T}(k)$ for the three
approximations for a larger reheating temperature, $T_{\rm reh}=10^{12}\,{\rm
GeV}$, assuming a mass of $10\,{\rm keV}$. These curves do not overlap with
the WDM bound, and they differ slightly between each other, albeit the
agreement between the numerical and fit cases is still excellent.

####  3.1.3 Relic density and phenomenology

We now discuss the phenomenological implications of a lower bound on a light
DM particle produced from inflaton decay. Given a reheating temperature and a
DM mass, the normalization of the distribution function is determined by the
value of the present DM fraction $\Omega_{\chi}=\rho_{\chi}/\rho_{c}$, where
$\rho_{c}\simeq 1.05\times 10^{-5}h^{2}{\rm GeV\,cm}^{-3}$ is the present
critical density of the Universe [80]. Integration of (3.16) at $t\gg t_{\rm
reh}$ gives

$n_{\chi}(t)\;\simeq\;0.70\pi^{2}{\rm Br}_{\chi}g_{*s}^{\rm
reh}\left(\frac{T_{\rm
reh}}{m_{\phi}}\right)^{4}\,\left(\frac{a_{0}}{a(t)}\right)^{3}T_{\star}^{3}\,,$
(3.20)

which in turn yields

$\Omega_{\chi}h^{2}\;\simeq\;0.1\left(\frac{{\rm Br}_{\chi}}{5.5\times
10^{-4}}\right)\left(\frac{m_{\text{DM}}}{1\,{\rm
MeV}}\right)\left(\frac{T_{\rm reh}}{10^{10}\,{\rm
GeV}}\right)\left(\frac{3\times 10^{13}\,{\rm GeV}}{m_{\phi}}\right)\,.$
(3.21)

Combining the bounds on the DM mass (3.19) and on the relic abundance (3.21),
the following constraint can be derived for the branching ratio of the decay
of the inflaton into dark matter,

${\rm Br}_{\chi}\;\lesssim\;1.5\times 10^{-4}\,\left(\frac{g_{*s}^{\rm
reh}}{106.5}\right)^{1/3}\left(\frac{3\,{\rm keV}}{m_{\rm
WDM}}\right)^{4/3}\,.$ (3.22)

Note the universality of this bound: it is independent of the inflaton mass
and the reheating temperature. As mentioned earlier, such a limit will apply
even in the absence of tree-level couplings between the inflaton and DM.
Assuming a dominant fermionic decay channel of the inflaton, with these decay
products in turn coupled to DM through an effective interaction of the
following form,

$\mathcal{L}\;=\;y\phi\bar{f}f+\frac{1}{\Lambda^{2}}\bar{f}f\bar{\chi}\chi\,,$
(3.23)

(which could arise from the exchange of a massive field with mass
$\sim\Lambda$), a non-vanishing decay rate for the
$\phi\rightarrow\bar{\chi}\chi$ process is induced at 1-loop [46],

$\Gamma_{\phi\rightarrow\bar{\chi}\chi}\;\simeq\;\frac{y^{2}}{128\pi^{5}}\left(1+\frac{\pi^{2}}{4}\right)\frac{m_{\phi}^{5}}{\Lambda^{4}}\,,$
(3.24)

corresponding to ${\rm
Br}_{\chi}=\frac{1}{16\pi^{4}}(1+\frac{\pi^{2}}{4})(\frac{m_{\phi}}{\Lambda})^{4}$.
Substitution into (3.22) reveals that

$\Lambda\;\gtrsim\;2\,m_{\phi}\left(\frac{106.5}{g_{*s}^{\rm
reh}}\right)^{1/12}\left(\frac{m_{\rm WDM}}{3\,{\rm keV}}\right)^{1/3}\,,$
(3.25)

a condition consistent with the form of the effective action (3.23), assumed
to be valid at all times during reheating.

We finish this section by emphasizing that the bounds (3.19) and (3.22) apply
for the perturbative decay of the inflaton $\phi$ while it oscillates about a
quadratic minimum. A different production mechanism, e.g. through perturbative
decay in a non-quadratic potential [81, 82, 83], or via non-adiabatic particle
production [84, 85, 86, 87, 88, 89], will lead to a different constraint on
${\rm Br}_{\chi}$.

###  3.2 Moduli decays

The inflaton is not necessarily the only scalar condensate that can decay in
the early Universe. In many BSM constructions, notably supersymmetric and
string SM extensions, a plethora of weakly-interacting unstable scalar fields,
collectively known as moduli, arise [90, 91, 92, 93, 94, 95]. During
inflation, these moduli can be excited away from the minima of their
potential, resulting in a posterior roll towards these minima. Depending on
the initial misalignment, and the masses of the moduli, the subsequent
oscillations about the minima may eventually dominate the energy density of
the Universe. The decay of these fields would then reheat the Universe at
temperatures below the inflationary reheating temperature, diluting any relics
produced earlier (such as DM) and the baryon asymmetry. This process would
also lead to deviations from the standard Big Bang Nucleosynthesis (BBN),
which is strongly constrained by the data, unless $T_{\rm reh}\gtrsim 1\,{\rm
MeV}$ [96, 97].

If a modulus $Z$ has a non-vanishing branching ratio to DM, the Ly-$\alpha$
bounds derived in the previous section can be mapped to its decay (provided
that $m_{Z}\gg m_{\rm DM}$) simply by replacing the inflaton mass and
reheating temperatures with their corresponding modulus values. In particular,
for the mass bound, we can write

$m_{\rm DM}\;\gtrsim\;12.6\,{\rm GeV}\left(\frac{m_{\rm WDM}}{3\,{\rm
keV}}\right)^{4/3}\left(\frac{106.75}{g_{*s}^{{\rm
reh},Z}}\right)^{1/3}\left(\frac{m_{Z}}{10\,{\rm
TeV}}\right)\left(\frac{1\,{\rm MeV}}{T_{{\rm reh},Z}}\right)\,,$ (3.26)

while (3.22) remains unchanged, except for the replacement $g_{*s}^{\rm
reh}\rightarrow g_{*s}^{{\rm reh},Z}$. Note that for moduli with masses
$m_{Z}\gtrsim 100\,{\rm TeV}$, the lower bound on the DM mass is $\gtrsim
100\,{\rm GeV}$, on the range of electroweak-scale DM candidates such as the
lightest neutralino. Moreover, the late decay of $Z$ would ensure that the
non-thermal phase space distribution (3.18) remains imprinted into this relic.
This is due to the fact that most of the DM is produced around $T_{{\rm
reh},Z}$, well below the corresponding thermal decoupling (freeze-out)
temperature. Fig. 6 shows the limit (3.26) in the mass vs. reheating
temperature plane, excluding the model-dependent $T_{{\rm reh},Z}>m_{Z}$
region. Note the wide range of values for $m_{\rm DM}$.

Figure 6: Ly-$\alpha$ constraint on the DM mass, as a function of the modulus
mass $m_{Z}$ and reheating temperature, in the case in which the oscillation
of $Z$ dominates the energy density of the Universe, leading to entropy
production upon its decay. The gray region corresponds to $T_{{\rm
reh},Z}>m_{Z}$, where in-medium and/or non-perturbative effects may determine
the decay of $Z$. At $m_{Z}=3\times 10^{13}\,{\rm GeV}$ the inflaton decay
scenario is recovered.

We emphasize that this kind of constraint must be accounted for in any
discussion regarding DM production in non-standard thermal histories, with an
intermediate matter-dominated epoch between the end of reheating and BBN [98,
99, 100, 101, 102, 103, 104]. For a sufficiently large branching ratio of $Z$
to $\chi$, this non-thermal production can dominate over DM freeze-out, which
would have occurred during the modified expansion history. We finally mention
that it is typical of the decay of a modulus into dark matter to occur in two
stages, $Z\rightarrow A\rightarrow\chi$, where $A$ is an intermediate unstable
particle, such as the gravitino. This scenario is studied in detail in Section
4.2. As we show there, although the phase space distribution of $A$ and $\chi$
differ noticeably in their shape, the rescaled bound on $m_{\rm DM}$ is only
corrected by an $\mathcal{O}(1)$ factor in some regimes.

Our main focus in this section is instead stabilized moduli: scalar
condensates that oscillate and subsequently decay in the early Universe, while
never dominating the energy budget of the Universe [105, 106, 107, 108, 109,
110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123]. Without
modifying inflation [124], this is typically achieved in model-building by
introducing additional interactions that rise the mass of the modulus,
increasing its decay rate, and by decreasing the amount of initial
misalignment.

It is important to realize that for a subdominant decaying scalar, the post-
inflationary background dynamics will be determined by either the oscillating
inflaton, or by its redshifting relativistic decay products. Therefore, it is
necessary to distinguish between three different scenarios: (a) the modulus
begins oscillating and decays during reheating, (b) the modulus begins
oscillating during reheating, but decays during radiation domination, or (c)
the modulus oscillates and decays during radiation domination. We now proceed
to determine the phase space distribution in all three cases, to subsequently
determine the Ly-$\alpha$ bounds and the corresponding phenomenologies. As we
discuss below, the observed DM abundance can be obtained from the decay of a
stabilized modulus when its energy density is much smaller than that of
radiation.

####  3.2.1 DM phase space distribution

Case a: Oscillation and decay during reheating
($m_{Z}>\Gamma_{Z}>\Gamma_{\phi}$)

We begin by studying the scenario in which the field $Z$ begins its
oscillations during the matter-dominated reheating, and fully decays before
the end of reheating. Given that we follow the decay of a classical
condensate, its distribution function will be of the form
$f_{Z}(k,t)=(2\pi)^{3}n_{Z}(t)\delta^{(3)}(\boldsymbol{k})$, where $n_{Z}$ is
the modulus number density (see Eq. (A.3)), and hence the DM distribution will
be given by (3.5), upon replacing $\phi\rightarrow Z$. The solution of Eq.
(3.6), necessary to determine the cosmic-time dependence of $f_{\chi}$, can be
found in a straightforward way, and is given by $\hat{t}=(2p/m_{Z})^{3/2}t$.
Moreover, the number density of the decaying $Z$ is found by integration of
(3.7), again replacing $\phi\rightarrow Z$,

$n_{Z}(t)\;=\;\frac{\rho_{\rm osc}}{m_{Z}}\left(\frac{a(t)}{a_{\rm
osc}}\right)^{-3}e^{-\Gamma_{Z}(t-t_{\rm osc})}\,.$ (3.27)

Here the subindex ‘osc’ refers to the beginning of the oscillation of $Z$,
which occurs at $t_{\rm osc}\simeq\frac{3}{2}H_{\rm osc}\simeq m_{Z}$.
Assuming, as we did for the inflaton, a quadratic minimum for the potential of
$Z$, we can write $\rho_{\rm osc}\simeq\frac{1}{2}m_{Z}^{2}Z_{0}^{2}$, where
$Z_{0}$ denotes the value of $Z$ at the initial misalignment. Straightforward
substitution gives then

$f_{\chi}(p,t)\;\simeq\;\frac{12\pi^{2}{\rm
Br}_{\chi}}{g_{\chi}(\Gamma_{Z}t)}\left(\frac{Z_{0}\Gamma_{Z}}{m_{Z}^{2}}\right)^{2}\left(\frac{m_{Z}}{2p}\right)^{3/2}e^{-(\Gamma_{Z}t)(2p/m_{Z})^{3/2}+(\Gamma_{Z}/m_{Z})}\theta(m_{Z}/2-p)\,.$
(3.28)

In analogy to the inflaton case, we estimate the decoupling time to be
$t=\Gamma_{Z}^{-1}$. The effect of any subsequent production is to populate
the exponential tail of the distribution. Hence, in what follows we evaluate
the distribution at this decoupling time, and disregard the effect of the
Heaviside function. Moreover, we will always work in the limit when
$\Gamma_{Z}\ll m_{Z}$, as is the case even for stabilized moduli.

To evolve the distribution at later times we make use of the decoupled-regime
solution (A.6). Note that in order to apply it we need to account for the
redshift that occurs from the decay of $Z$ to the end of reheating, and the
subsequent redshift from the end of reheating to present times. Since

$\frac{a(t)}{a_{\rm dec}}\;=\;\frac{a(t)/a_{0}}{a_{\rm
dec}/a_{0}}\;=\;\frac{a(t)/a_{0}}{(a_{\rm dec}/a_{\rm reh})(a_{\rm
reh}/a_{0})}\;\simeq\;\frac{a(t)}{a_{0}}\left(\frac{g_{*s}^{\rm
reh}}{g_{*s}^{0}}\right)^{1/3}\left(\frac{T_{\rm
reh}}{T_{0}}\right)\left(\frac{\Gamma_{Z}}{\Gamma_{\phi}}\right)^{2/3}\,,$
(3.29)

we can finally write, at late times,

$f_{\chi}(p,t)\,\mathop{}\\!\mathrm{d}^{3}\boldsymbol{p}\;\simeq\;\frac{16\pi^{2}{\rm
Br}_{\chi}}{g_{\chi}}\left(\frac{Z_{0}}{m_{Z}}\right)^{2}\left(\frac{\Gamma_{Z}}{m_{Z}}\right)^{2}\left(\frac{a(t)}{a_{0}}\right)^{3}T_{\star,a}^{3}\bar{f}_{{\rm
M},a}(q)\,\mathop{}\\!\mathrm{d}^{3}\boldsymbol{q}\,,$ (3.30)

where

$T_{\star,a}\;=\;\frac{m_{Z}}{2T_{\rm reh}}\left(\frac{g_{*s}^{0}}{g_{*s}^{\rm
reh}}\right)^{1/3}\left(\frac{\Gamma_{\phi}}{\Gamma_{Z}}\right)^{2/3}T_{0}\,,$
(3.31)

and

$\bar{f}_{{\rm M},a}(q)\;=\;\frac{3}{4}q^{-3/2}e^{-q^{3/2}}\,.$ (3.32)

Figure 7: The rescaled distribution functions $\bar{f}_{{\rm M},i}$,
$i=\\{a,b,c\\}$, defined in (3.32), (3.40) and (3.44), as a functions of the
rescaled momentum $q$, for DM produced from modulus decay. For the case
$m_{Z}>\Gamma_{\phi}>\Gamma_{Z}$, the reheating-radiation domination
transition scale $q_{Z}$ has been chosen here to be $q_{Z}=1/4$, and both the
analytical approximation and a numerical solution are shown.

Fig. 7 shows the form of this rescaled distribution (blue, solid curve). The
low momentum power-law dependence and the exponential tail are evident.
Clearly, this distribution is of the form (2.27) with $\alpha=-\gamma=3/2$ and
$\beta=1$.

Case b: Oscillation during reheating, decay after reheating
($m_{Z}>\Gamma_{\phi}>\Gamma_{Z}$)

Let us now consider the case for which $Z$ starts oscillating during
reheating, and its decay is not completed until the subsequent radiation
domination. It is crucial to notice that when this occurs there are two
possible solutions for Eq. (3.6),

$\hat{t}\;\simeq\;\begin{cases}t\left(\dfrac{2p}{m_{Z}}\right)^{2}\,,&p>p_{\rm
reh}\\\\[10.0pt] t_{\rm reh}\left(\dfrac{t}{t_{\rm
reh}}\right)^{3/4}\left(\dfrac{2p}{m_{Z}}\right)^{3/2}\,,&p<p_{\rm
reh}\,,\end{cases}\,,\qquad p_{\rm reh}\equiv\frac{m_{Z}}{2}\left(\frac{t_{\rm
reh}}{t}\right)^{1/2}\,.$ (3.33)

Here we have assumed for simplicity a sharp transition from matter to
radiation domination at $t_{\rm reh}$, with $a\propto t^{2/3}$ in the former
case and $a\propto t^{1/2}$ in the later case. This approximation necessarily
leads to a discontinuity in the Hubble parameter, which will translate into a
discontinuity in the distribution function $f_{\chi}$. This is nothing but an
artifact of our approximations, and it has minimal phenomenological
consequences as we will show below.

For $p>p_{\rm reh}$ we have $\hat{t}>t_{\rm reh}$. In this case we write the
number density of $Z$ as follows,

$\displaystyle n_{Z}(\hat{t})\;$ $\displaystyle\simeq\;\frac{\rho_{\rm
osc}}{m_{Z}}\left(\frac{a_{\rm reh}}{a_{\rm
osc}}\right)^{-3}\left(\frac{a(\hat{t})}{a_{\rm
reh}}\right)^{-3}e^{-\Gamma_{Z}(\hat{t}-t_{\rm osc})}$
$\displaystyle\simeq\;\frac{1}{2}m_{Z}Z_{0}^{2}\left(\Gamma_{\phi}t\right)^{-3/2}\left(\frac{\Gamma_{\phi}}{m_{Z}}\right)^{2}\left(\frac{m_{Z}}{2p}\right)^{3}e^{-(\Gamma_{Z}t)(2p/m_{Z})^{2}}\,,$
(3.34)

and

$H(\hat{t})\;\simeq\;\frac{1}{2t}\left(\frac{m_{Z}}{2p}\right)^{2}\,.$ (3.35)

On the other hand, if $p<p_{\rm reh}$, $\hat{t}<t_{\rm reh}$. Therefore,

$\displaystyle n_{Z}(\hat{t})\;$ $\displaystyle\simeq\;\frac{\rho_{\rm
osc}}{m_{Z}}\left(\frac{a(t)}{a_{\rm
osc}}\right)^{-3}e^{-\Gamma_{Z}(\hat{t}-t_{\rm osc})}$
$\displaystyle\simeq\;\frac{1}{2}m_{Z}Z_{0}^{2}\left(\Gamma_{\phi}t\right)^{-3/2}\left(\frac{\Gamma_{\phi}}{m_{Z}}\right)^{2}\left(\frac{m_{Z}}{2p}\right)^{3}\exp\left[-\left(\Gamma_{\phi}t\right)^{3/4}\left(\frac{\Gamma_{Z}}{\Gamma_{\phi}}\right)\left(\frac{2p}{m_{Z}}\right)^{3/2}\right]\,,$
(3.36)

and

$H(\hat{t})\;\simeq\;\frac{2}{3}\Gamma_{\phi}\left(\Gamma_{\phi}t\right)^{-3/4}\left(\frac{m_{Z}}{2p}\right)^{3/2}\,.$
(3.37)

By substituting into (3.5) and evaluating at $t_{\rm dec}=\Gamma_{Z}^{-1}$ we
obtain the distribution at decoupling. Moreover, noting that in this case the
redshift occurs in the absence of intermediate entropy production, we can
finally write the form of the distribution at late times in the following
simplified way,

$f_{\chi}(p,t)\,\mathop{}\\!\mathrm{d}^{3}\boldsymbol{p}\;\simeq\;\frac{16\pi^{2}{\rm
Br}_{\chi}}{g_{\chi}}\left(\frac{Z_{0}}{m_{Z}}\right)^{2}\left(\frac{\Gamma_{\phi}}{m_{Z}}\right)^{2}\left(\frac{\Gamma_{Z}}{\Gamma_{\phi}}\right)^{3/2}\left(\frac{a(t)}{a_{0}}\right)^{3}T_{\star,b}^{3}\bar{f}_{{\rm
M},b}(q)\,\mathop{}\\!\mathrm{d}^{3}\boldsymbol{q}\,,$ (3.38)

with

$T_{\star,b}\;=\;\frac{m_{Z}}{2T_{\rm dec}}\left(\frac{g_{*s}^{0}}{g_{*s}^{\rm
reh}}\right)^{1/3}T_{0}\,,$ (3.39)

where $T_{\rm dec}=(45/(2\pi^{2}g_{*s}^{\rm
dec}))^{1/4}(\Gamma_{Z}M_{P})^{1/2}$ denotes the background temperature at the
moment of decay, and

$\bar{f}_{{\rm
M},b}(q)\;=\;\begin{cases}q^{-1}e^{-q^{2}}\,,&q>q_{Z}\\\\[10.0pt]
\dfrac{3}{4}q_{Z}^{1/2}q^{-3/2}e^{-q_{Z}^{1/2}q^{3/2}}\,,&q<q_{Z}\end{cases}\,,\qquad
q_{Z}\;\equiv\;\left(\frac{\Gamma_{Z}}{\Gamma_{\phi}}\right)^{1/2}\,.$ (3.40)

This rescaled distribution is shown in Fig. 7 for $q_{Z}=1/4$. The analytical
expression (3.40) is shown as the light green dot-dashed curve. It shows the
different scaling with $q$ for $q>q_{Z}$ and $q<q_{Z}$, with a jump at
$q=q_{Z}$. As we mention above, this discontinuity is an artifact of our
approximations, demonstrated by the dark green, dotted curve in this same
figure, which shows the fully numerical solution, which interpolates smoothly
between the two regimes. Note that for $q_{Z}\sim 1$ the fitting function
(2.27) fails to accurately describe the distribution. Nevertheless, for
$q_{Z}\ll 1$, it accurately describes the DM phase space distribution for any
$q\sim\mathcal{O}(1)$, with $\alpha=-\beta=-1$ and $\gamma=2$.

Case c: Oscillation and decay during radiation domination
($\Gamma_{\phi}>m_{Z}>\Gamma_{Z}$)

For the last case we assume that the beginning of the oscillation of $Z$ is
delayed beyond the end of reheating, due to a rapidly decaying inflaton, a
relatively light $Z$, or a combination of both. The absence of a matter-
radiation crossover during oscillations, and of an intermediate entropy
production regime, make this analysis straightforward. The solution of (3.6)
is simply given by $\hat{t}=t(2p/m_{Z})$, and from it we obtain the following
expressions for the number density in $Z$,

$n_{Z}(\hat{t})\;\simeq\;\frac{1}{2}m_{Z}Z_{0}^{2}(m_{Z}t)^{3/2}\left(\frac{m_{Z}}{2p}\right)^{2}e^{-(\Gamma_{Z}t)(2p/m_{Z})^{2}}\,,$
(3.41)

and the Hubble parameter,

$H(\hat{t})\;=\;\frac{1}{2t}\left(\frac{m_{Z}}{2p}\right)^{2}\,.$ (3.42)

Substitution into (3.5) and (A.6)

$f_{\chi}(p,t)\,\mathop{}\\!\mathrm{d}^{3}\boldsymbol{p}\;\simeq\;\frac{16\pi^{2}{\rm
Br}_{\chi}}{g_{\chi}}\left(\frac{Z_{0}}{m_{Z}}\right)^{2}\left(\frac{\Gamma_{Z}}{m_{Z}}\right)^{3/2}\left(\frac{a(t)}{a_{0}}\right)^{3}T_{\star,c}^{3}\bar{f}_{{\rm
M},c}(q)\,\mathop{}\\!\mathrm{d}^{3}\boldsymbol{q}\,,$ (3.43)

with $T_{\star,c}=T_{\star,b}$ and

$\bar{f}_{{\rm M},c}(q)\;=\;q^{-1}e^{-q^{2}}\,.$ (3.44)

The resulting distribution is trivially of the form (2.27), and is shown in
Fig. 7 as the red, dashed curve.

####  3.2.2 Power spectrum and Ly-$\alpha$ constraints

The analytical determination of the phase space distributions in all cases
allows us to map the WDM Ly-$\alpha$ constraints to the production of DM from
moduli decay. The main hurdle consists in the evaluation of the second moment
of the distribution in the case when the oscillation and the decay of $Z$
occur in different epochs,

$\langle
q^{2}\rangle\;=\;\begin{cases}\Gamma(7/3)\,,&m_{Z}>\Gamma_{Z}>\Gamma_{\phi}\\\\[10.0pt]
e^{-q_{Z}^{2}}(1+q_{Z}^{2})-q_{Z}^{4}E_{-4/3}(q_{Z}^{2})+\dfrac{\Gamma(7/3)}{q_{Z}^{2/3}}\,,&m_{Z}>\Gamma_{\phi}>\Gamma_{Z}\,,\\\\[10.0pt]
1\,,&\Gamma_{\phi}>m_{Z}>\Gamma_{Z}\,.\end{cases}$ (3.45)

Here $E_{n}(x)$ denotes the exponential integral function. Nevertheless, we
find the following to be a good approximation,

$m_{\rm DM}\;\gtrsim\;3.78\,{\rm keV}\,\left(\frac{m_{\rm WDM}}{3\,{\rm
keV}}\right)^{4/3}\left(\frac{g_{*s}^{0}}{g_{*s}^{\rm
reh}}\right)^{1/3}\frac{m_{Z}}{T_{\rm
dec}}\times\begin{cases}\sqrt{\Gamma(7/3)}q_{Z}^{-4/3}\,,&\Gamma_{Z}>\Gamma_{\phi}\,,\\\\[10.0pt]
1\,,&\Gamma_{\phi}>\Gamma_{Z}\,.\end{cases}$ (3.46)

As expected, the limit on the DM mass is weakened if $Z$ decays during
reheating, relative to $Z$ decay during radiation domination. In this case, DM
is cooled down in two stages: from the redshift from $t_{\rm dec}$ to $t_{\rm
reh}$ and from the subsequent redsift from the end of inflation to the present
epoch.

Fig. 8 shows the transfer function for stabilized modulus decay compared to
WDM with $m_{\rm WDM}=1$ and $3\,{\rm keV}$, for the three cases discussed in
this section. The overlap between NCDM and WDM is good for all shown scales,
albeit a slight shift can be observed for $m_{\rm WDM}=1\,{\rm keV}$. As Fig.
3 shows, the relative difference is in all cases $\lesssim 1\%$. The DM masses
are taken from (3.46), where the modulus mass and decay temperature are in
turn chosen to be $m_{Z}\simeq 3\times 10^{6}\,{\rm GeV}$ and $T_{\rm
dec}\simeq 1\,{\rm GeV}$ for case (a), $m_{Z}\simeq 5\times 10^{7}\,{\rm GeV}$
and $T_{\rm dec}\simeq 800\,{\rm GeV}$ for case (b), and $m_{Z}\simeq 3\times
10^{8}\,{\rm GeV}$ and $T_{\rm dec}\simeq 10^{5}\,{\rm GeV}$ for case (c).
These values are motivated by our discussion of the phenomenology of a
strongly stabilized Polonyi-like modulus, in Sec. 3.2.3. For these choices of
the $Z$ mass and decay temperature both the Lyman-$\alpha$ bound and the
closure fraction bound $\Omega_{\chi}h^{2}\simeq 0.1$ are saturated (see Fig.
9).

Figure 8: Linear transfer function for the scenario where DM is produced by
moduli decay, for the cases (a), (b) and (c) described in Section 3.2.1. The
transfer function for the Warm Dark Matter case is shown in gray and black
dashed lines, with $m_{\text{DM}}=1,\,3$ keV, respectively, for comparison.
The numerical values chosen for $m_{\rm DM}$ are estimated from Eq. (3.46),
with $m_{Z}$ and $T_{\rm dec}$ given by the values that saturate the
Ly-$\alpha$ and abundance constraints for the strongly stabilized Polonyi
scenario discussed in Section 3.2.3, shown as stars in Fig. 9 (see text for
details).

####  3.2.3 Relic density and phenomenology

We now consider the possible phenomenological consequences of the Ly-$\alpha$
bound on $m_{\rm DM}$ found above. We first determine the DM relic abundance
from stabilized moduli decays. Integration of Eqs. (3.30), (3.38) and (3.43)
provides the following expression for the late-time DM number density,

$n_{\chi}(t_{0})\;=\;4{\rm
Br}_{\chi}\left(\frac{Z_{0}}{m_{Z}}\right)^{2}\times\begin{cases}\left(\dfrac{\Gamma_{Z}}{m_{Z}}\right)^{2}T_{\star,a}^{3}\,,&m_{Z}>\Gamma_{Z}>\Gamma_{\phi}\\\\[10.0pt]
\left(\dfrac{\Gamma_{\phi}}{m_{Z}}\right)^{2}\left(\dfrac{\Gamma_{Z}}{\Gamma_{\phi}}\right)^{3/2}T_{\star,b}^{3}\,,&m_{Z}>\Gamma_{\phi}>\Gamma_{Z}\,,\\\\[10.0pt]
\left(\dfrac{\Gamma_{Z}}{m_{Z}}\right)^{3/2}T_{\star,c}^{3}\,,&\Gamma_{\phi}>m_{Z}>\Gamma_{Z}\,.\end{cases}$
(3.47)

As mentioned above, the discontinuity in the phase space distribution for
$\chi$ is not inherited by the number density, justifying our approximations.
We emphasize that our results are valid only if the field $Z$ does not
dominate the energy budget of the Universe at any time. For the first
scenario, decay before reheating, this is ensured for $Z_{0}\ll M_{P}$, since
if $\rho_{\rm osc}<\rho_{\phi}(t_{\rm osc})$ then it will continue being so
until the decay of $Z$. For the other two cases we must ensure that the energy
density in radiation, $\rho_{r}$, is always greater than $\rho_{Z}$. Since the
oscillating modulus redshifts more slowly than the background radiation, it is
sufficient to enforce this condition at $Z$-decay. With $\rho_{Z}(t_{\rm
dec})\simeq\frac{1}{2}m_{Z}^{2}Z_{0}^{2}(a_{\rm osc}/a_{\rm dec})^{3}$ and
$\rho_{r}(t_{\rm dec})\simeq\Gamma_{\phi}^{2}M_{P}^{2}(a_{\rm reh}/a_{\rm
dec})^{4}$, we can evaluate the scale factors explicitly to obtain that

$\frac{\rho_{Z}(t_{\rm dec})}{\rho_{r}(t_{\rm
dec})}\;\simeq\frac{1}{2}\left(\frac{Z_{0}}{M_{P}}\right)^{2}\times\begin{cases}\left(\dfrac{\Gamma_{\phi}}{\Gamma_{Z}}\right)^{1/2}\,,&m_{Z}>\Gamma_{\phi}>\Gamma_{Z}\\\\[10.0pt]
\left(\dfrac{m_{Z}}{\Gamma_{Z}}\right)^{1/2}\,,&\Gamma_{\phi}>m_{Z}>\Gamma_{Z}\,.\end{cases}$
(3.48)

This ratio must be $<1$ if our present analysis is to be valid. Otherwise, a
$Z$-dominated epoch occurs, and the bound (3.26) applies. Note that for a
branching ratio ${\rm Br}_{\chi}=1$, the condition $\rho_{Z}\ll\rho_{r}$ is
necessary to obtain the observed DM abundance. For example, saturating the
Ly-$\alpha$ bound (3.46) one obtains $\Omega_{\chi}h^{2}\sim
270(\rho_{Z}/\rho_{r})$ for cases b and c.

We now consider as a proof-of-concept example a particular realization of
modulus stabilization, corresponding to a strongly stabilized Polonyi
field111111The Polonyi field, if left unstabilized, is an example of a
problematic modulus for BBN that can arise in $\mathcal{N}=1$ supergravity
[125, 126, 127, 128, 129]. This field, responsible for the breaking of
supersymmetry, communicates with the SM through Planck-suppressed
interactions. It is also relatively light: its mass of the order of the
gravitino mass, which in turn is parametrically related to the scale at which
supersymmetry is broken. Moreover, typically its initial misalignment is
$\mathcal{O}(M_{P})$. in $\mathcal{N}=1$ supergravity, stabilized by the non-
minimal addition to the Kähler potential $\Delta
K=-(Z\bar{Z})^{2}/\Lambda_{Z}^{2}$ [107, 108, 109, 110, 111, 118]. For our
purposes it is sufficient to note the following values of the Polonyi modulus
mass, its misalignment, and its decay rate

$m_{Z}\;=\;\sqrt{12}\,m_{3/2}\left(\frac{M_{P}}{\Lambda_{Z}}\right)\,,\qquad
Z_{0}\;=\;\frac{\Lambda_{Z}^{2}}{\sqrt{6}M_{P}}\,,\qquad\Gamma_{Z}\;=\;\frac{3\sqrt{3}m_{3/2}^{3}M_{P}^{3}}{\pi\Lambda_{Z}^{5}}\,.$
(3.49)

Here $m_{3/2}\gtrsim\mathcal{O}(10\,{\rm TeV})$ is the gravitino mass for this
particular case of gravity-mediated supersymmetry breaking, and
$\Lambda_{Z}\ll M_{P}$. Hence, the entropy production problem is averted by
simultaneously increasing the $Z$ mass well above the electroweak scale, by
reducing the misalignment to deep sub-Planckian values, and by enhancing the
decay rate. The dominant decay channel of $Z$ is to two gravitinos, which then
subsequently decay into the lightest neutralino. Although in this example this
decay chain implies that the (rescaled) DM distribution will not be exactly
given by the $\bar{f}_{{\rm M},i}(q)$, the scaling of the Ly-$\alpha$
constraint will be maintained up to $\mathcal{O}(1)$ corrections (see Section
4.2).

Figure 9: Allowed range for $\Lambda_{Z}$ as a function of
$d_{\phi}^{-2/3}m_{3/2}/m_{\phi}$, where
$\Gamma_{\phi}=d_{\phi}^{2}m_{\phi}^{3}/M_{P}^{2}$, for the stabilized modulus
defined by (3.49). Shown are the regions excluded by $Z$-domination (entropy
production) and by the Ly-$\alpha$ constraint, assuming ${\rm Br}_{\chi}=1$.
The allowed parameter space is divided into the regions where $Z$ oscillates
and decays after reheating (left), where it begins oscillations during
reheating, and decays after reheating (middle), and where it oscillates and
decays during reheating (right). The orange curve corresponds to
$\Omega_{\chi}h^{2}=0.1$ for $m_{\rm DM}=100\,{\rm GeV}$. Above it, DM is
overproduced. The stars correspond to the points selected to construct the
transfer functions shown in Fig. 8. Where necessary, the gravitino mass is
chosen $m_{3/2}=10^{-13}M_{P}$. See [118] for further details.

Fig. 9 shows the allowed parameter space for $\Lambda_{Z}$ as a function of
the quantity $d_{\phi}^{-2/3}m_{3/2}/m_{\phi}$, where
$\Gamma_{\phi}=d_{\phi}^{2}m_{\phi}^{3}/M_{P}^{2}$. Here
$d_{\phi}\lesssim\mathcal{O}(10^{-1})$ includes the inflaton-matter (or
radiation) couplings and the phase space factors of the width. This
parametrization is chosen to coincide with that of [118], and is inspired by
the Planck-suppressed decays which are a generic feature of supersymmetric
reheating (see e.g. [130, 119]). As it can be seen, when the stabilization
scale is close to the Planck scale, the modulus ceases to be strongly
stabilized, and it dominates the energy budget of the Universe after
inflation. This is averted for

$\Lambda_{Z}\;\lesssim\;M_{P}\times\begin{cases}1.5\left(\dfrac{m_{3/2}^{3}}{M_{P}^{2}\Gamma_{\phi}}\right)^{1/13}\,,&m_{Z}>\Gamma_{\phi}>\Gamma_{Z}\,,\\\\[10.0pt]
1.4\left(\dfrac{m_{3/2}}{M_{P}}\right)^{1/6}\,,&\Gamma_{\phi}>m_{Z}>\Gamma_{Z}\,.\end{cases}$
(3.50)

In this figure we have also shown the domain restricted by Ly-$\alpha$
observations. We observe that it extends the disallowed region (due to entropy
production) by about an order of magnitude in $\Lambda_{Z}$. Its boundary, and
the orange line for which the observed DM abundance is obtained for $m_{\rm
DM}=100\,{\rm GeV}$, are determined through the following expression,

$\displaystyle\Omega_{\chi}h^{2}\;\simeq\;0.1$
$\displaystyle\left(\frac{106.75}{g_{*s}^{\rm
reh}}\right)^{1/4}\left(\frac{m_{\chi}}{100\,{\rm GeV}}\right)$
$\displaystyle\times\begin{cases}\left(\dfrac{\Lambda_{Z}}{6.2\times
10^{14}\,{\rm
GeV}}\right)^{9/2}\left(\dfrac{10^{-13}M_{P}}{m_{3/2}}\right)^{-1/2}\,,&\Gamma_{\phi}>m_{Z}\,,\\\\[10.0pt]
d_{\phi}\left(\dfrac{\Lambda_{Z}}{2.4\times 10^{15}\,{\rm
GeV}}\right)^{5}\left(\dfrac{m_{\phi}}{3\times 10^{13}\,{\rm
GeV}}\right)^{3/2}\left(\dfrac{10^{-13}M_{P}}{m_{3/2}}\right)\,,&m_{Z}>\Gamma_{\phi}\,.\end{cases}$
(3.51)

In the parameter range shown in the figure, the Ly-$\alpha$ and DM abundance
constraints are simultaneously saturated for

$m_{\rm DM}\;\simeq\;3.5\,{\rm MeV}\,\left(\frac{m_{\rm WDM}}{3\,{\rm
keV}}\right)^{4/3}\times\begin{cases}\left(\dfrac{1.3\times
10^{-5}}{d_{\phi}}\right)^{3/13}\left(\dfrac{10^{-13}M_{P}}{m_{3/2}}\right)^{2/13}\,,&m_{Z}>\Gamma_{\phi}>\Gamma_{Z}\,,\\\\[10.0pt]
1\,,&\Gamma_{\phi}>m_{Z}>\Gamma_{Z}\,,\end{cases}$ (3.52)

assuming ${\rm Br}_{\chi}=1$. For
$d_{\phi}^{-2/3}m_{3/2}/m_{\phi}=\\{10^{-1},10^{-6},10^{-10}\\}$, $m_{\rm
DM}\simeq\\{3.99\,{\rm GeV},71.3\,{\rm MeV},\allowbreak 3.70\,{\rm MeV}\\}$,
c.f. Fig. 8. We finish by noting that for this particular stabilization
scenario, the Ly-$\alpha$ constraint is irrelevant compared to the requirement
that $\Omega_{\chi}h^{2}\simeq 0.1$ assuming electroweak-scale LSP masses.
Nevertheless, the power spectrum bound may be relevant for alternative
constructions in which the modulus mass and the DM mass are independent.

##  4 Freeze-in via decay

In the previous section we considered the production of DM from the decay of
the spatially homogeneous condensate. We now extend our discussion to decays
of particles with distributions populated above the zero-momentum mode.
Specifically, we will determine the phase space distribution and the mass
lower bound for DM produced from the decay of a thermalized relic, and from
the decay of a non-thermalized inflaton decay product. As in all cases, we
will assume that DM interactions are sufficiently suppressed to prevent it
from reaching kinetic and/or chemical equilibrium. For this reason we dub this
scenario freeze-in through decays [6].

###  4.1 Thermal decay

####  4.1.1 DM phase space distribution

Let us first consider the decay of a population of particles in thermal
equilibrium, which decays during radiation domination totally or partially
into DM. For definiteness we will assume again that the unstable particle,
denoted here by $A$, decays to DM, $\chi$, via a two-body channel,
$A\rightarrow\chi+\psi$.

The integration of the corresponding collision term can be performed in
complete analogy to the inflaton decay scenario (see Eq. (3.1.1)). Noting in
particular that, for a two-body decay, the unpolarized amplitude squared is
determined solely by the masses of the initial and final state particles, we
can write

$\displaystyle\mathcal{C}[f_{\chi}(p,t)]\;$
$\displaystyle=\;\frac{|\mathcal{M}|^{2}_{A\rightarrow\chi\psi}}{2p_{0}}\int\frac{\mathop{}\\!\mathrm{d}^{3}{\boldsymbol{k}}}{(2\pi)^{3}2k_{0}}\frac{g_{\psi}\mathop{}\\!\mathrm{d}^{3}{\boldsymbol{p}}_{\psi}}{(2\pi)^{3}2p_{\psi}^{0}}(2\pi)^{4}\delta^{(4)}(k-p-
p_{\psi})f_{A}(k_{0})$ $\displaystyle=\;\frac{{\rm
Br}_{\chi}\Gamma_{A}m_{A}}{p_{0}\sqrt{p_{0}^{2}-m_{\text{DM}}^{2}}}\int_{k_{-}}^{k_{+}}\mathop{}\\!\mathrm{d}k_{0}\,f_{A}(k_{0})\,,$
(4.1)

where

$\displaystyle 2m_{\text{DM}}^{2}k_{\pm}\;=\;p_{0}$
$\displaystyle(m_{A}^{2}+m_{\text{DM}}^{2}-m_{\psi}^{2})$
$\displaystyle\pm\sqrt{(p_{0}^{2}-m_{\text{DM}}^{2})(m_{A}^{4}+m_{\text{DM}}^{4}+m_{\psi}^{4}-2m_{\text{DM}}^{2}m_{\psi}^{2}-2m_{\text{DM}}^{2}m_{A}^{2}-2m_{\psi}^{2}m_{A}^{2})}\,.$
(4.2)

Note that up to this point no assumptions have been made regarding the form of
$f_{A}$. For our exploration of the decay of a thermalized relic $A$ into DM,
we can assume that $m_{A}\gg m_{\text{DM}},m_{\psi}$, and substitute a thermal
Bose-Einstein (BE) of Fermi-Dirac (FD) form for $f_{A}$,

$f_{A}(k_{0})\;=\;\frac{1}{e^{k_{0}/T}\pm 1}\,.$ (4.3)

Substitution into (4.1) yields the following collision term,

$\displaystyle\mathcal{C}[f_{\chi}(p,t)]\;$ $\displaystyle\simeq\;\frac{{\rm
Br}_{\chi}\Gamma_{A}m_{A}}{p^{2}}\int_{p+\frac{m_{A}^{2}}{4p}}^{\infty}\frac{\mathop{}\\!\mathrm{d}k_{0}}{e^{k_{0}/T}\pm
1}$ $\displaystyle=\;(\pm)\frac{{\rm
Br}_{\chi}\Gamma_{A}m_{A}T}{p^{2}}\ln\left[1\pm\exp\left(-\frac{p}{T}-\frac{m_{A}^{2}}{4pT}\right)\right]\,.$
(4.4)

Disregarding the inverse decay process, and recalling the relation between
time and temperature during radiation domination,

$H\;=\;\left(\frac{\pi^{2}g_{*\rho}(T)}{90}\right)^{1/2}\frac{T^{2}}{M_{P}}\;\simeq\;\frac{1}{2t}\,,$
(4.5)

the solution of the transport equation (3.1) is a straightfoward application
of the freeze-in solution (A.4). After some algebraic manipulation, the DM
phase space distribution can be cast in the following form [131, 132]

$\displaystyle f_{\chi}\left(p,T\right)\,=\,$ $\displaystyle(\pm){\rm
Br}_{\chi}\frac{\Gamma_{A}T^{2}M_{\rm
Pl}}{p^{2}m_{A}^{2}}\left(\dfrac{90}{\pi^{2}}\right)^{1/2}g_{*s}^{2/3}(T)\int_{0}^{m_{A}/T}\mathop{}\\!\mathrm{d}x\,x^{2}\,g_{*s}^{-2/3}(m_{A}/x)\,g_{*\rho}^{-1/2}(m_{A}/x)$
$\displaystyle\times\left(1-\dfrac{1}{3}\dfrac{\mathop{}\\!\mathrm{d}\log
g_{*s}}{\mathop{}\\!\mathrm{d}\log
x}\right)\ln\left[1\pm\exp\left(-\dfrac{p}{T}\left(\dfrac{g_{*s}(m_{A}/x)}{g_{*s}(T)}\right)^{\frac{1}{3}}-\frac{x^{2}T}{4p}\left(\dfrac{g_{*s}(T)}{g_{*s}(m_{A}/x)}\right)^{\frac{1}{3}}\right)\right]\,.$
(4.6)

Such expression is valid up to the decoupling temperature $T>T_{\rm dec}\sim
m_{A}$ below which the dark matter production from the thermal bath is
negligible.

Figure 10: The rescaled distribution function $\bar{f}_{{\rm TD}}$, defined in
(4.7), as a function of the rescaled momentum $q=p/T$, assuming $T\ll m_{A}\ll
T_{0}$. Solid: the numerically computed phase space distributions for a
fermionic (blue) or bosonic (red) decaying thermalized particle. Dashed: the
phenomenological fits (4.8).

A closed form for $f_{\chi}$ for either bosonic or fermionic $A$ is not
available, and (4.6) must be integrated numerically. These distributions are
presented in Fig. 10 in the limit when $T\ll m_{A}\ll T_{\rm reh}$, by
neglecting the temperature evolution of the effective degrees of freedom
during production, in terms of the rescaled distribution

$\bar{f}_{{\rm TD}}(q)\;\equiv\;\sqrt{\frac{g_{*s}^{\rm dec}}{90}}\frac{\pi
m_{A}^{2}}{{\rm Br}_{\chi}\Gamma_{A}M_{P}}f_{\chi}(q)\,.$ (4.7)

Here $q=p/T$, noting that (4.5) can be extended up to recombination, where
$g_{*s}\simeq g_{*s}^{0}$. The continuous red (blue) curve corresponds to a
decaying fermion (boson) $A$. It is worth noting that the difference between
the two curves is relatively small, which suggests that a phenomenological
Maxwell-Boltzmann-like fit could describe these distributions. Indeed, Fig. 10
also shows two dashed curves which correspond to the following fitting
functions,

$\bar{f}_{{\rm
TD}}(q)\;\simeq\;q^{-1/2}e^{-q}\times\begin{cases}3.38\,,&\text{FD}\,,\\\
3.77\,,&\text{BE}\,.\end{cases}$ (4.8)

Save for the fitting factors, the functional form for this expression may
trivially be obtained from (4.6) in the Maxwell-Boltzmann limit, for which
$\ln\left[1\pm\exp\left(-\frac{p}{T}-\frac{x^{2}T}{4p}\right)\right]\rightarrow\pm\exp\left(-\frac{p}{T}-\frac{x^{2}T}{4p}\right)$
[133, 134, 72]. Worth noting is the mapping of the exponential tail from the
thermalized progenitor $A$ to the daughter particles. Nevertheless, the low-
momentum behavior is different, manifesting the lack of thermal equilibrium in
the $\chi$ sector. This distribution is of the form (2.27), with $\gamma=1$.

####  4.1.2 Power spectrum and Ly-$\alpha$ constraints

The fact that the phase space distribution of $\chi$ is quasi-thermal suggests
that the power spectrum should match the one of WDM. Fig. 11 attests the
reliability of this matching. The leftmost set of curves shows the transfer
functions for the thermal decay cases with BE or FD initial states with masses
determined by Eq. (2.25), which in this case corresponds to the following
rescaled bound,

Figure 11: Linear transfer function for the scenario where DM is produced by
decay of a thermalized particle (denoted by $A$ in the main text), assuming a
Fermi-Dirac (FD), Bose-Einstein (BE) or a fitting phase space distribution as
described in Section 3.1, by taking the mass estimated in Eq. (4.9). The
transfer function for the WDM case is shown for comparison in a black dashed
line. Also depicted here are the transfer functions for FD, BE and fitted
phase space distribution (4.8) with identical masses
$m_{\text{DM}}=\leavevmode\nobreak\ 3\ \text{keV}$.

$m_{\text{DM}}\,\gtrsim\,\,\left(\dfrac{m_{\text{WDM}}}{3\leavevmode\nobreak\
\text{keV}}\right)^{4/3}\left(\dfrac{106.75}{g_{*s}(m_{A})}\right)^{1/3}\times\begin{cases}7.51\leavevmode\nobreak\
\text{keV}\,,\quad&{\rm FD}\,,\\\ 7.32\leavevmode\nobreak\
\text{keV}\,,\quad&{\rm BE}\,,\\\ 7.43\leavevmode\nobreak\
\text{keV}\,,\quad&{\rm Fit}\,.\end{cases}$ (4.9)

Here ‘fit’ stands for both the FD and BE approximations (4.8), which differ
only by a $q$-independent numerical factor. The overlap of these transfer
functions with the WDM result is evident in the whole range of scales shown in
the figure, the relative deviation being $\simeq 1\%$ at $k_{1/2}^{\rm WDM}$
(see Fig. 3). In Fig. 11 we also show the form of $\mathcal{T}(k)$ if we
consider a smaller DM mass, and ignore the difference in statistics. In this
case, all three curves shift to the left, as expected, but the difference
between them remains small. As mentioned earlier, this is the result of the
relatively minimal dependence of $f_{\chi}$ on the spin of the decaying
particle $A$.

####  4.1.3 Relic density and phenomenology

In addition to the power spectrum constraint on the mass discussed above, one
must address the limit from the DM abundance which determines the
normalization of the $\chi$ distribution function. Integration of $f_{\chi}$
gives the following expression for the DM number density at late times, $T\ll
T_{\rm reh}$,

$n_{\chi}(T)\;\simeq\;\sqrt{90}\frac{g_{\chi}{\rm
Br}_{\chi}\Gamma_{A}M_{P}T^{3}}{2\pi^{3}m_{A}^{2}}g_{*s}(T)\left(\dfrac{1}{g_{*s}^{\rm
dec}}\right)^{3/2}\times\begin{cases}4.58\,,\quad&{\rm FD}\,,\\\
4.89\,,\quad&{\rm BE}\,.\end{cases}$ (4.10)

Correspondingly,

$\Omega_{\chi}h^{2}\;\simeq\;0.12\,{\rm
Br}_{\chi}\,\left(\dfrac{g_{\chi}}{2}\right)\left(\frac{106.5}{g_{*s}^{\rm
dec}}\right)^{3/2}\left(\frac{m_{\rm DM}}{6\,{\rm
keV}}\right)\left(\frac{\Gamma_{A}}{10^{-14}\,{\rm
GeV}}\right)\left(\frac{1\,{\rm
TeV}}{m_{A}}\right)^{2}\times\begin{cases}1.17\,,\quad&{\rm FD}\,,\\\
1.02\,,\quad&{\rm BE}\,.\end{cases}$ (4.11)

Except for the number of degrees of freedom, which we consistently normalize
to the SM value, the normalizations chosen in the previous equation are
inspired by the decay of thermalized supersymmetric particles into light DM
candidates, such as the Higgsino $\rightarrow$ axino + Higgs production
process in $R$-parity violating DFSZ models [135, 136], for which

$\Gamma(\tilde{H}\rightarrow\tilde{a}+H)\;=\;\frac{1}{8\pi}\left(\frac{\mu}{f_{a}}\right)^{2}\mu\,,$
(4.12)

with the $\mu$-term parameter $\mu\sim 500\,{\rm GeV}$, and the Peccei-Quinn
scale $f_{a}\sim 10^{10}\,{\rm GeV}$. Similarly to the inflaton decay case, a
mass-independent constraint on the branching ratio to DM from the decay of the
thermalized $A$ could be derived. Nevertheless, this bound would not be
universal, as the mass and width of $A$ are model dependent, as opposed to the
inflaton decay case (see Eq. (3.22)).

###  4.2 Non-thermal decay

####  4.2.1 DM phase space distribution

Let us now assume that the particle $A$ whose decay produces the DM interacts
very weakly with the SM and was produced via inflaton decay, but does not
reach thermal equilibrium. Unlike in the previously studied thermal case, this
particle cannot be assumed to be produced abundantly in the thermal plasma
during the decay of the latter, Therefore, in principle the imprint that its
decay leaves on its phase space distribution must be taken into account.

Disregarding the effect of Bose enhancement/Pauli blocking, and the inverse
decay process, the Boltzmann equation satisfied by this non-thermal unstable
relic is given by [137]

$\frac{\partial f_{A}}{\partial t}-Hp\frac{\partial f_{A}}{\partial
p}\;=\;-\frac{m_{A}\Gamma_{A}}{\sqrt{m_{A}^{2}+p^{2}}}f_{A}\,.$ (4.13)

This equation can be exactly solved in the relativistic and non-relativistic
regimes. In both cases the decay of $A$ proceeds exponentially in time. For
this reason we will be content to approximate the evolution of $f_{A}$ as that
of a free-streaming particle until its sudden decay, which occurs at

$t_{\rm
dec}\;\simeq\;\begin{cases}\Gamma_{A}^{-1}\,,&\dfrac{\Gamma_{A}}{H_{A}}\ll
1\,,\\\\[10.0pt] \left(\dfrac{m_{\phi}\langle
q_{A}\rangle}{2m_{A}\Gamma_{A}\Gamma_{\phi}^{1/2}}\right)^{2/3}\,,&\dfrac{\Gamma_{A}}{H_{A}}\gg
1\,.\end{cases}$ (4.14)

Here $H_{A}$ denotes the Hubble parameter at the time when $A$ becomes non-
relativistic. We have estimated the effective lifetime as the inverse of the
mean $f_{A}$ prefactor in the right-hand side of (4.13) [138].

With the previous arguments in mind, for $t<t_{\rm dec}$ we write the
collision term for $\chi$ (4.1) as

$\displaystyle\mathcal{C}[f_{\chi}(p,t)]\;=\;\frac{4\pi^{4}g_{*s}^{\rm
reh}{\rm Br}_{\chi}{\rm Br}_{A}\Gamma_{A}m_{A}}{5g_{A}p^{2}}$
$\displaystyle\left(\frac{T_{\rm
reh}}{m_{\phi}}\right)^{4}\left(\frac{m_{\phi}}{2}\right)\left(\frac{a_{\rm
reh}}{a(t)}\right)$
$\displaystyle\times\int_{\left|\frac{2p}{m_{\phi}}\frac{a(t)}{a_{\rm
reh}}-\frac{m_{A}^{2}}{2pm_{\phi}}\frac{a(t)}{a_{\rm
reh}}\right|}^{\infty}\frac{z\,\mathop{}\\!\mathrm{d}z}{\sqrt{z^{2}+\left(\frac{2m_{A}a(t)}{m_{\phi}a_{\rm
reh}}\right)^{2}}}\,\bar{f}_{\rm R}\left(z\right)\,,$ (4.15)

where the distribution for inflaton decay products $\bar{f}_{\text{R}}$, given
in terms of the 3D momentum magnitude, was defined in (3.16). In this
expression ${\rm Br}_{A}$ stands for the branching ratio of the decay from
inflaton to $A$. Substitution into the general freeze-in solution (A.4) gives

$\displaystyle f_{\chi}(p,t_{\rm dec})\;=\;\;$
$\displaystyle\frac{8\pi^{4}g_{*s}^{\rm reh}{\rm Br}_{\chi}{\rm
Br}_{A}\Gamma_{A}m_{A}}{5g_{A}m_{\phi}}\left(\frac{T_{\rm
reh}}{m_{\phi}}\right)^{4}q_{\rm dec}^{-2}$ $\displaystyle\times\int_{t_{\rm
reh}}^{t_{\rm
dec}}\mathop{}\\!\mathrm{d}t^{\prime}\,\frac{a(t^{\prime})}{a_{\rm
reh}}\int_{\left|q_{\rm dec}-\frac{1}{q_{\rm
dec}}\left(\frac{m_{A}}{m_{\phi}}\frac{a(t^{\prime})}{a_{\rm
reh}}\right)^{2}\right|}^{\infty}\frac{z\,\mathop{}\\!\mathrm{d}z}{\sqrt{z^{2}+\left(\frac{2m_{A}a(t^{\prime})}{m_{\phi}a_{\rm
reh}}\right)^{2}}}\,\bar{f}_{\rm R}(z)\,,$ (4.16)

where $q_{\rm dec}=(2p/m_{\phi})(a_{\rm dec}/a_{\rm reh})$. The ratio

$\displaystyle\frac{m_{A}a(t)}{m_{\phi}a_{\rm
reh}}\;\propto\;\frac{m_{A}}{\langle p\rangle}\,,$ (4.17)

quantifies how relativistic the distribution for $A$ is at a given moment of
time. In particular, we define

$\mathcal{R}\;\equiv\;\frac{m_{A}a_{\rm dec}}{m_{\phi}a_{\rm
reh}}\;=\;\left(\frac{g_{*s}^{\rm reh}}{g_{*s}^{\rm
dec}}\right)^{1/3}\frac{m_{A}T_{\rm reh}}{m_{\phi}T_{\rm
dec}}\;=\;\begin{cases}\left(\dfrac{2H_{A}}{\Gamma_{A}}\right)^{1/2}\gg
1\quad&\text{for}\quad\dfrac{\Gamma_{A}}{H_{A}}\ll 1\,,\\\\[10.0pt]
\left(\langle q_{A}\rangle\dfrac{3H_{A}}{2\Gamma_{A}}\right)^{1/3}\ll
1\quad&\text{for}\quad\dfrac{\Gamma_{A}}{H_{A}}\gg 1\,.\end{cases}\,.$ (4.18)

Extending the solution past $t_{\rm dec}$ we can write

$\displaystyle
f_{\chi}(p,t)\,\mathop{}\\!\mathrm{d}^{3}\boldsymbol{p}\;=\;\frac{24\pi^{3}\sqrt{10g_{*s}^{\rm
reh}}{\rm Br}_{\chi}{\rm
Br}_{A}\Gamma_{A}M_{P}}{5g_{A}m_{A}^{2}}\left(\frac{T_{\rm
reh}}{m_{\phi}}\right)^{2}\mathcal{F}(q,\mathcal{R})\left(\frac{a_{0}}{a(t)}\right)^{3}T_{\star}^{3}\mathop{}\\!\mathrm{d}^{3}\boldsymbol{q}\,,$
(4.19)

where

$\mathcal{F}(q,\mathcal{R})\;=\;q^{-2}\int_{0}^{\mathcal{R}}\mathop{}\\!\mathrm{d}y\,y^{2}\int_{\left|q-\frac{y^{2}}{q}\right|}^{\infty}\frac{z\,\mathop{}\\!\mathrm{d}z}{\sqrt{q^{2}+4y^{2}}}\,\bar{f}_{\rm
R}(z)\;\simeq\;\begin{cases}\bar{f}_{\rm D,NR}(q)\,,&\mathcal{R}\gg
1\,,\\\\[10.0pt] \dfrac{\mathcal{R}^{3}}{3}\bar{f}_{\rm
D,R}(q)\,,&\mathcal{R}\ll 1\,.\end{cases}$ (4.20)

Here $q$ and $T_{\star}$ are the same as in (3.16) and (3.17). The rescaled
distributions $\bar{f}_{\rm D,NR}$ and $\bar{f}_{\rm D,R}$ can be computed by
making use of the fit approximation (3.18) for $\bar{f}_{\rm R}$. We obtain

$\displaystyle\bar{f}_{\rm D,NR}(q)\;$
$\displaystyle\simeq\;0.36\,q^{-1}\left[0.43\,q\,\Gamma\left(\frac{1}{4},0.19\,q^{2}\right)-\Gamma\left(\frac{3}{4},0.19\,q^{2}\right)+2\,\Gamma\left(\frac{3}{4}\right)\right]\theta\big{(}\mathcal{R}-q\big{)}\,,$
(4.21) $\displaystyle\bar{f}_{\rm D,R}(q)\;$
$\displaystyle=\;q^{-2}\int_{q}^{\infty}\mathop{}\\!\mathrm{d}z\,\bar{f}_{\rm
R}(z)\;\simeq\;1.06\,q^{-2}\,\Gamma\left(-\frac{1}{4},0.74\,q^{2}\right)\,,$
(4.22)

where $\Gamma(a,x)$ denotes the upper incomplete gamma function.

Figure 12: The rescaled distribution function $\bar{f}_{{\rm D,NR}}$, defined
in (4.21), as a function of the rescaled momentum $q$ and the order parameter
$\mathcal{R}=(m_{A}/m_{\phi})(\Gamma_{\phi}/\Gamma_{A})^{1/2}$. Solid, black:
numerically computed phase space distribution. Dashed, orange: the fit
approximation (4.21).

The DM phase space distribution corresponding to the decay of a non-
relativistic particle $A$ is shown in Fig. 12 as a function of $q$ and
$\mathcal{R}>1$. The solid black line shows the result of the numerical
integration of (4.20). The distribution grows with an almost linear universal
envelope, independent of the decay rate of $A$, until $q\sim\mathcal{R}$, at
which point the distribution sharply decreases. This non-universality of the
cutoff prevents us from constructing a reasonable fit approximation of the
form (2.27) for generic values of $\mathcal{R}$. In the same figure, the
orange dashed lines show the analytical approximation (4.21), which as can be
seen is equivalent to imposing a hard cutoff at $q=\mathcal{R}$ on the
universal envelope.121212The numerical distribution can be well fitted by
substituting the $\theta$ function in Eq. (4.21) by a logistic function.

Fig. 13 shows the numerically computed relativistic distribution $\bar{f}_{\rm
D,R}$ as the solid black curve, and the analytical approximation given by
(4.22) as the orange, dashed curve. In the same figure a ‘fit’ approximation
of the form (2.27) is also shown. This approximation is obtained by mimicking
the asymptotic behavior of the gamma function at large and small $q$, while
preserving the normalization, and is given by

$\displaystyle\bar{f}_{\rm ND}(q)\;$
$\displaystyle\approx\;2.19q^{-5/2}e^{-0.74q^{2}}\,.$ (4.23)

It is worth noting that in this case the Gaussian tail is of the same form as
that of the parent unstable particle. It is important to emphasize that this
distribution is obtained in the limit $\mathcal{R}\rightarrow 0$, as we
discuss below.

Figure 13: The rescaled distribution function $\bar{f}_{{\rm D,R}}$, defined
in (4.22), as a function of the rescaled momentum $q$. Solid, black:
numerically computed phase space distribution. Dashed, orange: the analytical
approximation (4.22). Dashed-dotted, blue: the fit approximation (4.23).

Fig. 14 shows the form of the function $\mathcal{F}(q,\mathcal{R})$, defined
in Eq. (4.20), for several values of $\mathcal{R}$, ranging from $10^{-2}$ to
10. Here we can appreciate the transition between the relativistic and non-
relativistic decay cases. In all cases the phase space distribution peaks at
$q\simeq\mathcal{R}$, with a positive skew for a relativistic $A$, and a
negative skew for non-relativistic $A$. For $\mathcal{R}<1$ the analytical
approximation (4.22) describes well the exact distribution for
$q\gtrsim\mathcal{R}$. For $\mathcal{R}>1$, the non-relativistic approximation
(4.21) is in turn a good fit for the exact distribution for $q\gtrsim 1/2$.

Figure 14: The function $\mathcal{F}(q,\mathcal{R})$, defined in (4.20), as a
function of the rescaled momentum $q$ and the decay parameter $\mathcal{R}$.
Solid: numerically computed distributions. Dashed: the non-relativistic
analytical approximation (4.21) for $\mathcal{R}\gg 1$. Dashed-dotted: the
relativistic analytical approximation (4.22).

####  4.2.2 Power spectrum and Ly-$\alpha$ constraints

For the distributions that we have derived, we can make use of (2.25) to
determine the rescaling of the bound on the DM mass. For the case of a
relativistic (R) decay we find that

$\displaystyle m_{\text{DM}}\,\gtrsim\,\,$
$\displaystyle\left(\dfrac{m_{\text{WDM}}}{3\leavevmode\nobreak\
\text{keV}}\right)^{4/3}\left(\dfrac{106.75}{g_{*s}^{\rm reh}}\right)^{1/3}$
$\displaystyle\times\left(\dfrac{m_{\phi}}{3\times 10^{13}\leavevmode\nobreak\
\text{GeV}}\right)\left(\dfrac{10^{10}\leavevmode\nobreak\
\text{GeV}}{T_{\text{reh}}}\right)\times\begin{cases}1.23\leavevmode\nobreak\
\text{MeV}\,,\leavevmode\nobreak\ &{\rm Numerical}\leavevmode\nobreak\
{\rm(R)}\,,\\\ 1.26\leavevmode\nobreak\ \text{MeV}\,,\leavevmode\nobreak\
&{\rm Analytical}\leavevmode\nobreak\ {\rm(R)}\,,\\\ 2.19\leavevmode\nobreak\
\text{MeV}\,,\leavevmode\nobreak\ &{\rm Fit}\leavevmode\nobreak\
{\rm(R)}\,,\\\ \end{cases}$ (4.24)

while for the non-relativistic (NR) case,

$\displaystyle m_{\text{DM}}\,\gtrsim\,\,$
$\displaystyle\left(\dfrac{m_{\text{WDM}}}{3\leavevmode\nobreak\
\text{keV}}\right)^{4/3}\left(\dfrac{106.75}{g_{*s}^{\rm reh}}\right)^{1/3}$
$\displaystyle\times\left(\dfrac{m_{\phi}}{3\times 10^{13}\leavevmode\nobreak\
\text{GeV}}\right)\left(\dfrac{10^{10}\leavevmode\nobreak\
\text{GeV}}{T_{\text{reh}}}\right)\left(\dfrac{\mathcal{R}}{6}\right)\times\begin{cases}16.1\leavevmode\nobreak\
\text{MeV}\,,\leavevmode\nobreak\ &{\rm Analytical}\leavevmode\nobreak\
{\rm(NR)}\,,\\\ 16.2\leavevmode\nobreak\ \text{MeV}\,,\leavevmode\nobreak\
&{\rm Numerical}\leavevmode\nobreak\ {\rm(NR)}\,.\end{cases}$ (4.25)

The $m_{\text{DM}}\propto\mathcal{R}$ behavior is only correct for large
$\mathcal{R}\gg 1$ but remains a reasonable approximation for
$\mathcal{R}\sim\mathcal{O}(1-10)$. We note here that the lower bound on the
NCDM mass can be many orders of magnitude larger than the corresponding WDM
bound, and it increases as the reheating temperature is decreased. This is
expected, as in this case the parent particle is produced from inflaton decay
(see Sec. 3.1.2). The difference between the numerical and analytical results
is minimal, consistent with the agreement between both curves in Fig. 12 and
Fig. 13. However, the fit approximation for the relativistic case provides a
relatively poor approximation to the bound, overestimating it by a factor of
$\sim 1.8$.

Figure 15: Linear transfer function for DM produced by the decay of a non-
thermalized non-relativistic particle. We show here the numerical results for
two sets of cosmological parameters: $T_{\rm reh}=10^{12}\,{\rm GeV}$ and
$m_{\rm WDM}=1\,{\rm keV}$, and $T_{\rm reh}=10^{10}\,{\rm GeV}$ and $m_{\rm
WDM}=3\,{\rm keV}$, making use of the rescaled bound (4.25). For comparison we
also show the transfer function for the corresponding WDM cases. Figure 16:
Linear transfer function for DM produced by the decay of a relativistic non-
thermalized particle. We show here the numerical, analytical and fit
approximations discussed in the text, for two sets of cosmological parameters:
$T_{\rm reh}=10^{12}\,{\rm GeV}$ and $m_{\rm WDM}=1\,{\rm keV}$, and $T_{\rm
reh}=10^{10}\,{\rm GeV}$ and $m_{\rm WDM}=3\,{\rm keV}$, making use of the
rescaled bound (4.24). For comparison we also show the transfer function for
the corresponding WDM cases.

Figs. 15 and 16 show the results of the numerical evaluation of the transfer
functions with CLASS [44, 45], and their comparison with the WDM
case.131313For the relativistic case, the small disagreement between the
numerical transfer function with values from Eq. (4.24) and the corresponding
WDM spectrum is also attributed to the sharp drop of the phase space
distribution for $q<\mathcal{R}\ll 1$, akin to a low-momentum cutoff. Such a
cutoff results in a loss of numerical precision if reasonable computation
times are required. For the two sets of curves shown in each figure, we use
the rescaled Ly-$\alpha$ bound (4.24) or (4.25). For the leftmost set we take
$T_{\rm reh}=10^{12}\,{\rm GeV}$ and $m_{\rm WDM}=1\,{\rm keV}$, while for the
rightmost set we consider $T_{\rm reh}=10^{10}\,{\rm GeV}$ and $m_{\rm
WDM}=3\,{\rm keV}$. For the decay of a non-relativistic particle, a comparison
is made between the three different choices for $\mathcal{R}=2,6$ and $10$.
Note the overlap between the three curves, with a relative difference of $\sim
1\%$ (see Fig. 3, where the relative difference is plotted as a function of
$k$ for $\mathcal{R}=6$). For the decay of a relativistic particle, the
agreement between the numerical and analytical results can be immediately
appreciated, as well as the difference between these and the result of using
the fit approximation (4.23) for the DM distribution.141414The analytical
expression of Eq. (4.21) is not represented in Fig. 15 as the sharp
$\theta$-function cannot be handled properly with CLASS as it requires the
distribution function to smoothly decrease at large $q$. Even more evident
though is the difference of the NCDM transfer functions with respect to the
one for WDM, of around $10\%$ at $k_{1/2}^{\rm WDM}$, c.f. Fig. 3. For the
relativistic case, the distribution $f_{\chi}$ has a very non thermal shape,
monotonically decreasing with $p$, resulting in a power spectrum that,
although not too dissimilar from the WDM case, exhibits in the figure an
appreciable difference from it.

####  4.2.3 Relic density and phenomenology

The present relic abundance of DM is obtained from integration of (4.19). To
do this we make use of the (numerical) result

$\int_{0}^{\infty}\mathop{}\\!\mathrm{d}q\,q^{2}\bar{f}_{\rm
D,NR}(q)\;\simeq\;0.4\,\mathcal{R}^{2}\,.$ (4.26)

At $t\gg t_{\rm dec}$ the number density has the form

$n_{\chi}(t)\;\simeq\;g_{*s}^{0}{\rm Br}_{\chi}{\rm
Br}_{A}\left(\frac{g_{\chi}}{g_{A}}\right)\left(\frac{T_{\rm
reh}}{m_{\phi}}\right)\left(\frac{a_{0}}{a(t)}\right)^{3}T_{0}^{3}\times\begin{cases}\left(\dfrac{g_{*s}^{\rm
reh}}{g_{*s}^{\rm dec}}\right)^{1/6}\,,&\mathcal{R}\gg 1\,,\\\\[10.0pt]
\left(\dfrac{g_{*s}^{\rm reh}}{g_{*s}^{\rm
dec}}\right)^{1/4}\,,&\mathcal{R}\ll 1\,.\end{cases}$ (4.27)

Note that both expressions agree up to a different power of the number of
relativistic degrees of freedom. This agreement is to be expected, as the
total number of the decaying particle $A$ and its decay product must be a
constant in a comoving volume. Considering for definiteness the case of a
relativistic decaying particle, we determine that the present abundance is
given by

$\displaystyle\Omega_{\chi}h^{2}\;\simeq\;0.1\,{\rm Br}_{\chi}$
$\displaystyle\left(\frac{{\rm Br}_{A}}{5.5\times
10^{-4}}\right)\left(\frac{g_{\chi}}{g_{A}}\right)\left(\dfrac{g_{*s}^{\rm
reh}}{g_{*s}^{\rm dec}}\right)^{1/4}\left(\frac{m_{\rm DM}}{1\,{\rm
MeV}}\right)\left(\frac{T_{\rm reh}}{10^{10}\,{\rm
GeV}}\right)\left(\frac{3\times 10^{13}\,{\rm GeV}}{m_{\phi}}\right)\,.$
(4.28)

As expected, $\Omega_{\chi}$ is independent of the properties of $A$, and
corresponds simply to a re-scaling by degrees of freedom of the inflaton decay
result (3.21).

Given this result, a universal lower bound on ${\rm Br}_{\chi}$ can be
obtained, in full analogy with the inflaton decay scenario. Let us discuss it
in the context of a specific model. Consider the decay chain inflaton
$\rightarrow$ gravitino $\rightarrow$ LSP (lightest supersymmetric particle),
which is generically present in supersymetric models of inflation with
supersymmetry breaking mediated gravitationally [139, 140, 141, 142,
143].151515In typical gauge-mediation scenarios, the gravitino can be very
light, $m_{3/2}\sim{\rm keV}$ and is produced through thermal freeze-out, thus
being an example of WDM [144, 145, 146, 147, 39]. Assuming a minimal
supersymmetric extension of the Standard Model (MSSM), the decay rate of the
spin-3/2 gravitino is [148]

$\Gamma_{3/2}\;=\;\frac{193}{384\pi}\frac{m_{3/2}^{3}}{M_{P}^{2}}\,.$ (4.29)

Generically, ${\rm Br}_{\rm LSP}=\mathcal{O}(1)$. Substitution into (4.24) and
(4.28) leads to the following absolute constraints on the branching ratio of
the decay of the inflaton into gravitinos, independent of the DM mass: For
non-relativistic decaying particles, $T_{\rm reh}\gg 10^{5}\,{\rm
GeV}(m_{3/2}/10\,{\rm TeV})^{1/2}$ and

${\rm Br}_{3/2}\;\lesssim\;1.3\times 10^{-8}\left(\dfrac{3\,{\rm keV}}{m_{\rm
WDM}}\right)^{4/3}\left(\frac{m_{\phi}}{3\times 10^{13}\,{\rm
GeV}}\right)\left(\frac{10^{10}\,{\rm GeV}}{T_{\rm
reh}}\right)\left(\frac{m_{3/2}}{10\,{\rm TeV}}\right)^{1/2}\,,$ (4.30)

while for relativistic decaying ones, $T_{\rm reh}\ll 10^{5}\,{\rm
GeV}(m_{3/2}/10\,{\rm TeV})^{1/2}$ and

${\rm Br}_{3/2}\;\lesssim\;1.2\times 10^{-3}\left(\dfrac{3\,{\rm keV}}{m_{\rm
WDM}}\right)^{4/3}\,.$ (4.31)

In this (MSSM) scenario, the excluded DM masses span a phenomenologically
interesting region in the parameter space of the model, as shown in Fig. 17.
The exclusion region corresponds to

$m_{\rm LSP}\;\lesssim\;\begin{cases}86\,{\rm GeV}\left(\dfrac{m_{\rm
WDM}}{3\,{\rm keV}}\right)^{4/3}\left(\dfrac{10\,{\rm
TeV}}{m_{3/2}}\right)^{1/2}\,,&T_{\rm reh}\gg 10^{5}\,{\rm
GeV}\left(\dfrac{m_{3/2}}{10\,{\rm TeV}}\right)^{1/2}\,,\\\\[10.0pt] 95\,{\rm
GeV}\left(\dfrac{m_{\rm WDM}}{3\,{\rm
keV}}\right)^{4/3}\left(\dfrac{10^{5}\,{\rm GeV}}{T_{\rm
reh}}\right)\,,&T_{\rm reh}\ll 10^{5}\,{\rm GeV}\left(\dfrac{m_{3/2}}{10\,{\rm
TeV}}\right)^{1/2}\,.\end{cases}$ (4.32)

These bounds are or the order of the electroweak scale, and are comparable to
collider and direct detection limits [149, 150]. Note that for a model-fixed
LSP mass, the Lyman-$\alpha$ constraint puts a bound on the inflaton-matter
couplings. For $m_{\rm LSP}\gtrsim 100\,{\rm GeV}$, $T_{\rm reh}\lesssim
100\,{\rm TeV}$ are excluded.

A straightforward computation shows that independently of the mean momentum of
the decaying gravitino, the decay occurs at temperatures at which the LSP can
be safely assumed to be decoupled from the thermal plasma, and hence preserves
its non-equilibrium phase space distribution.

Figure 17: Ly-$\alpha$ constraint on the LSP mass, as a function of
temperature, in the case of production through the decay chain inflaton
$\rightarrow$ gravitino $\rightarrow$ LSP. For $T_{\rm reh}>10^{5}\,{\rm
GeV}(m_{3/2}/10\,{\rm TeV})^{1/2}$ the decay of the gravitino occurs when it
is non-relativistic. For $T_{\rm reh}<10^{5}\,{\rm GeV}(m_{3/2}/10\,{\rm
TeV})^{1/2}$, the gravitino is relativistic at the moment of decay.

##  5 Ultraviolet freeze-in via scatterings

In this section we consider the production of light DM from scatterings in the
primordial plasma. We will restrict ourselves to $2\rightarrow 2$ processes,
for which the integrated effective cross section is assumed to be of the form

$\sigma(s)\;=\;\frac{s^{\frac{n}{2}}}{{}^{n+2}}\,,$ (5.1)

where $n$ is an integer and $s$ is the Mandelstam variable, related at high
energies with the center of mass energy $E$ by $\sqrt{s}=E$. Although for non-
negative $n$ this cross section naively violates unitarity at high energies,
we assume that it merely corresponds to the low-energy effective description
of a UV-complete model. The energy scale $\Lambda$ can be thought to be
parametrically related to the mass of a heavy mediator. The suppression by
$\Lambda$ guarantees that the primordial abundance is determined by forward
processes (plasma $\rightarrow$ DM) rather than by annihilations. Therefore,
Pauli-blocking/Bose-enhancement for $\chi$ can be safely disregarded, and in
the absence of other interaction channels, $\chi$ never reaches thermal
equilibrium with the plasma. Thus, freeze-in is realized [6, 151].

Assuming no post-reheating entropy production (that is, a standard thermal
history), particle production is dominated by temperatures $T\geq T_{\rm reh}$
if $n>-1$. This is referred to as ultraviolet (UV) freeze-in [152, 153, 49].
Moreover, for $n>-1$, and for a sufficiently large reheating temperature, we
can safely assume that both the parent scatterers and the produced DM are
ultrarelativistic at the time of production,161616This justifies disregarding
any dependence on thresholds. For $n\leq-1$, the masses of the scatterers play
an important role to determine the lower bound on $m_{\rm DM}$ [71, 72]. if
the former are in thermal equilibrium. If the parent scatterers are not in
equilibrium at production time, the condition that their masses are $\ll
m_{\phi}$ suffices. Here we will consider both scenarios.

In order to evaluate the necessary collision terms for thermal and non-thermal
production, we need to make assumptions regarding the form of the scattering
amplitude. Its dependence on the angles (or Mandelstam variables $s,t,u$)
involved in the scattering varies between different microscopic descriptions
of the process. We will assume that for the scattering process
$A(k)+B(\tilde{k})\rightarrow\chi(p)+\psi(\tilde{p})$, the mean, unpolarized
squared scattering amplitude can be parametrized in the following way,

$|\mathcal{M}|^{2}\;=\;16\pi\frac{s^{\frac{n}{2}+1}}{\Lambda^{n+2}}\,.$ (5.2)

Integration with respect to the two-particle phase space recovers (5.1). For a
different combination of $s,t,u$, our results will generically only differ by
numerical factors, which can be absorbed into the value of
$\Lambda$.171717Exceptions include those cases in which finite-temperature in-
medium effects are necessary to regulate infrared divergences, which arise
from the exchange of massless mediators. Thermal axion production and
gravitino production in low-scale supersymmetry are included in these cases
[154, 155, 156, 157, 158, 159].

Under the freeze-in assumption, and with the square amplitude given by (5.2),
the collision term for the production of $\chi$ can be written as follows,

$\mathcal{C}[f_{\chi}]\;=\;\frac{16\pi
g_{A}g_{B}g_{\psi}}{\Lambda^{n+2}2p_{0}}\int\frac{\mathop{}\\!\mathrm{d}^{3}\tilde{\boldsymbol{p}}}{2(2\pi)^{3}\tilde{p}_{0}}\frac{\mathop{}\\!\mathrm{d}^{3}\boldsymbol{k}}{2(2\pi)^{3}k_{0}}\frac{\mathop{}\\!\mathrm{d}^{3}\tilde{\boldsymbol{k}}}{2(2\pi)^{3}\tilde{k}_{0}}(2\pi)^{4}\delta^{(4)}(p+\tilde{p}-k-\tilde{k})s^{\frac{n}{2}+1}f_{A}(k_{0})f_{B}(\tilde{k}_{0})\,.$
(5.3)

The integration of this collision term for arbitrary $f_{A,B}$ can be easily
done following the steps of [155, 160]. We detail these steps in Appendix A.2.
As a result, we obtain Eq. (A.13), which will be the starting point of our
discussion of thermal and non-thermal UV freeze-in.

###  5.1 Thermal freeze-in

We begin by applying the general solution (A.13) to the production of DM from
thermalized scatterers, i.e. with Fermi-Dirac or Bose-Einstein distributions.
As stated earlier in this section, we focus on UV freeze-in, for which the
bulk of the DM relic abundance is produced during reheating. Thermal
production during reheating is the dominant production channel in the absence
of significant direct inflaton $\rightarrow$ DM decays for $-1<n\leq 2$ in
(5.1). Moreover, for higher $n$, thermal production can dominate over non-
thermal effects if the parent scatterers that couple to the dark sector are