###### Abstract

Why does dark matter (DM) live longer than the age of the Universe? Here we
study a novel sub-eV scalar DM candidate whose stability is due to the Pauli
exclusion of its fermionic decay products. We analyze the stability of the DM
condensate against decays, scatterings (i.e., evaporation), and parametric
resonance, delineating the viable parameter regions in which DM is
cosmologically stable. In a minimal scenario in which the scalar DM decays to
a pair of new exotic fermions, we find that scattering can populate an
interacting thermal dark sector component to energies far above the DM mass.
This self-interacting dark radiation may potentially alleviate the Hubble
tensions. Furthermore, our scenario can be probed through precise measurements
of the halo mass function or the masses of dwarf spheroidal galaxies since
scattering prevents the DM from becoming too dense. On the other hand, if the
lightest neutrino stabilizes the DM, the cosmic neutrino background (C$\nu$B)
can be significantly altered from the $\Lambda$CDM prediction and thus be
probed in the future by C$\nu$B detection experiments.

TU-1168, PITT-PACC-2403

Cosmic Stability of Dark Matter
from Pauli Blocking

Brian Batell1 and Wen Yin2,3,4,5

1 Pittsburgh Particle Physics, Astrophysics, and Cosmology Center, Department
of Physics and Astronomy, University of Pittsburgh, Pittsburgh, USA
2 Department of Physics, Tohoku University, Miyagi 980-8578, Japan
3 Department of Physics, Tokyo Metropolitan University, Tokyo, 192-0397, Japan
4 Department of Physics, Faculty of Science, The University of Tokyo,
Bunkyo-ku, Tokyo 113-0033, Japan
5 Department of Physics, KAIST, Daejeon, Korea

##  1 Introduction

The stability of dark matter (DM) on cosmic timescales provides an important
clue to its fundamental nature. For example, it may suggest DM is endowed with
a conserved charge, has a small mass, and/or has suppressed interactions,
leading to a small or vanishing decay rate. Understanding the possible reasons
behind DM stability is extremely important because they may suggest novel
correlated observational signatures, as is well illustrated by the cases of DM
stabilized by symmetry (e.g., WIMP DM) or by small mass/suppressed
interactions (e.g., ALP DM).

In this paper, we explore a novel explanation for the stability of DM in our
Universe that exploits the Pauli exclusion principle. In a simple model of
sub-eV scalar DM with sizable Yukawa couplings to even lighter or massless
exotic fermions, the DM proper lifetime in vacuum is much shorter than the age
of the Universe. However, the fermionic states accessible in the decay are
rapidly filled, thereby Pauli blocking the decay and effectively rendering DM
stable on cosmological timescales. We carefully study the stability of the
scalar DM condensate by taking account of decay, scattering, and parametric
resonance effects, which turn out to be important during different eras of the
DM density evolution. In particular, we highlight the importance of DM
scattering, including the self-annihilation of condensate particles and their
scattering with thermalized scalars, which proceeds via the same coupling that
mediates the DM decay. Such scattering processes can lead to the fast
evaporation of the condensate and thus significantly constrain this scenario
at large couplings. Nevertheless, we demonstrate that over a broad range of DM
masses and at moderate values of the coupling constant, a viable scenario of
DM stabilization due to Pauli blocking can be realized.

For relatively sizable Yukawa couplings, our scenario leads to a novel late
time cosmology. The decay and scattering of DM can lead to the formation of a
quasi-thermal component of fermions and scalars that behaves as dark radiation
with a temperature far above the mass of the DM. The resulting dark radiation
can have a substantial abundance and, furthermore, is self-interacting since
it attains kinetic equilibrium via the decay and inverse decay of the
relativistic dark sector plasma. The deviation of the effective neutrino
number at recombination can be as large as $\mathcal{O}(0.1)$ with only small
corrections to the DM abundance relative to $\Lambda$CDM. Very interestingly,
such a scenario may have important implications for the more than $5\sigma$
tension between local and global measurements of the Hubble constant [1, 2, 3,
4, 5, 6, 7]. (See also Refs. [8, 9], and reviews [10, 11, 12, 13] for the
Hubble tension.) Previous studies have shown that the Hubble tensions may be
alleviated due to the presence of self-interacting dark radiation [14, 15, 12,
16]. Furthermore, even if the dark sector only couples to the SM via gravity,
the resulting DM halo density profiles will feature a cutoff since the
formation of dense regions of DM will be obstructed by fast DM evaporation.
Thus our scenario may prefer cored DM density profiles with a cutoff depending
on the masses and coupling strength. Precisely simulating structure formation
and measuring the small-scale structure, e.g., via halo-mass function
measurements in Vera Rubin Observatory [17] or 21 cm line observation [18],
can thus provide probes of this scenario.

An obvious candidate for the fermions in our scenario are the Standard Model
(SM) neutrinos, which leads to several additional consequences in comparison
to the case of an exotic BSM fermion. In this case, the composition of the
cosmic neutrino background (C$\nu$B) may be altered if the coupling of the DM
to neutrinos is not diagonal in the neutrino mass eigenbasis. This is because
the heavier neutrinos may decay into DM and a lighter neutrino. The produced
DM, which is energetic and thus not Pauli blocked, can further decay into
lighter neutrinos. The cascade decays result in a C$\nu$B that is composed
entirely of the lightest neutrinos, some of which are boosted. Interestingly,
in the case of normal ordering our scenario predicts up to a factor 1-10
enhancement of the boosted neutrino interaction rate to an electron by the
neutrino decays with an energy of $\mathcal{O}(0.01){\rm\,eV}$. In addition,
either in the normal or inverted ordering, the DM may further produce the
lightest neutrinos via the evaporation or decay while still maintaining
stability of the DM condensate. For certain regions of parameter space,
although the typical energy of each neutrino is smaller than
$\mathcal{O}(0.01){\rm\,eV}$, this contribution can enhance the C$\nu$B
interaction rate by many orders of magnitude. Thus, the imprint of DM
interactions on the C$\nu$B may potentially be probed in future experiments
such as PTOLEMY [19, 20, 21].

Before moving to the main part of the paper, let us mention a few relevant
studies. Ref. [22] investigated the possibility of scalar DM decaying to
neutrinos and the impact of Pauli blocking. However, this work did not
consider the effect of scattering, which turns out to be important even for
modest Yukawa couplings. The cosmological production of fermions due to scalar
field dynamics, either from broad parametric resonance [23, 24] or
perturbative decays [25, 26, 27], has been widely investigated in the context
of inflationary reheating or preheating. While distinct from our scenario,
light dark sectors in the form of a nearly degenerate Fermi sea has been
explored in Refs. [28, 26, 29, 30].

This paper is organized as follows. In the next section, we discuss our
theoretical setup and in the perturbative regime. We formulate a system of
Boltzmann equations, including quantum statistical effects, which allows us to
account for the Pauli blocking effect and investigate the stability of the DM
condensate. (Some technical details are presented in Appendices A and B). In
Sec. 3, we investigate the minimal scenario in which DM is stabilized by
exotic BSM fermions. We study the cases of DM annihilation with and without a
chirality suppression, the production of dark radiation and the DM density
evolution, the potential implications for the Hubble tensions, the cutoff of
the present local DM density, the broad parametric resonance regime, and
conclude with a summary of the cosmology of this scenario. In Sec. 4 we
explore the scenario in which the C$\nu$B stabilizes DM and the resulting
implications for the evolution and detection of the the C$\nu$B. (A gauge-
invariant UV model is discussed in Appendix C) In Sec. 5 we discuss some
generalizations to other models and signatures. The final section is devoted
to the conclusions and outlook.

##  2 Setup and Boltzmann equations

###  2.1 Setup

To clarify our idea, we introduce a minimal model containing a real scalar
field $\phi$ with mass $m_{\phi}$, which will be identified as the DM, and a
Majorana fermion $\psi$ with mass $m_{\psi}$. Throughout this work we will
consider $m_{\psi}<m_{\phi}/2$, and for the moment we neglect the effect of
the fermion mass.111We note that fermionic bound states [31, 22, 32, 33] do
not form in our scenario since the range of the Yukawa force is much smaller
than the fermion de Broglie wavelength.  Here and in Sec. 3 we assume that
$\psi$ is a BSM fermion with negligible couplings to SM particles, while in
Sec. 4 we will examine the case in which $\psi$ are the SM neutrinos. The
scalar-fermion interaction is represented by the Lagrangian

${\cal L}\supset-\frac{y}{2}\phi\bar{\psi}\psi\,,$ (1)

where $y$ is a real Yukawa coupling.222 We have also considered the case of a
pseudoscalar coupling, see Appendix A. The dynamics do not change much, but
the annihilation of $\phi$ condensates proceeds in the $D$-wave, and the dark
radiation production is similar to the chirality-suppressed case examined in
Sec. 3.1 even if $m_{\psi}=0$.  As a consequence of this interaction, $\phi$
decays to a $\psi$ pair in the vacuum.

Suppose that at the redshift $z=z_{\rm ini}$, there is a population of non-
relativistic $\phi$ particles with number density $n_{\phi}\neq 0$ and
negligible population of $\psi$ particles, $n_{\psi}=0$. Then, $\phi$ decays
via the process $\phi\to\psi\psi$, with decay rate

$\Gamma_{\phi\to\psi\psi}\simeq\frac{y^{2}m_{\phi}}{16\pi}.$ (2)

For simplicity we neglect the mass of $\psi$ here and in the following (see
Appendix A for the formula without this approximation). If the decay is fast
enough, the states of $\psi$ in the spherical shell of momentum around
$E_{\psi}\simeq m_{\phi}/2$ are soon occupied. Strictly speaking, there are
also effects due to scattering processes, e.g., $\phi\phi\to\psi\psi$, which
will turn out to be important and will be the main focus of this paper. The
width of the shell, $\Delta E$, depends on the velocity dispersion of $\phi$,
for which we may safely assume $\Delta E\ll m_{\phi}/2$ if $\phi$ is cold
enough. For concreteness, we assume that $\phi$ is produced through the
misalignment mechanism and that $\Delta E$ is negligible [34, 35,
36].333Strictly speaking, the momentum distribution has a width of $Qm_{\phi}$
in the misalignment case [23, 24], where $Q$ is defined in Eq.(9).  Our setup
will work with other DM production mechanisms as well.

Then the backreaction, i.e., Pauli blocking and inverse decay, becomes
important. By accounting for the fermion phase space distribution,
$f_{\psi}[E_{\psi}]$, the effective decay rate of the DM is corrected as (see
also Eq. (37) below)

$\Gamma^{\rm
eff}_{\phi\to\psi\psi}\simeq\frac{y^{2}m_{\phi}}{16\pi}(1-2f_{\psi}[m_{\phi}/2]).$
(3)

Consequently $f_{\psi}$ cannot increase beyond $f_{\psi}=1/2$ around which the
$\phi$ decay is forbidden, i.e., the decay is Pauli blocked [22].

Since the Universe expands, the energy of the Fermi shell decreases as
$E_{\psi}\propto a^{-1}.$ As the initial shell shrinks, the $\phi$ decay
continually replenishes the shell of $E_{\psi}\simeq m_{\phi}/2$. As a
consequence, a distribution of $f_{\psi}$ at $z<z_{\rm ini}$ is formed as

$f_{\psi}=\frac{1}{2}\Theta\left[\frac{m_{\phi}}{2}-p\right]\times\Theta\left[p-\frac{m_{\phi}}{2}\frac{(1+z)}{(1+z_{\rm
ini})}\right]\to\frac{1}{2}\Theta\left[\frac{m_{\phi}}{2}-p\right],$ (4)

where $\Theta$ is a step function. For a sufficiently long period, $z\ll
z_{\rm ini}$, the $\psi$ form a Fermi sea, characterized by a Fermi sphere of
energy $E_{F}\simeq m_{\phi}/2$ as shown in the r.h.s. of Eq. (4).

Even though the DM is continuously decaying to fill the Fermi sea, the energy
transferred from the DM condensate to the Fermi sea is negligible if444 The
motion and gravitational potential, e.g. during structure formation, may
slighly alter the momenta of $\psi$. As long as the effects do not change the
$\psi$ phase distribution by a substantial amount, our conclusions do not
change as long as Eq.(5) is satisfied.

$\rho_{\phi}\gg\frac{g_{\psi}}{2}\times\frac{E_{F}^{4}}{8\pi^{2}}\simeq
g_{\psi}\frac{m_{\phi}^{4}}{256\pi^{2}},$ (5)

where the r.h.s is the energy density of the Fermi sphere and $g_{\psi}=2$
denotes the spin degree of freedom of $\psi.$ For $\phi$ to account for the
cosmologically stable DM, this relation should be at least satisfied in the
present Universe. By taking $\rho_{\phi}=\rho_{\rm DM}\simeq(2\,{\rm
meV})^{4},$ which is the DM energy density today [37], we obtain the following
constraint on the DM mass in our scenario:

$m_{\phi}\ll 0.01{\rm\,eV}.$ (6)

Therefore, what we will mainly discuss is a sub-eV DM candidate. In the
following, we will see that additional reactions beyond decays can lead to the
formation of a quasi-thermal component of scalars and fermions with
temperature much higher than $m_{\phi}$, which, however, may still be
subdominant to the DM.

So far we have discussed the stability in the leading order approximation of
$y$. As we have mentioned several times, a higher-order effect, i.e.,
$2$-to-$2$ scatterings with rates scaling as $\mathcal{O}(y^{4})$, is quite
important for discussing the stability of the DM condensate. If either the
decay or scattering (including DM annihilation) were absent and the condition
Eq.(5) is parametrically satisfied, there would be both comoving number and
energy conservations555In the ${1\leftrightarrow 2}$ process,
$n_{\phi}+n_{\psi}/2$ is conserved, and in the $2\leftrightarrow 2$ process
$n_{\phi}+n_{\psi}$ is conserved with $n_{\psi}$ being the fermion+antifermion
number.  which stabilize the DM condensate.

The annihilation/scattering rate of a plasma $\phi$ or $\psi$ particle of
energy $E\gg m_{\phi}$ with a condensate $\phi$, neglecting the Pauli blocking
effect, is given by

$v\sigma_{0}n_{\phi}\approx\frac{y^{4}}{\pi m_{\phi}E}n_{\phi}.$ (7)

where $\sigma_{0}$ is scattering cross section. For now, let us consider the
regime

$\Gamma_{\phi\to\psi\psi}\frac{m_{\phi}}{E}\gg v\sigma_{0}n_{\phi}\gg H,$ (8)

where the l.h.s is the (inverse) decay rate of the relativistic $(\psi)\phi$
of energy $E$. The first inequality leads to the condition

$Q\equiv\frac{y^{2}\rho_{\phi}}{m_{\phi}^{4}}\ll 1.$ (9)

This corresponds to the narrow resonance regime in the context of parametric
resonance [23, 24]. We will come back to the broad resonance regime, which can
be important at early times, in Sec.3.4. In the narrow resonance regime, one
can study the system, including the effects of both decays/inverse decays and
scatterings, by employing Boltzmann equations.

Before moving on, we remark that light scalar field will generically receive
radiative corrections to its potential which are sensitive to UV physics.
These corrections could be small if the dark scalar and fermion have only weak
couplings to UV degrees of freedom. Alternatively, one can consider $\phi$ to
be a pseudo-Goldstone boson with interaction
$\partial_{\mu}\phi\overline{\psi}\gamma^{\mu}\gamma^{5}\psi$ instead of
Eq.(1), in which case the corrections to the potential are protected by a
shift symmetry. Additionally, the free energy density of fermions sources a
contribution to the effective potential of the scalar field. While this
contribution is always subdominant in comparison to the scalar mass term for
the minimal model considered here and in Sec. 4, it can be relevant at early
times when the fermion is the SM neutrino due to the presence of the C$\nu$B,
see Sec. 5 for further discussion.

###  2.2 Boltzmann Equations

We now investigate the impact of scattering between the dark sector particles,
which can cause the DM condensate to evaporate and lead to the formation of a
thermalized dark sector component that serves as dark radiation. To
investigate the stability of the DM, estimate its lifetime, and discern the
properties of the dark radiation component, we will use the formalism of
Boltzmann equations. As will be argued in Sec. 3, in large regions of
parameter space it is expected that the system reaches a quasi-equilibrium
state, which motivates us to employ a specific ansatz for the form of the
phase space distribution functions that separates the condensate and the
thermalized components. For the latter, we take for the dark radiation $\psi$
and $\phi$ particles the standard equilibrium Fermi-Dirac and Bose-Einstein
distributions, which should provide a qualitatively correct description as the
system thermalizes.

We wish to follow the evolution of the phase space distributions of the scalar
and fermion fields, $f_{\phi}[p_{\phi},t]$ and $f_{\psi}[p_{\psi},t]$. These
are governed by the Boltzmann equations,

$\displaystyle\frac{\partial f_{\phi}[p_{\phi},t]}{\partial
t}-p_{\phi}H\frac{\partial f_{\phi}[p_{\phi},t]}{\partial p_{\phi}}$
$\displaystyle=C^{\phi}[p_{\phi},t],$ (10) $\displaystyle\frac{\partial
f_{\psi}[p_{\psi},t]}{\partial t}-p_{\psi}H\frac{\partial
f_{\psi}[p_{\psi},t]}{\partial p_{\psi}}$
$\displaystyle=C^{\psi}[p_{\psi},t],$ (11)

where $C^{\phi}$ and $C^{\psi}$ are the collision terms to be described in
detail below. Following the motivation outlined above, we make the following
ansatz for the phase space distribution functions:

$\displaystyle f_{\phi}[p_{\phi},t]$
$\displaystyle=\hat{f}_{\phi}[p_{\phi},t]+f_{\phi}^{\rm th}[p_{\phi},t],$ (12)
$\displaystyle f_{\psi}[p_{\psi},t]$ $\displaystyle=f_{\psi}^{\rm
th}[p_{\psi},t],$ (13)

where

$\displaystyle\hat{f}_{\phi}[p_{\phi},t]$
$\displaystyle=(2\pi)^{3}\hat{n}_{\phi}[t]\delta^{3}(\vec{p}_{\phi})=(2\pi)^{3}\hat{n}_{\phi}[t]\left[\frac{\delta(p_{\phi})}{4\pi
p_{\phi}^{2}}\right],$ (14) $\displaystyle f_{\phi}^{\rm th}[p_{\phi},t]$
$\displaystyle=\left\\{\exp\left[\frac{E_{\phi}-\mu_{\phi}[t]}{T[t]}\right]-1\right\\}^{-1},$
(15) $\displaystyle f_{\psi}^{\rm th}[p_{\psi},t]$
$\displaystyle=\left\\{\exp\left[\frac{E_{\psi}-\mu_{\psi}[t]}{T[t]}\right]+1\right\\}^{-1},$
(16)

with $E_{i}=\sqrt{p_{i}^{2}+m_{i}^{2}}$ for $i=\phi,\psi$. In particular, we
have separated the scalar phase-space distribution into condensate (denoted by
the hat) and thermal components, and we have assumed the thermal components
are well-described by the standard equilibrium forms with temperature $T[t]$
and chemical potential $\mu_{i}[t]$. We note that $\mu_{i}[t]$ is not assumed
to take any specific value associated with chemical/kinetic equilibrium.
Rather, we will follow its dynamical evolution from generic initial
conditions. As we will show below $\mu_{\phi}[t]=m_{\phi}$ is the condition
for the system to reach kinetic equilibrium.

It will be convenient to follow the integrated quantities such as the number
density, energy density, and pressure, defined as

$\displaystyle n_{i}[t]$
$\displaystyle=g_{i}\\!\int\\!\frac{d^{3}p_{i}}{(2\pi)^{3}}\,f_{i}[p_{i},t],$
(17) $\displaystyle\rho_{i}[t]$
$\displaystyle=g_{i}\\!\int\\!\frac{d^{3}p_{i}}{(2\pi)^{3}}\,E_{i}\,f_{i}[p_{i},t],$
(18) $\displaystyle P_{i}[t]$
$\displaystyle=g_{i}\\!\int\\!\frac{d^{3}p_{i}}{(2\pi)^{3}}\,\frac{p_{i}}{3E_{i}}\,f_{i}[p_{i},t],$
(19)

for $i=\phi,\psi$. Here $g_{i}$ counts the number of internal degrees of
freedom for species $i$ ($g_{\phi}=1$, $g_{\psi}=2$). In fact, for the scalar
$\phi$ we can isolate the condensate and thermal components by defining a type
of projection operator $\hat{P}$ that excludes a small region around
$p_{\phi}=0$ in the $p_{\phi}$ integral. In particular, we define

$\displaystyle\hat{P}\\!\int\\!\frac{d^{3}p_{\phi}}{(2\pi)^{3}}$
$\displaystyle\equiv\frac{1}{(2\pi)^{3}}\int
d\Omega_{\phi}\int_{\epsilon}^{\infty}p_{\phi}^{2}\,dp_{\phi},$ (20)
$\displaystyle[1-\hat{P}]\\!\int\\!\frac{d^{3}p_{\phi}}{(2\pi)^{3}}$
$\displaystyle\equiv\frac{1}{(2\pi)^{3}}\int
d\Omega_{\phi}\int_{0}^{\epsilon}p_{\phi}^{2}\,dp_{\phi},$ (21)

where $\epsilon$ is a small parameter. We then have the relations

$\displaystyle[1-\hat{P}]\\!\int\\!\frac{d^{3}p_{\phi}}{(2\pi)^{3}}\,f_{\phi}[p_{\phi},t]$
$\displaystyle=[1-\hat{P}]\\!\int\\!\frac{d^{3}p_{\phi}}{(2\pi)^{3}}\,\hat{f}_{\phi}[p_{\phi},t]=\hat{n}_{\phi}[t],$
(22)
$\displaystyle\hat{P}\\!\int\\!\frac{d^{3}p_{\phi}}{(2\pi)^{3}}\,f_{\phi}[p_{\phi},t]$
$\displaystyle=\hat{P}\\!\int\\!\frac{d^{3}p_{\phi}}{(2\pi)^{3}}\,f_{\phi}^{\rm
th}[p_{\phi},t]=n_{\phi}^{\rm th}[t].$ (23)

Using the projectors, the Boltzmann equations for the number densities can
then be written as

$\displaystyle\dot{\hat{n}}_{\phi}+3H\hat{n}_{\phi}$ $\displaystyle={\cal
C}^{\phi_{c}},$ (24) $\displaystyle\dot{n}_{\phi}^{\rm th}+3Hn_{\phi}^{\rm
th}$ $\displaystyle={\cal C}^{\phi},$ (25)
$\displaystyle\dot{n}_{\psi}+3Hn_{\psi}$ $\displaystyle={\cal C}^{\psi},$ (26)

where the ${\cal C}$s on the r.h.s represent the corresponding integrated
collision terms, which will be discussed below. Similarly, we can write the
Boltzmann equations for the energy densities as

$\displaystyle\dot{\hat{\rho}}_{\phi}+3H(\hat{\rho}_{\phi}+\hat{P}_{\phi})$
$\displaystyle={\cal E}^{\phi_{c}},$ (27) $\displaystyle\dot{\rho}_{\phi}^{\rm
th}+3H(\rho_{\phi}^{\rm th}+P_{\phi}^{\rm th})$ $\displaystyle={\cal
E}^{\phi},$ (28) $\displaystyle\dot{\rho}_{\psi}+3H(\rho_{\psi}+P_{\psi})$
$\displaystyle={\cal E}^{\psi},$ (29)

where the ${\cal E}$s on the r.h.s represent the corresponding integrated
collision terms. Summing the three equations above, we can write a Boltzmann
equation for the total energy density:

$\displaystyle\dot{\rho}_{\rm tot}+3H(\rho_{\rm tot}+P_{\rm tot})$
$\displaystyle={\cal E}^{\phi_{c}}+{\cal E}^{\phi}+{\cal E}^{\psi}=0,$ (30)

where $\rho_{\rm tot}\equiv\hat{\rho}_{\phi}+\rho_{\phi}^{\rm
th}+\rho_{\psi}$, $P_{\rm tot}\equiv\hat{P}_{\phi}+P_{\phi}^{\rm
th}+P_{\psi}$, and the last equality follows from energy conservation. In
particular, for the condensate we have

$\hat{\rho}_{\phi}=m_{\phi}\hat{n}_{\phi},~{}~{}~{}~{}~{}~{}~{}\hat{P}_{\phi}=0.$
(31)

On the other hand, the thermal $\phi$ and $\psi$ components are typically
relativistic, leading to

$P_{\phi}^{\rm th}\simeq\rho_{\phi}^{\rm
th}/3,~{}~{}~{}~{}~{}~{}~{}P_{\psi}\simeq\rho_{\psi}/3.$ (32)

We now turn to the collision terms for the various reactions.

####  2.2.1 Decay/inverse decay $\phi\leftrightarrow\psi\psi$

We first consider the decay/inverse decay process

$\phi(p_{\phi})\leftrightarrow\psi(p_{1})+\psi(p_{2}).$ (33)

The unintegrated $\phi$ collision term
$C^{\phi}_{\phi\leftrightarrow\psi\psi}$ is

$\displaystyle C^{\phi}_{\phi\leftrightarrow\psi\psi}$
$\displaystyle=-\frac{1}{S_{\psi}}\frac{1}{g_{\phi}}\frac{1}{2E_{\phi}}\sum_{\rm
spins}\int\\!d\Pi_{1}\,d\Pi_{2}\,(2\pi)^{4}\delta^{4}(p_{\phi}-p_{1}-p_{2})|{\cal
M}_{\phi\rightarrow\psi\psi}|^{2}$ (34)
$\displaystyle\times\left\\{f_{\phi}[p_{\phi}](1-f_{\psi}[p_{1}])(1-f_{\psi}[p_{2}])-(1+f_{\phi}[p_{\phi}])f_{\psi}[p_{1}]f_{\psi}[p_{2}]\right\\}.$

Here $1/S_{\psi}=1/2$ is a symmetry factor accounting for identical $\psi$
particles, $g_{\phi}=1$, and

$\displaystyle d\Pi_{i}\equiv\frac{d^{3}p_{i}}{(2\pi)^{3}2E_{i}}.$ (35)

We now integrate over the $\phi$ phase space to derive the integrated
collision terms, using the projection operators $[1-\hat{P}]$ and $\hat{P}$ to
isolate the condensate and thermal components, respectively. For the
condensate term, we obtain

$\displaystyle{\cal C}^{\phi_{c}}_{\phi_{c}\leftrightarrow\psi\psi}$
$\displaystyle=g_{\phi}[1-\hat{P}]\\!\int\\!\frac{d^{3}p_{\phi}}{(2\pi)^{3}}C^{\phi}_{\phi\leftrightarrow\psi\psi}$
(36) $\displaystyle=-\frac{1}{S_{\psi}}\sum_{\rm
spins}[1-\hat{P}]\\!\int\\!d\Pi_{\phi}\,d\Pi_{1}\,d\Pi_{2}\,(2\pi)^{4}\delta^{4}(p_{\phi}-p_{1}-p_{2})|{\cal
M}_{\phi\rightarrow\psi\psi}|^{2}$
$\displaystyle\times\hat{f}_{\phi}[p_{\phi}](1-f_{\psi}[p_{1}]-f_{\psi}[p_{2}]).$

All integrals can be carried out analytically, leading to the simple
expression

$\displaystyle{\cal
C}^{\phi_{c}}_{\phi_{c}\leftrightarrow\psi\psi}=-\Gamma_{\phi\to\psi\psi}\left(1-2f_{\psi}[E_{\psi}=m_{\phi}/2]\right)\hat{n}_{\phi},$
(37)

where $\Gamma_{\phi\to\psi\psi}$ is given in (2) (See Appendix A or Eq. (105)
for the mixed real and pseudo-scalar coupling case and full mass dependence).
We note that the collision term (37) reproduces the effective $\phi$
condensate decay rate advertised above in Eq.(3).

Next, for the thermal component we obtain

$\displaystyle{\cal C}^{\phi}_{\phi\leftrightarrow\psi\psi}$
$\displaystyle=g_{\phi}\hat{P}\\!\int\\!\frac{d^{3}p_{\phi}}{(2\pi)^{3}}C^{\phi}_{\phi\leftrightarrow\psi\psi}$
(38) $\displaystyle=-\frac{1}{S_{\psi}}\sum_{\rm
spins}\hat{P}\\!\int\\!d\Pi_{\phi}\,d\Pi_{1}\,d\Pi_{2}\,(2\pi)^{4}\delta^{4}(p_{\phi}-p_{1}-p_{2})|{\cal
M}_{\phi\rightarrow\psi\psi}|^{2}$ $\displaystyle\times\left\\{f^{\rm
th}_{\phi}[p_{\phi}](1-f_{\psi}[p_{1}])(1-f_{\psi}[p_{2}])-(1+f^{\rm
th}_{\phi}[p_{\phi}])f_{\psi}[p_{1}]f_{\psi}[p_{2}]\right\\}.$

We can reduce the collision term to the following double integral expression:

$\displaystyle{\cal C}^{\phi}_{\phi\leftrightarrow\psi\psi}$
$\displaystyle=-\frac{m_{\phi}\Gamma_{\phi\to\psi\psi}}{2\pi^{2}\beta_{\psi}}\int_{m_{\phi}}^{\infty}\\!dE_{\phi}\int_{E_{1}^{-}}^{E_{1}^{+}}\\!dE_{1}\times\left\\{f^{\rm
th}_{\phi}[p_{\phi}](1-f_{\psi}[p_{1}]-f_{\psi}[p_{2}])-f_{\psi}[p_{1}]f_{\psi}[p_{2}]\right\\},$
(39)

where the integration limits $E_{1}^{\pm}$ are (see Appendix B)

$\displaystyle E_{1}^{\pm}$
$\displaystyle=\frac{1}{2}(E_{\phi}\pm\beta_{\psi}p_{\phi}).$ (40)

The phase space distributions here are regarded as functions of the energies
$E_{i}=\sqrt{p_{i}^{2}+m_{i}^{2}}$, and the energy conservation condition
$E_{\phi}=E_{1}+E_{2}$ implies
$f_{\psi}[p_{2}]=f_{\psi}[E_{2}=E_{\phi}-E_{1}]$.

Following similar steps, the integrated collision term ${\cal
C}^{\psi}_{\phi\leftrightarrow\psi\psi}$ can be related to the $\phi$
collision terms above as

$\displaystyle{\cal
C}^{\psi}_{\phi\leftrightarrow\psi\psi}=-N_{\psi}\,[\,{\cal
C}^{\phi_{c}}_{\phi_{c}\leftrightarrow\psi\psi}+{\cal
C}^{\phi}_{\phi\leftrightarrow\psi\psi}\,],$

where the factor $N_{\psi}=2$ is present because two $\psi$ particles are lost
in the reaction.

We mention that ${\cal C}^{\phi_{c}}_{\phi_{c}\leftrightarrow\psi\psi},-{\cal
C}^{\psi}_{\phi_{c}\leftrightarrow\psi\psi}\propto
e^{m_{\phi}/2T}-e^{\mu_{\psi}/T}$ biases $\mu_{\psi}\to m_{\phi}/2$, at which
point the DM condensate decay is Pauli blocked. Furthermore, ${\cal
C}^{\phi}_{\phi\leftrightarrow\psi\psi},-{\cal
C}^{\psi}_{\phi\leftrightarrow\psi\psi}\propto
e^{\mu_{\phi}/T}-e^{2\mu_{\psi}/T}$ biases $\mu_{\phi}\to 2\mu_{\psi}$.
Therefore the chemical/kinetic equilibrium for decays/inverse decays is
achieved for $\mu_{\psi}=\mu_{\phi}/2=m_{\phi}/2.$

####  2.2.2 Annihilation/inverse annihilation
$\phi\phi\leftrightarrow\psi\psi$

Next, we consider the process

$\phi(p_{\phi})+\phi(p_{2})\leftrightarrow\psi(p_{1})+\psi(p_{2}).$ (41)

The unintegrated $\phi$ collision term
$C^{\phi}_{\phi\phi\leftrightarrow\psi\psi}$ is

$\displaystyle C^{\phi}_{\phi\phi\leftrightarrow\psi\psi}$
$\displaystyle=-\frac{1}{S_{\psi}}\frac{1}{g_{\phi}}\frac{1}{2E_{\phi}}\sum_{\rm
spins}\int\\!d\Pi_{2}\,d\Pi_{3}\,d\Pi_{4}\,(2\pi)^{4}\delta^{4}(p_{\phi}+p_{2}-p_{3}-p_{4})|{\cal
M}_{\phi\phi\rightarrow\psi\psi}|^{2}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle\times\left\\{f_{\phi}[p_{\phi}]f_{\phi}[p_{2}](1-f_{\psi}[p_{3}])(1-f_{\psi}[p_{4}])-(1+f_{\phi}[p_{\phi}])(1+f_{\phi}[p_{2}])f_{\psi}[p_{3}]f_{\psi}[p_{4}]\right\\}.$
(42)

Next, we integrate over the $\phi$ phase space to obtain the integrated
collision terms, again employing the projection operators $[1-\hat{P}]$ and
$\hat{P}$ to isolate the condensate and thermal components, respectively. We
also include multiplicative factors $N_{\phi}=2$ (two $\phi$ particles are
lost) and a factor $1/S_{\phi}=1/2$ (identical $\phi$ particles). There are
two contributions to the condensate collision term. First, we have a term
accounting for the annihilation of two condensate particles:

$\displaystyle{\cal C}^{\phi_{c}}_{\phi_{c}\phi_{c}\leftrightarrow\psi\psi}$
$\displaystyle=\frac{N_{\phi}}{S_{\phi}}g_{\phi}[1-\hat{P}]_{\phi}\\!\int\\!\frac{d^{3}p_{\phi}}{(2\pi)^{3}}[1-\hat{P}]_{2}C^{\phi}_{\phi\phi\leftrightarrow\psi\psi}$
(43) $\displaystyle=-\frac{N_{\phi}}{S_{\phi}S_{\psi}}\sum_{\rm
spins}[1-\hat{P}]_{\phi}[1-\hat{P}]_{2}\int\\!d\Pi_{\phi}\,d\Pi_{2}\,d\Pi_{3}\,d\Pi_{4}$
$\displaystyle\times(2\pi)^{4}\delta^{4}(p_{\phi}+p_{2}-p_{3}-p_{4})|{\cal
M}_{\phi\phi\rightarrow\psi\psi}|^{2}\,\hat{f}_{\phi}[p_{\phi}]\hat{f}_{\phi}[p_{2}](1-f_{\psi}[p_{3}]-f_{\psi}[p_{4}]),$

where the subscript on the projector refers to the $\phi$ momentum. All
integrals can be carried out analytically, leading to the simple expression

$\displaystyle{\cal
C}^{\phi_{c}}_{\phi_{c}\phi_{c}\leftrightarrow\psi\psi}=-\frac{N_{\phi}}{S_{\phi}S_{\psi}}\frac{1}{32\pi
m_{\phi}^{2}}\sum_{\rm spins}|{\cal
M}_{\phi\phi\rightarrow\psi\psi}|^{2}\bigg{|}_{\vec{p}_{\phi}=\vec{p}_{2}=0}\left[1-\frac{m_{\psi}^{2}}{m_{\phi}^{2}}\right]^{1/2}\\!\\!\left(1-2f_{\psi}[E_{\psi}=m_{\phi}]\right)\hat{n}^{2}_{\phi}.$
(44)

There is also a collision term accounting for the annihilation of a condensate
scalar with a thermal one:

$\displaystyle{\cal C}^{\phi_{c}}_{\phi_{c}\phi\leftrightarrow\psi\psi}$
$\displaystyle=\frac{N_{\phi}}{S_{\phi}}g_{\phi}[1-\hat{P}]_{\phi}\\!\int\\!\frac{d^{3}p_{\phi}}{(2\pi)^{3}}\hat{P}_{2}\,C^{\phi}_{\phi\phi\leftrightarrow\psi\psi}$
(45) $\displaystyle=-\frac{N_{\phi}}{S_{\phi}S_{\psi}}\sum_{\rm
spins}[1-\hat{P}]_{\phi}\,\hat{P}_{2}\,\int\\!d\Pi_{\phi}\,d\Pi_{2}\,d\Pi_{3}\,d\Pi_{4}\,(2\pi)^{4}\delta^{4}(p_{\phi}+p_{2}-p_{3}-p_{4})\,$
$\displaystyle\times|{\cal
M}_{\phi\phi\rightarrow\psi\psi}|^{2}\,\hat{f}_{\phi}[p_{\phi}]\left\\{f_{\phi}^{\rm
th}[p_{2}](1-f_{\psi}[p_{3}]-f_{\psi}[p_{4}])-f_{\psi}[p_{3}]f_{\psi}[p_{4}]\right\\},~{}~{}~{}~{}~{}$

We can reduce the collision term to the following double integral expression:

$\displaystyle{\cal C}^{\phi_{c}}_{\phi_{c}\phi\leftrightarrow\psi\psi}$
$\displaystyle=-\frac{N_{\phi}}{S_{\phi}S_{\psi}}\frac{\hat{n}_{\phi}}{64\pi^{3}m_{\phi}}\sum_{\rm
spins}\int_{m_{\phi}}^{\infty}\\!dE_{2}\int_{E_{3}^{-}}^{E_{3}^{+}}\\!dE_{3}\,|{\cal
M}_{\phi\phi\rightarrow\psi\psi}|^{2}$
$\displaystyle\times\left\\{f_{\phi}^{\rm
th}[p_{2}](1-f_{\psi}[p_{3}]-f_{\psi}[p_{4}])-f_{\psi}[p_{3}]f_{\psi}[p_{4}]\right\\},$
(46)

where the integration limits $E_{3}^{\pm}$ are given (see Appendix B):

$\displaystyle E_{3}^{\pm}$
$\displaystyle=\frac{1}{2}(E_{2}+m_{\phi}\pm\kappa_{\psi}\,p_{2}),~{}~{}~{}~{}\kappa_{\psi}=\left(1-\frac{2m_{\psi}^{2}}{m_{\phi}(E_{2}+m_{\phi})}\right)^{1/2}.$
(47)

Again, the phase space distributions here are regarded as functions of the
energies $E_{i}=\sqrt{p_{i}^{2}+m_{i}^{2}}$, and the energy conservation
condition $m_{\phi}+E_{2}=E_{3}+E_{4}$ implies
$f_{\psi}[p_{4}]=f_{\psi}[E_{4}=m_{\phi}+E_{2}-E_{3}]$.

Next, it is straightforward to see that the integrated collision term of the
thermal $\phi$ component, which arises from the annihilation with a condensate
scalar, is identical to ${\cal
C}^{\phi_{c}}_{\phi_{c}\phi\leftrightarrow\psi\psi}$ considered above:

$\displaystyle{\cal C}^{\phi}_{\phi\phi_{c}\leftrightarrow\psi\psi}={\cal
C}^{\phi_{c}}_{\phi_{c}\phi\leftrightarrow\psi\psi}.$ (48)

We also note that there is a collision term involving the annihilation of two
thermal $\phi$ particles, ${\cal C}^{\phi}_{\phi\phi\leftrightarrow\psi\psi}$.
We neglect the effect of this term since it is subdominant to the other terms.

Finally, we consider the integrated $\psi$ collision term ${\cal
C}^{\psi}_{\phi\phi\leftrightarrow\psi\psi}$, which can be related to the
$\phi$ collision terms as follows:

$\displaystyle{\cal C}^{\psi}_{\phi\phi\leftrightarrow\psi\psi}$
$\displaystyle=-\left[{\cal
C}^{\phi_{c}}_{\phi_{c}\phi_{c}\leftrightarrow\psi\psi}+{\cal
C}^{\phi_{c}}_{\phi_{c}\phi\leftrightarrow\psi\psi}+{\cal
C}^{\phi}_{\phi\phi_{c}\leftrightarrow\psi\psi}+{\cal
C}^{\phi}_{\phi\phi\leftrightarrow\psi\psi}\right]$ (49)

We mention that ${\cal
C}^{\phi_{c}}_{\phi_{c}\phi\leftrightarrow\psi\psi},-{\cal
C}^{\psi}_{\phi_{c}\phi\leftrightarrow\psi\psi}\propto
e^{(m_{\phi}+\mu_{\phi})/T}-e^{2\mu_{\psi}/T}$ biases
$\mu_{\psi}\to(\mu_{\phi}+m_{\phi})/2$, and ${\cal
C}^{\phi_{c}}_{\phi_{c}\phi_{c}\leftrightarrow\psi\psi},-{\cal
C}^{\psi}_{\phi_{c}\phi_{c}\leftrightarrow\psi\psi}\propto
e^{m_{\phi}/T}-e^{\mu_{\psi}/T}$ biases $\mu_{\psi}\to\mu_{\phi}$. Therefore
the chemical/kinetic equilibrium for the annihilation processes is realized
when $\mu_{\phi}=\mu_{\psi}=m_{\phi}.$

####  2.2.3 Scattering $\phi\psi\leftrightarrow\phi\psi$

Next, we consider the scattering process

$\phi(p_{\phi})+\psi(p_{2})\leftrightarrow\phi(p_{3})+\psi(p_{4}).$ (50)

The unintegrated $\phi$ collision term
$C^{\phi}_{\phi\psi\leftrightarrow\phi\psi}$ is

$\displaystyle C^{\phi}_{\phi\psi\leftrightarrow\phi\psi}$
$\displaystyle=-\frac{1}{g_{\phi}}\frac{1}{2E_{\phi}}\sum_{\rm
spins}\int\\!d\Pi_{2}\,d\Pi_{3}\,d\Pi_{4}\,(2\pi)^{4}\delta^{4}(p_{\phi}+p_{2}-p_{3}-p_{4})|{\cal
M}_{\phi\psi\rightarrow\phi\psi}|^{2}~{}~{}~{}~{}~{}~{}~{}~{}$
$\displaystyle\times\left\\{f_{\phi}[p_{\phi}]f_{\psi}[p_{2}](1+f_{\phi}[p_{3}])(1-f_{\psi}[p_{4}])-(1+f_{\phi}[p_{\phi}])(1-f_{\psi}[p_{2}])f_{\phi}[p_{3}]f_{\psi}[p_{4}]\right\\}.$
(51)

Expanding out the phase space factors, we find several terms. The first one
involves two $\phi$ condensate distributions, corresponding to the process
$\phi_{c}\psi\leftrightarrow\phi_{c}\psi$. It is given by

$\hat{f}_{\phi}[p_{\phi}]\hat{f}_{\phi}[p_{3}]\left(f_{\psi}[p_{2}]-f_{\psi}[p_{4}]\right)~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}(\phi_{c}\psi\leftrightarrow\phi_{c}\psi).$
(52)

It is easily seen that this term vanishes since the delta functions from the
condensate distributions enforce $\vec{p}_{\phi}=\vec{p}_{3}=0$, and therefore
the energy conserving delta function enforces $E_{2}=E_{4}$, implying
$\left(f_{\psi}[p_{2}]-f_{\psi}[p_{4}]\right)=0$. Thus, Eq. (52) vanishes.

The second class of terms encountered in the phase space factor involve one
$\phi$ condensate and one thermal $\phi$, e.g., for the process
$\phi_{c}\psi\leftrightarrow\phi\psi$ we find

$\hat{f}_{\phi}[p_{\phi}]\left\\{f_{\phi}^{\rm
th}[p_{3}]\left(f_{\psi}[p_{2}]-f_{\psi}[p_{4}]\right)+f_{\psi}[p_{2}]\left(1-f_{\psi}[p_{4}]\right)\right\\}~{}~{}~{}~{}~{}~{}~{}~{}~{}(\phi_{c}\psi\leftrightarrow\phi\psi),$
(53)

and an analogous expression for the process
$\phi\psi\leftrightarrow\phi_{c}\psi$. Plugging in the explicit forms of the
equilibrium distributions and using energy conservation, we find

$\hat{f}_{\phi}[p_{\phi}]\left\\{f_{\phi}^{\rm
th}[p_{3}]\left(f_{\psi}[p_{2}]-f_{\psi}[p_{4}]\right)+f_{\psi}[p_{2}]\left(1-f_{\psi}[p_{4}]\right)\right\\}\propto
e^{(E_{2}+\mu_{\phi})/T}-e^{(E_{2}+m_{\phi})/T}.$ (54)

Notice that once $\mu_{\phi}=m_{\phi}$ is established, this factor also
vanishes. However the relevant collision term is always slower than the other
collision terms, and we therefore neglect it in the following.

Finally, there is a term corresponding to the scattering involving two thermal
$\phi$ particles, $\phi\psi\leftrightarrow\phi\psi$. This contribution
vanishes for our ansatz. Motivated by the equilibrium of this number-
conserving process, we consider kinetic equilibrium of the system to be
achieved when the condition $\mu_{\phi}[t]=m_{\phi}$ is satisfied.

##  3 DM stability with exotic fermions

In this section, we study the stability of the $\phi$ DM condensate against
decay and scattering processes in a simple model with a new exotic fermion
$\psi$. Let us first describe the approach to thermalization in this scenario.
Suppose that after the $\phi$ starts oscillating, in addition to the
condensate, there is plasma with number density $\mathcal{O}(m_{\phi}^{3})$
and a typical particle momentum of $p_{\rm typ}\sim m_{\phi}$. As will be
discussed below in Sec. 3.4, this is expected to arise at early times due to
broad parametric resonance production of $\psi$. It is known that for $p_{\rm
typ}\gtrsim m_{\phi}$ in such a system, particle decays lead to the momenta
being fully randomized on the timescale of
$\sim\left(\frac{y^{2}m_{\phi}}{16\pi}(\frac{m_{\phi}}{p_{\rm
typ}})^{n}\right)^{-1}$, with positive power $n$; see, e.g., Refs. [38, 39,
40, 41, 42] for investigations of such anisotropic plasma as well as further
discussion below in Sec. 4.1. Therefore, at the beginning stages, the
thermalization rate is parametrically $\sim\frac{y^{2}m_{\phi}}{16\pi}$, which
is model-independent. Thus the distribution functions for the plasma
components follow the quasi-equilibrium forms [42], justifying the ansatz
Eq.(12) and the theoretical setup described in the previous section, at least
for timescales longer than $(\frac{y^{2}m_{\phi}}{16\pi})^{-1}$. At this
timescale, the system consists of the nearly homogeneous and isotropic dark
radiation plasma and the DM condensate, which then evolve according to the
Boltzmann equations from the previous section.

###  3.1 Case with chirality suppression

We first consider the stability of the $\phi$ condensate, the emergence of the
thermalized $\phi$, $\psi$ dark radiation component, and its phenomenological
implications in the case that the self-annihilation of the $\phi$ condensate
is chirality suppressed. We therefore take the limit $m_{\psi}\to 0$ in this
subsection.666We note that qualitatively similar results will apply for the
case of a purely pseudoscalar coupling, in which case the $\phi$ annihilation
is velocity suppressed; see Appendix A for further details. In the presence of
the chirality suppression, the annihilation of $\phi$ is dominated by high
momentum modes with $p_{\phi}\gtrsim m_{\phi}.$ In particular, the
annihilation between two condensate particles, $\phi_{c}\phi_{c}\to\psi\psi$
is negligible, i.e. ${\cal
C}^{\phi_{c},\psi}_{\phi_{c}\phi_{c}\to\psi\psi}\rightarrow 0$. For
simplicity, we neglect the effect of the Hubble expansion in this subsection,
$H=0$.

The numerical solution to the Boltzmann equations, Eqs.(24), (25), (26) and
(30), are presented in Fig. 1. The initial conditions are taken to be
$n_{\phi}[0]/m_{\phi}^{3}=10^{8},T[0]=3/2m_{\phi},\mu_{\psi}[0]=\mu_{\phi}[0]=0$.
We fix $y=10^{-8}$ and neglect the Hubble expansion. We take account of the
annihilation and decay collision terms, ${\cal
C}_{\phi\phi_{c}\leftrightarrow\psi\psi}^{\phi,\phi_{c},\psi},{\cal
C}_{\phi/\phi_{c}\leftrightarrow\psi\psi}^{\phi,\phi_{c},\psi}$, respectively,
and neglect the others, which are sub-leading. For the relevant matrix element
expressions, see Appendix A. To study the late time evolution of the DM
condensate and the dark sector thermalization, we expand $\mu_{\phi}$ and
$\mu_{\psi}$ around $0$ in the collision terms of ${\cal
C}^{\phi_{c},\phi,\psi}_{\phi_{c}\phi\leftrightarrow\psi\psi}$ and ${\cal
C}^{\phi,\psi}_{\phi\leftrightarrow\psi\psi}$. This is done up to the third
power in both $\mu_{\phi}~{}{\rm and}~{}\mu_{\psi}$, which is a good
approximation in the regime $T\gtrsim\mu_{\rm\psi},\mu_{\phi}.$ This Taylor
expansion is also carried out for $n_{\phi}^{\rm th},\rho_{\phi}^{\rm
th},n_{\psi},~{}{\rm and}~{}\rho_{\psi}$. On the other hand, we use the full
expression for ${\cal C}^{\phi_{c},\psi}_{\phi_{c}\leftrightarrow\psi\psi}$ as
it can be evaluated exactly analytically.

Fig. 1: Stability and thermalization in the model of scalar DM $\phi$
interacting with exotic fermions $\psi$, with $m_{\psi}=0$. We take the
initial conditions
$\hat{n}_{\phi}[0]/m_{\phi}^{3}=10^{8},T[0]=3/2m_{\phi},\mu_{\psi}[0]=\mu_{\phi}[0]=0$
and fix $y=10^{-8}$. We also show relevant timescales corresponding to various
interactions in vertical green dashed lines (see the main text for details).

Let us now describe the evolution of the system using semi-analytic arguments.
Given the ansatz of the initial distribution, the decay of the $\phi$
condensate rapidly produces $\psi$ particles with momenta $m_{\phi}/2$. This
process is quickly saturated due to Pauli blocking, $\mu_{\psi}\rightarrow
m_{\phi}/2$, as indicated by the vertical green dashed line labeled
$\phi_{c}\leftrightarrow\psi\psi$ in Fig. 1. The system subsequently
approaches kinetic equilibrium defined at a rate of777Strictly speaking, the
collision term (39) with our ansatz gives the r.h.s to be
$-C_{\phi\leftrightarrow\psi\psi}^{\phi}/n_{\phi}^{\rm
th}\sim{m_{\phi}\Gamma_{\phi}}\frac{(2\mu_{\psi}-\mu_{\phi})}{T^{2}}$, with
$n_{\phi}^{\rm th}\sim T^{3}/\pi^{2}.$

$\Gamma_{\rm kin}\equiv\frac{y^{2}m^{2}_{\phi}}{16\pi T}\,.$ (55)

This represents the typical interaction rate of the splitting process of the
ambient plasma $\phi\leftrightarrow\psi\psi$. The typical momentum of the
fermion is taken to be of order $T$, and as such the Pauli blocking effect is
not important and is thus omitted in the rate estimate, Eq. (55). Although
realistically, the evolution for timescales smaller than $1/\Gamma_{\rm
kin}|_{T\sim m_{\phi}}$ depends in detail on the initial distribution, as we
have explained, when the timescale is longer than $1/\Gamma_{\rm kin}|_{T\sim
m_{\phi}}$, the evolution of bulk quantities should not depend much on the
initial distribution with the number density and energy density fixed. Thus
our estimates should be robust for timescales longer than $1/\Gamma_{\rm
kin}|_{T\sim m_{\phi}}$. As a consistency check we have verified numerically
that the late time evolution seen in Fig. 1 is obtained for a variety of
initial conditions for $\mu_{\psi},\mu_{\phi},T$.

Therefore, we expect that Fig. 1 provides a good description of the evolution
of the system to the right of the vertical green dashed line labeled
$\phi\leftrightarrow\psi\psi$. One observes from the figure that
$1/\Gamma_{\rm kin}|_{T\sim m_{\phi}}$, which corresponds to the conventional
DM lifetime in vacuum, does not represent the lifetime of the DM condensate in
our scenario, a consequence of Pauli blocking.

Next, the temperature of the plasma starts to grow at the timescale888More
precisely, the timescale for the onset of temperature growth, i.e., when
$\Delta T/T\sim{\cal O}(1)$, can be estimated from Eq.(59) and is in agreement
with the timescale Eq.(56) for $T\sim m_{\phi}$.

$\Delta t\sim\Gamma_{\rm
scat}^{-1}\equiv\left(\frac{y^{4}\hat{n}_{\phi}}{32\pi
Tm_{\phi}}\right)^{-1},$ (56)

corresponding to the inverse of the characteristic interaction rate of the
$\phi\phi_{c}\to\psi\psi$ scattering process, as indicated by the green
vertical line labeled $\phi\phi_{c}\to\psi\psi$ in Fig. 1. The reason the
temperature does not grow for $t\ll\Delta t$ is the number conservation
associated with the rapid $1\to 2$ process, corresponding to the chemical
equilibrium condition

$\mu_{\phi}=2\mu_{\psi}=m_{\phi}\,.$ (57)

At early times, while the $2\leftrightarrow 2$ process is inefficient, the
temperature does not grow because number conservation and energy conservation
fix the temperature. This is seen in Fig. 1 to the left of the
$\phi\phi_{c}\to\psi\psi$ line. We also note that $\Delta t\gg 1/\Gamma_{\rm
kin}$ is typically satisfied as a consequence of Eq.(9). Thus Eq.(57) is
satisfied in the whole range right to $\phi\leftrightarrow\psi\psi$ line
including the $t\gg\Delta t$ region.

When $t\gtrsim\Delta t$, the temperature starts to grow. This is because both
the decay and annihilation processes are efficient, and the particle number is
no longer conserved. This contradicts Eq.(57). Thus the condensate starts to
evaporate. Since the chemical potential is fixed to be (57) by the fastest
reaction, the evaporation is partially Pauli blocked. Indeed the evaporation
happens at the timescale of

$\Gamma_{\rm th}\sim\frac{1}{\Delta t}\frac{m_{\phi}}{T}.$ (58)

The suppression by $m_{\phi}/T$ is because the $2\to 2$ annihilation would be
blocked (or in kinetic equilibrium) if $\mu_{\psi}=\mu_{\phi}=m_{\phi}$ were
exact. Thus the collision term of the annihilation, i.e., Eq.(45), is
suppressed by
$\mathcal{O}((\mu_{\psi}-m_{\phi})/T,(\mu_{\phi}-m_{\phi})/T)\sim m_{\phi}/T$
from Eq.(57).

The plasma temperature growth can be studied by solving the following
Boltzmann equation,999More precisely, one can estimate the r.h.s from (2.2.2)
to be
$C_{\phi_{c}\phi\leftrightarrow\psi\psi}^{\phi_{c}}m_{\phi}\sim\frac{T^{3}}{\pi^{2}}[\frac{y^{4}\hat{n}_{\phi}}{32\pi
Tm_{\phi}}\frac{m_{\phi}}{T}]\log{\frac{T}{m_{\phi}}}\times m_{\phi}$.

${d\over dt}{\left(\rho_{\phi}^{\rm
th}+\rho_{\psi}\right)}=(g_{\phi}+7/8g_{\psi})\frac{\pi^{2}}{30}{d\over
dt}{T^{4}}\sim\Gamma_{\rm th}n_{\phi}^{\rm th}m_{\phi},$ (59)

with $n_{\phi}^{\rm th}\sim T^{3}/\pi^{2}.$ The last term is understood as the
number changing rate per unit volume via $\phi\phi_{c}\to\psi\psi$ reaction by
multiplying the transferred energy, $m_{\phi}$, per each collision. This can
be also checked from the behavior of the collision term of Eq.(2.2.2).
Integrating Eq.(59) until
$(g_{\phi}+7/8g_{\psi})\frac{\pi^{2}}{30}{T^{4}[t_{\rm
DM}]}\sim\hat{n}_{\phi}[0]m_{\phi}$ allows us to estimate the lifetime of
$\phi$ condensate:

$t_{\rm
DM}\sim\left[(g_{\phi}+\tfrac{7}{8}g_{\psi})\frac{\pi^{2}}{30}\right]^{1/4}\frac{128\pi^{3}}{3y^{4}\left(m_{\phi}\hat{n}_{\phi}\right)^{1/4}}\sim\frac{1300}{y^{4}\left(m_{\phi}\hat{n}_{\phi}\right)^{1/4}},$
(60)

where we have neglected the change of $\hat{n}_{\phi}$ and taken
$g_{\phi}=1,g_{\psi}=2.$ The timescale agrees well with the numerical
estimation for the benchmark in Fig. 1.

By requiring the lifetime to be longer than the age of the Universe, $13.8$
Gyr, and by assuming $\hat{n}_{\phi}\simeq\rho_{\rm DM}/m_{\phi}$ with
$\rho_{\rm DM}$ being the current (global) DM density, we obtain the condition
for $\phi$ as the DM,

$|y|\lesssim 2\times 10^{-7}.$ (61)

Strictly speaking, we need to take account of the redshift in the estimation,
which would not change the bound much. The inclusion of the Hubble expansion
will be discussed next in Sec. 3.2.

We can also estimate the dark radiation at the recombination era, $z_{\rm
rec}\sim 1100$, from

$\frac{(g_{\phi}+\tfrac{7}{8}g_{\psi})\pi^{2}}{30}T^{4}\sim\left(\frac{\Gamma_{\rm
th}m_{\phi}n_{\phi}^{\rm th}}{H}\sim\rho_{\rm
DM}\left.\frac{y^{4}T}{32\pi^{3}H}\right)\right|_{z=z_{\rm rec}}.$ (62)

Applying the constraint from Eq.(61), we obtain

$\Delta N_{\rm eff}[z_{\rm rec}]\lesssim 0.03\times\left(\frac{y}{2\times
10^{-7}}\right)^{16/3}.$ (63)

Here

$\Delta N_{\rm eff}[z]\equiv\frac{4}{7}\frac{30\rho_{\rm
DS}[z]}{\pi^{2}T_{\nu}[z]^{4}},$ (64)

with $T_{\nu}$ being the neutrino temperature.

Thus, the amount of dark radiation produced is relatively small, and the
thermal components of the dark sector are unlikely to impact the evolution of
the cosmic history given that $\phi$ lives much longer than the age of the
Universe. This is the case for $Q\ll 1$101010When $Q\lesssim 1$ but it is not
very small, there may be a similar effect as in the chirality unsuppressed
case $m_{\psi}\neq 0$ studied in Sec. 3.2 since there is an induced effective
fermion mass of $m_{\psi}^{\rm eff}\sim\sqrt{Q}m_{\phi}\sim
y\sqrt{\hat{n}_{\phi}/m_{\phi}}\lesssim m_{\phi}.$ One can find the effective
mass by taking the time average of $(y\phi[t])^{2}$, which is identified as
$(m_{\psi}^{\rm eff})^{2}$. This may provide an efficient dark radiation
production mechanism soon before the recombination.  and $m_{\psi}=0$ (or with
$m_{\psi}\neq 0$ and a $\phi\bar{\psi}\gamma^{5}\psi$ coupling, see Appendix
A). We will see in the next subsection that a qualitatively different
situation arises for $m_{\psi}\neq 0$, where the evolution of energy densities
can be quite different from that in the $\Lambda$CDM.

###  3.2 Case without chirality suppression

Let us now turn on the mass term of the fermion,

$\delta{\cal L}\supset-\frac{m_{\psi}}{2}\bar{\psi}\psi,$ (65)

with $m_{\psi}<m_{\phi}/2$ being the aforementioned mass parameter of the
fermion. The main difference from the massless case is that the $\phi$
condensate can self-annihilate, $\phi_{c}\phi_{c}\to\psi\psi$. Then we have
the collision term from Eq.(44),

${\cal
C}^{\phi_{c}}_{\phi_{c}\phi_{c}\leftrightarrow\psi\psi}\simeq-\frac{2y^{4}m^{2}_{\psi}}{\pi
m_{\phi}^{4}}\hat{n}_{\phi}[t]^{2}(1-2f_{\psi}[m_{\phi}]),$ (66)

where we have taken the leading term in the Taylor series in powers of
$m_{\psi}/m_{\phi}$.

The numerical result for the evolution of $\mu_{i},T,\hat{n}_{\phi}$,
including the collision term (44), is presented in Fig. 2 where we take
$m_{\psi}=m_{\phi}/50$. We do not expand ${\cal
C}^{\phi_{c},\psi}_{\phi_{c}\phi_{c}\leftrightarrow\psi\psi}$ in terms of
$T,\mu_{\phi},\mu_{\psi}$ as it may be evaluated exactly analytically. The
other parameters are the same as in Fig.1.

The time to the onset of temperature growth is shorter than in the chirality
suppressed case due to the enhanced condensate annihilation. However, the
duration of the temperature growth is longer than in the chirality suppressed
case. This is due to the Pauli blocking factor $1-2f_{\psi}[m_{\phi}]\sim
m_{\phi}/4T$ with $\mu_{\psi}=m_{\phi}/2$, which leads to a decrease in the
annihilation rate as $T$ increases. In contrast, in the previous chirality
suppressed case, $-{d\over dt}\log\hat{n}_{\phi}\propto T$ increases when $T$
increases. By including the contributions from both
$\phi_{c}\phi_{c}\to\psi\psi$ and $\phi_{c}\phi\to\psi\psi$ (see Eq.(60)), the
lifetime of the DM can be estimated as

$t_{\rm
DM}^{-1}\simeq\max{[0.2\frac{y^{4}m^{2}_{\psi}}{m_{\phi}^{4}}{\left(m_{\phi}\hat{n}_{\phi}\right)^{3/4}},~{}~{}0.0008y^{4}\left(m_{\phi}\hat{n}_{\phi}\right)^{1/4}]}$
(67)

where we have solved ${d\over dt}{\left(\rho_{\phi}^{\rm
th}+\rho_{\psi}\right)}\sim{\cal C}^{\phi_{c}}_{\phi_{c}\phi_{c}}\times
m_{\phi}$ for the first component and again defined the lifetime from the
condition $\left(\rho_{\phi}^{\rm
th}+\rho_{\psi}\right)=m_{\phi}\hat{n}_{\phi}[0]$ and neglected the decrease
of $\hat{n}_{\phi}$ in the estimation.

This form agrees well with the numerical solution to the Boltzmann equations
presented in Fig. 2. When $\phi_{c}\phi_{c}\to\psi\psi$ is the dominant
scattering process, the constraint, $t_{\rm DM}>13.8$ Gyr, turns out to be

$\frac{m_{\psi}}{m_{\phi}}\lesssim
0.8\left(\frac{10^{-9}}{y}\right)^{2}\left(\frac{m_{\phi}}{10^{-6}{\rm\,eV}}\right).$
(68)

A more stringent bound will be obtained in the following.

Fig. 2: Same as Fig.1 except that we take $m_{\psi}=m_{\phi}/50$.

Intriguingly, the slow thermalization provides a production mechanism for
abundant dark radiation in the early Universe. The dark radiation contribution
can lead to slight deviations from the $\Lambda$CDM predictions, while
maintaining the cosmic stability of the DM condensate. To study this
quantitatively, let us consider the Boltzmann equations including the effect
of the Hubble expansion. To simplify the discussion let us focus on timescales
longer than the time to reach kinetic equilibrium (right of the green dashed
line labeled $\phi\leftrightarrow\psi\psi$ in Fig. 2), i.e. $\Gamma_{\rm
kin}>H$, and let us first focus on the evolution when $T\gg m_{\phi}$ is
satisfied and the $\phi_{c}\phi_{c}\to\psi\psi$ contribution to scattering
dominates. Then Eq.(57) is satisfied.

The integrated Boltzmann equation for the $\phi$ condensate number density is
obtained from Eqs.(24) and (57) as

$\dot{\hat{n}}_{\phi}+3H\hat{n}_{\phi}\simeq-\frac{2y^{4}m_{\psi}^{2}\hat{n}^{2}_{\phi}}{\pi
m_{\phi}^{4}}(1-2f_{\psi}[m_{\phi}])\simeq-\frac{y^{4}m_{\psi}^{2}\hat{n}_{\phi}^{2}}{2\pi
m^{3}_{\phi}T},$ (69)

while that for the thermal component of dark sector energy density is found
from Eqs.(28), (29), and (57) to be

$\dot{\rho}_{\rm DS}+4H\rho_{\rm
DS}\simeq\frac{y^{4}m_{\psi}^{2}\hat{n}_{\phi}^{2}}{2\pi m_{\phi}^{2}T}.$ (70)

Note that (70) also follows straightforwardly from energy conservation.
Namely, the decrease of a $\phi$ particle in the condensate corresponds to the
increase of the radiation energy by $m_{\phi}$. Here $\rho_{\rm
DS}\simeq(g_{\phi}+\tfrac{7}{8}g_{\psi})\pi^{2}T^{4}/30=({11}/{4})\pi^{2}T^{4}/30$
is the dark radiation of relativistic $\phi$ and $\psi$ (with $T\gg
m_{\phi}$). Although $\phi_{c}\phi_{c}\to\psi\psi$ is the reaction that
transfers energy from the condensate to the fermions, the $\phi,\psi$
radiation readily shares the transferred energy since the system is in kinetic
equilibrium. The r.h.s of both equations are suppressed by $m_{\phi}/T$ thanks
to the Pauli exclusion principle.

By solving the equations one can estimate the final temperature of radiation
at different epochs. Before showing the numerical result, we will present an
approximate analytic calculation of the evolution of the dark radiation and
the resulting constraints. In our analytical estimate, we will neglect the
decrease of $\hat{n}_{\phi}a^{3}$ with time.

The radiation generated in a Hubble time $1/H$ can be estimated from Eq.(70)
and is given by $\delta\rho_{\rm
DS}a^{4}\sim\frac{y^{4}m_{\psi}^{2}\hat{n}_{\phi}^{2}}{2\pi
m_{\phi}^{2}TH}a^{4}$, with $a=1/(1+z)$ being the scale factor. By noting
$\hat{n}_{\phi}\propto a^{-3}$, we see that $a^{4}\rho_{\rm DS}$ is dominantly
generated at late times during the radiation- and matter-dominated eras.
During these eras, we can therefore obtain an estimate of the dark radiation
temperature from the equality $H\rho_{\rm
DS}\sim\frac{y^{4}m_{\psi}^{2}\hat{n}_{\phi}^{2}}{2\pi m_{\phi}^{2}T}$:

$T\sim\left.\left(\frac{15y^{4}m_{\psi}^{2}{\rho}_{\phi}^{2}}{(g_{\phi}+\frac{7}{8}g_{\psi})\pi^{3}Hm_{\phi}^{4}}\right)^{1/5}\right|_{z\gg
z_{\rm DE}},$ (71)

where $z_{\rm DE}\simeq 0.3$ represents the redshift at the dark energy-matter
equality. This temperature scales as $a^{-4/5}~{}{\rm and}~{}a^{-9/10}$ during
the radiation- and matter-dominated eras, respectively. Although the
production of the dark radiation is most efficient at late times during these
eras, $\rho_{\rm DS}/\rho_{\phi}$ is largest at early times in the radiation-
dominated era.111111The production of $\rho_{\rm DS}/\rho_{\phi}$ is dominant
at the early stage, say the onset of oscillation of $\phi$. However at that
epoch we will have $Q\gg 1$ and the discussion here breaks down (see Sec. 3.4
for the discussion in this regime.)  Also, the dark radiation scales
$\rho_{\rm DS}a^{4}\propto a^{4/5}~{}{\rm and}~{}a^{2/5}$ in the radiation-
and matter-dominated eras, respectively.

When dark energy is dominant, $z<z_{\rm DE}$, the produced $\delta\rho_{\rm
DS}a^{4}$ is smaller than that produced at $z=z_{\rm DE}$. Thus $\rho_{\rm
DS}$ for $z<z_{\rm DE}$ can be approximated from the redshift of $\rho_{\rm
DS}$ at $z=z_{\rm DE}$. Using Eq.(71) and accounting for the redshift, we
obtain the dark sector temperature at the present epoch:

$T[z=0]\sim\left.(1+z_{\rm
DE})^{-1}\left(\frac{15y^{4}m_{\psi}^{2}{\rho}_{\phi}^{2}}{(\frac{7}{8}g_{\psi}+g_{\phi})\pi^{3}Hm_{\phi}^{4}}\right)^{1/5}\right|_{z=z_{\rm
DE}}.$ (72)

By using the DM density scaled from today, $\rho_{\phi}=\rho_{\rm
DM}^{\rm\Lambda CDM}\equiv(1+z)^{3}\rho_{\rm DM}$ and the $\Lambda$CDM Hubble
parameter we obtain $T[z]$ and thus the dark radiation with a given set of $y$
and $m_{\phi}$. Interestingly, the resulting dark radiation does not depend on
the initial conditions.

In order to satisfy the bounds from the cosmic microwave background (CMB) and
baryon acoustic oscillation (BAO) observations [37, 16], we need $\Delta
N_{\rm eff}[z_{\rm rec}]\lesssim 1$.121212Here we use a relaxed bound rather
than $\Delta N_{\rm eff}\lesssim 0.3$[37] since our conclusions here are based
on order-of-magnitude estimates and since the dark radiation may not be free-
streaming as we will discuss (see also, e.g., Ref. [16] for discussion on the
relaxed bound for non-free-streaming dark radiation). This constraint turns
out to be

$\displaystyle y\sqrt{\frac{m_{\psi}}{m_{\phi}}}\lesssim 7\times
10^{-11}\sqrt{\frac{m_{\phi}}{10^{-6}{\rm\,eV}}}\text{ if }T\gg m_{\phi}\text{
at }z=z_{\rm rec}.$ (73)

On the other hand, $\Delta N_{\rm eff}[z_{\rm
rec}]=\mathcal{O}(0.01\text{-}1)$ may be tested in the future CMB and BAO
observations [43, 44, 45].

We solve the integrated Boltzmann equations (69) and (70) numerically for the
benchmark $y=10^{-9}$, $m_{\phi}=10^{-4}{\rm\,eV}$, and
$m_{\psi}=0.2m_{\phi}.$ The evolution of $\Delta N_{\rm eff}$ and
$1-\rho_{\phi}/\rho_{\rm DM}^{\rm\Lambda CDM}$, i.e., the deviation of the DM
density from the $\Lambda$CDM prediction as a function of the redshift, is
shown in Fig. 3 by the red dashed and solid lines, respectively. The initial
condition is taken to be $\rho_{\phi}[z_{\rm ini}]=\rho^{\Lambda\rm CDM}_{\rm
DM}$ and $T[z_{\rm ini}]=m_{\phi}$ at $z_{\rm ini}=10^{5}$.131313We take
$z=10^{5}$ as the initial condition because we would like to show the
evolution after $z=10^{4}-10^{5}$ at which the broad parametric resonance ends
(see the next subsection).  Here $\rho^{\Lambda\rm CDM}_{\rm DM}$ is the DM
density in the $\Lambda$CDM model. The temperature $T[z\ll z_{\rm ini}]$ does
not depend much on the initial condition given that the DM dominates. The
resulting $\Delta N_{\rm eff}[z=z_{\rm rec}]\simeq 0.33$ for this benchmark.
For comparison, we also consider the predictions of conventional decaying DM
models (i.e., negligible effects from Pauli blocking). We show the decaying DM
predictions for the same quantities by the black lines in Fig. 3, taking the
same initial conditions and fixing the decay rate to $\Gamma_{\rm
Decaying~{}DM}=0.1H_{0}$.141414One may consider exactly the same Lagrangian
Eq.(1) as a model of ordinary decaying DM. This is realized for large enough
$m_{\phi}$ so that the Pauli blocking is irrelevant and small enough $y$ so
that the decay rate in vacuum is suppressed. In other words, our findings may
be considered as a different parameter region of the models of the decaying DM
into fermions. We note that the decaying DM case is close to the bound on the
decay width given in [46, 47]. In this case, the DM density deviates by
$\mathcal{O}(10)\%$ from $\Lambda$CDM at the present epoch while the amount of
dark radiation is negligible (significant) at $z\sim z_{\rm rec}$ ($z\sim 0$).
In contrast, in our scenario, the DM density differs by $\mathcal{O}(1-10)\%$
from the $\Lambda$CDM model today while having substantial dark radiation near
recombination. Thus we expect that it is safer than the decaying DM with
$\Gamma_{\rm Decaying~{}DM}=0.1H_{0}.$ On the other hand, the production of
dark radiation in our scenario is slower than in the decaying DM scenario at
the earlier times while it ceases at late times $z\lesssim 1$. Thus the
numerical solution confirms the previous discussions. Indeed, an analytic
relation in our scenario can be obtained from the redshift dependence (see
Eq.(72)),

$\frac{\Delta N_{\rm eff}[z=0]}{\Delta N_{\rm eff}[z_{\rm
rec}]}=\left(\frac{H(z_{\rm rec})(1+z_{\rm DE})}{H(z_{\rm DE})(1+z_{\rm
rec})}\right)^{4/5}=\mathcal{O}(10).$ (74)

An analytic estimation gives this factor around 10, but numerically solving
the evolution equation gives around $20$. This is quite different from the
production of the dark radiation in the ordinary decaying DM scenario. This
slower DM “decay” and larger $\Delta N_{\rm eff}[z_{\rm rec}]$ may help to
alleviate the Hubble tensions, as will be discussed in Sec. 3.3.

Fig. 3:  The evolution of $\Delta N_{\rm eff}[z]$ and $1-\rho_{\phi}/\rho_{\rm
DM}^{\rm\Lambda CDM}$ (i.e., the deviation of the DM density from the
$\Lambda$CDM prediction) with redshift in our scenario (red dashed and solid
lines, respectively). The initial condition is set by
$\rho_{\phi}[z=10^{5}]=\rho^{\Lambda\rm CDM}_{\rm DM},T[z=10^{5}]=m_{\phi}$,
though the result is relatively insensitive to the dark sector temperature. We
have fixed $y=10^{-9}$, $m_{\phi}=10^{-4}{\rm\,eV}$, and
$m_{\psi}=m_{\phi}/5.$ We find $\Delta N_{\rm eff}[z_{\rm rec}]\sim 0.33.$ For
comparison, we also show the prediction for conventional decaying DM with
decay rate $\Gamma_{\rm Decaying~{}DM}=0.1H_{0}$ and the same initial
conditions (black lines).

When $T\gg m_{\phi}$ is not satisfied at some epoch, we obtain the stable DM
with the Fermi sea (4) discussed in Sec. 2.1 afterward. Whether $T\gg
m_{\phi}$ is satisfied or not depends not only on the parameters
$y,m_{\phi},m_{\psi}$ but also on $z$. Thus $T\gg m_{\phi}$ can be invalid
today but valid in the early Universe, implying that there is parameter region
$\Delta N_{\rm eff}[z_{\rm rec}]\simeq(0.01-1)$ in the early Universe while
$E_{F}=m_{\phi}/2$ today (right to the blue dashed line in Fig.4).

###  3.3 Predictions for the case without a chirality suppression

The slow but steady production of the dark radiation, a consequence of the
self-annihilation of the condensate DM in the chirality-unsuppressed scenario,
as well as the enhanced stability of the condensate, may have interesting
effects on cosmological evolution.

##### Preliminaries for alleviating Hubble tensions

If $\Delta N_{\rm eff}[z_{\rm rec}]\ll 1$, then our DM scenario reproduces the
ordinary $\Lambda$CDM model predictions. Nevertheless, when $\Delta N_{\rm
eff}[z_{\rm rec}]\sim\mathcal{O}(0.01-1)$ (corresponding to the purple shaded
region in Fig. 4), the model predicts a slight deviation from $\Lambda$CDM
which has the potential to alleviate the Hubble tensions. Let us qualitatively
describe the effects that should help to alleviate the Hubble tensions.

The analysis of the previous section and Fig.3 shows that our model realizes a
novel decaying DM plus dark radiation scenario in which the dark radiation has
self-interactions. It has been discussed that cold DM decaying into dark
radiation can partially alleviate the $H_{0}$ and $\sigma_{8}$-tensions [48,
49, 50]. However, this scenario faces severe constraints arising from the late
time changes to the DM density, which affect structure formation, enhance the
late integrated Sachs Wolfe effect, and suppress the gravitational lensing
effects [46, 47, 51, 52]. These effects leave observable signatures in BAO,
Sunyaev-Zel’dovich cluster counts, and CMB lensing. In contrast, the DM
condensate in our scenario has a suppressed “decay” rate at late times,
implying minimal impact on these observables such that those bounds should be
evaded (see Fig.3).

Furthermore, it is known that when $\Delta N_{\rm eff}[z_{\rm
rec}]=\mathcal{O}(0.1\text{-}1)$ around the time of recombination, the
$H_{0},M_{B}$-tensions can be relaxed [14, 15, 12, 16]. This is because the
extra radiation decreases the sound horizon size at the decoupling of CMB
photons, which decreases the angular distance and thus alleviates the
tensions. However, the amount of free-streaming dark radiation is severely
constrained by the CMB data since it increases the Silk-damping and neutrino
drag effects. On the other hand, if the dark radiation is self-interacting, it
clusters more than the free-streaming radiation on small scales, counteracting
the enhancement of the Silk damping effect. As a consequence, the bound on the
amount of dark radiation is relaxed if it is non-free-streaming [16]. The
efficient self-interactions of the dark radiation in our scenario may inhibit
free streaming.

Thus, for sizable couplings $y\sim{\cal O}(10^{-9})$, our scenario, which
features a slow dark matter condensate evaporation along with substantial dark
radiation with self interactions, contains the right elements to help to ease
the Hubble tensions. A rigorous assessment of this possibility is quite non-
trivial, requiring the use of a Boltzmann code to evolve the cosmological
perturbations in order to make accurate predictions for cosmological
observables, as well as a careful comparison of these predictions with various
observational datasets, and goes beyond the scope of this work.

##### Suppression of small scale structure

Another interesting prediction of our scenario is the absence of structures
with large DM densities, whose formation is inhibited by DM evaporation. For
instance, our scenario disfavors cuspy DM halo profiles for couplings that not
too small. One can naively estimate the maximal energy density of the DM by
assuming that the dark radiation is homogeneous and has a temperature of
$T[0]$. Then the maximal energy density persisting during the timescale
$\Delta t$ satisfies

$\frac{y^{4}m_{\psi}^{2}\rho_{\phi}^{\rm
max}}{2\pi\max{[T[0],m_{\phi}]}m^{4}_{\phi}}\times\Delta t\sim 1.$ (75)

By taking $\Delta t=1/H_{0}$ this gives

$\frac{\rho_{\phi}^{\rm max}}{{\rm\,GeV}{\rm
cm}^{-3}}\sim\max{[0.2\left(\frac{m_{\phi}}{1\,{\rm
meV}}\right)^{8/5}\left(\frac{10^{-9}}{y}\right)^{16/5}\frac{m_{\phi}^{8/5}}{m_{\psi}^{8/5}},1\left(\frac{m_{\phi}}{1\,{\rm
meV}}\right)^{3}\left(\frac{10^{-9}}{y}\right)^{4}\frac{m_{\phi}^{2}}{m_{\psi}^{2}}]}.$
(76)

where we have used Eq.(72). Note that $m_{\psi}/m_{\phi}\ll 1/2$ gives an
enhancement of the l.h.s. This naive estimation implies that the scenario may
not lead to cuspy DM halo profiles, at least for large enough couplings.

Given $\rho_{\phi}^{\rm max}$, we can estimate an upper bound on the DM mass
in a finite sphere with a radius $r$,

$M^{\rm max}[r]=2.9\times 10^{7}M_{\odot}\left(\frac{r}{300\rm
pc}\right)^{3}\frac{\rho_{\phi}^{\rm max}}{10{\rm\,GeV}{\rm cm}^{-3}},$ (77)

where $M_{\odot}$ denotes the solar mass. For instance, a conventional dwarf
spheroidal galaxy has a mass $M_{300}\sim 10^{7}M_{\odot}$ within $r<300$ pc
[53]. This may set a bound

$\rho^{\rm max}_{\phi}\gtrsim(1\text{-}10){\rm\,GeV}{\rm cm}^{-3}.$ (78)

While this is satisfied in most of the parameter regions, it may probe part of
the region when the dark radiation is abundant. A more robust bound, beyond
our preliminary estimate here, would require further detailed studies,
including simulations of structure formation in our scenario with relatively
large couplings. This could be particularly interesting given the anticipated
future improvements in precision measurements of the masses of dwarf
spheroidal galaxies as well as the (sub)halo mass function, e.g., using the
Vera Rubin Observatory [17].

###  3.4 Parametric resonance and DM production

So far, we have focused on the perturbative regime defined by Eq.(9), where we
can use the Boltzmann equations to describe the system. On the contrary, when

$Q\gg 1,$ (79)

$\psi$ acquires a large oscillating mass $M_{\rm eff}(t)\simeq y\phi{(t)}$,
the typical size of which is

$\bar{M}_{\rm eff}\sim y\phi_{\rm amp}=Q^{1/2}m_{\phi}\gg m_{\phi},$ (80)

where $\phi_{\rm amp}$ is the amplitude of the oscillating $\phi$ field. In
this case, a perturbative description breaks down, and we cannot naively use
the Boltzmann equation formalism. A similar system has been studied in the
context of preheating after inflation [54, 55, 23, 56, 24]. Here we follow the
analytical discussion for preheating given in [55] to quantify and describe
the dynamics.

One can consider the frequency of the fermion particle, $w_{k}$ of momentum
$k$ at time $t$

$w_{k}=\sqrt{k^{2}+y^{2}\phi[t]^{2}}.$ (81)

The particle picture of the fermion can be justified when the adiabatic
condition is satisfied

$\dot{w}_{k}\lesssim w_{k}^{2}.$ (82)

When this condition is violated with large $|\dot{\phi}|$ , efficient $\psi$
particle production occurs. Taking $\phi=\phi_{\rm amp}\cos[m_{\phi}t],$ this
violation can only occur at $\phi\sim 0$, or more precisely
$|\phi|\lesssim(m_{\phi}\phi_{\rm amp}/y)^{1/2}$. The momentum $k$ should also
be small for (82) to be violated, i.e.,

$k\lesssim k_{*}\simeq\sqrt{ym_{\phi}\phi_{\rm amp}}=Q^{1/4}m_{\phi}.$ (83)

As a result, a Fermi sphere with all modes less than $k_{*}$ is produced. The
Fermi momentum is thus $k_{F}\sim k_{*}$ [23, 24].

Since $k_{F}\ll\bar{M}_{\rm eff}$, the fermion is non-relativistic when the
particle picture is justified with (82). During most of the evolution (and
thus the time average of) the energy density is given by

$\rho_{\psi}\sim\bar{M}_{\rm
eff}n_{\psi}=\frac{g_{\psi}}{2\pi^{2}}Q^{5/4}m_{\phi}^{4}=\frac{g_{\psi}}{2\pi^{2}}y^{2}Q^{1/4}\rho_{\phi},$
(84)

where we used $n_{\psi}\simeq\frac{g_{\psi}}{2\pi^{2}}k_{F}^{3}.$ We get the
scaling of the energy density of Fermi sea as $\rho_{\psi}\propto a^{-15/4}.$

However, the simplification that $\phi$ is described as a classical
oscillating field should be treated with some care. For example, it is known
that if there is a large quartic interaction of $\phi$, ${\cal
L}\supset-\lambda\phi^{4}$, satisfying $y^{2}/\lambda\ll 1$, the system enters
into a turbulence regime. This is found by performing a lattice simulation
with a 2PI effective action in specific models [57, 58]. However, this should
not apply to our scenario if a small quartic coupling of typical radiative
size $\lambda\sim y^{4}/16\pi^{2}$ is present. The detailed numerical study of
the $Q\gg 1$ regime for the quadratic potential is beyond our scope. In the
following, we make use of the results of Refs. [23, 56, 24] with a quadratic
potential of $\phi$ in which back-reactions, such as the re-scattering effect,
are neglected.

$Q\propto a^{-3}$ in the expanding Universe. Thus $Q$ may be larger than unity
in the early Universe even if $Q\ll 1$ at the later time.

Fig. 4: Parameter region of DM stabilized by exotic fermions in
$m_{\phi}\text{-}y$ plane with $m_{\phi}=5m_{\psi}$. Above the red dashed
lines, broad parametric resonance is important at $z\leq z_{\rm rec}$ and
$z=0$ from bottom to top, respectively. Above the blue solid line, $\Delta
N_{\rm eff}(z_{\rm rec})\geq 1$ and is disfavored. In the purple region,
$\Delta N_{\rm eff}=0.01\text{-}1$ is predicted at $z=z_{\rm rec}$. Below the
blue dotted line, we get the Fermi sea with Eq.(4) in the present Universe.
Above it, the thermal component including the Fermi sea has a temperature
$T[0]$ much larger than $m_{\phi}/2$. The green region above the green solid
line is disfavored since the comoving DM density may change significantly at
$z<10^{6}$ (See Eq.(88)). We also show isocontours of the maximal local DM
density in the present Universe, which is a specific prediction of our
scenario, roughly estimated according to Eq.(78).

Due to the suppressed fermion energy density $\rho_{\psi}$ compared to that of
DM $\rho_{\phi}$ in the regime $Q\gg 1$, we can consider the conventional
misalignment mechanism [34, 35, 36] for the production of the $\phi$
abundance. That is, $\phi$ starts to oscillate around the potential minimum
when $m_{\phi}\sim H$, and subsequently, the coherently oscillating scalar
field forms a pressureless, nonrelatvistic fluid and thus serves as cold DM.
The onset of oscillation is at the photon temperature

$T^{\rm osc}_{\gamma}\sim
300{\rm\,GeV}\left(\frac{m_{\phi}}{10^{-4}{\rm\,eV}}\right)^{1/2}\left(\frac{100}{g_{\star}}\right)^{1/4},$
(85)

from $H\sim m_{\phi}.$ From Eq.(84) and $\rho_{\phi}=\rho_{\rm DM}^{\rm\Lambda
CDM}$, we get

$\rho_{\psi}\sim
10^{-10}\left(\frac{y}{10^{-9}}\right)^{5/2}\left(\frac{10^{-4}{\rm\,eV}}{m_{\phi}}\right)\left(\frac{T_{\gamma}}{300{\rm\,GeV}}\right)^{3/4}\rho_{\phi},$
(86)

and consequently $\rho_{\psi}<\rho_{\phi}$ is satisfied with
$T_{\gamma}<T_{\gamma}^{\rm osc}.$ Noting that $T_{\rm osc}\propto
m_{\phi}^{1/2},$ we always have $\rho_{\psi}<\rho_{\phi}$ assuming the
constraints Eqs.(6) and (61) are satisfied. The fermion energy density
$\rho_{\psi}$ is suppressed until $Q$ decreases to $\lesssim 1$ due to the
expansion of the Universe, and following this we should use the expression of
(71). This happens below

$T_{\gamma}^{\rm pert}\sim
5{\rm\,eV}\left(\frac{m_{\phi}}{10^{-4}{\rm\,eV}}\right)^{4/3}\left(\frac{10^{-9}}{y}\right)^{2/3}.$
(87)

Thus our previous discussion with Eq.(71) in the regime $z\lesssim
T_{\gamma}^{\rm pert}/T_{\gamma}^{\rm today}-1$ with $T_{\gamma}^{\rm today}$
being the photon temperature today, can be consistent with the misalignment
mechanism.151515 It is conceivable that a full numerical study of the broad
parametric resonance regime could reveal further nonperturbative effects which
cause the DM condensate to evaporate. We emphasize that our scenario should
still work in the perturbative regime $Q<1$ together with a mechanism for
lowering the initial redshift $z_{\rm ini}$ and dynamically delaying the onset
of the $\phi$ oscillations (see, e.g., Refs. [59, 60, 61, 62]). In particular,
a thermal potential may delay the onset of the $\phi$ oscillations, see, e.g.,
[63, 64, 65].

We must require that the DM does not evaporate too quickly after the end of
the broad parametric resonance due to the perturbative interactions. To this
end, we employ the condition

$\rho_{\rm DS}[z_{Q=1}]<\rho_{\rm DM}[z_{Q=1}]$ (88)

with $z_{Q=1}$ defined by $Q[z_{Q=1}]=1$, i.e., $T_{\gamma}=T_{\gamma}^{\rm
pert}$. Here $\rho_{\rm DS}$ is obtained using Eq.(71) evaluated in the $Q\ll
1$ region. We impose this constraint if $z_{Q=1}<10^{6}$. If this condition is
violated, too much of the DM energy density may be transferred into radiation
during structure formation, cf.[66, 67].161616If $z_{Q=1}\gg 10^{6}$, it
should be safe as long as the remnant at lower red shift is consistent with
the DM abundance. Namely in the earlier time, the DM abundance decreases like
in the WIMP case. This condition excludes the green region in Fig.4. However,
we emphasize that this constraint should be taken with caution. When $Q\sim 1$
the system is non-perturbative, i.e., many $\leftrightarrow$ many processes
may be as important as $1\leftrightarrow 2,2\leftrightarrow 2$ processes and
there is also an oscillating effective mass of $\psi$, which are not accounted
for in our Boltzmann equation approach. These nonperturbative effects could in
principle be studied using lattice methods.

###  3.5 Cosmology and summary

Let us summarize one of the possible cosmological histories for the scenario
in which DM is stabilized by exotic BSM fermions. When
$T_{\gamma}>T_{\gamma}^{\rm osc}$ (Eq.(85)), $\phi$ is frozen at some field
value, which may be generated stochastically due to inflationary quantum
fluctuations. As the temperature drops, $\phi$ eventually starts to oscillate
around its potential minimum when $T_{\gamma}\sim T_{\gamma}^{\rm osc}$, i.e.,
it is produced via the misalignment mechanism. At this moment, $Q\gg 1$ and
the system is in the regime of broad parametric resonance, resulting in the
production of a Fermi sea with energy density (84). Still, the energy density
of the $\phi$ condensate dominates over $\rho_{\psi}$, and $\phi$ is stable
[23, 56, 24]. As the universe expands, $Q\propto a^{-3}$ redshifts, and when
$Q\sim 1$, $k_{*}\sim m_{\phi}$ and $\rho_{\psi}\sim y^{2}\rho_{\phi}\sim
m_{\phi}^{4}$. Afterward, the perturbative description accounting for the
leading $1\leftrightarrow 2$ decay/inverse decay and $2\leftrightarrow 2$
scattering processes using Boltzman equations is justified. The results of the
Boltzmann equation analysis can be found in Figs. 1 and 2 with the initial
condition of $T\sim m_{\phi}$. The evolution of the system on timescales
longer than $1/\Gamma_{\rm kin}$ are relatively insensitive to the initial
conditions, and therefore, the imprints of the initial broad parametric
resonance phase are largely erased.

The parameter space for DM stabilized by exotic fermions in the $m_{\phi}-y$
plane is displayed in Fig. 4, where we take $m_{\psi}/m_{\phi}=1/5.$ In the
green region the DM density changes by an $\mathcal{O}(1)$ amount soon after
the broad parametric resonance ends, which takes place at $z<10^{6}$, and may
thus be excluded. The gray region above the blue solid line is disfavored
since $\Delta N_{\rm eff}(z_{\rm rec})>1.$ In the purple region we predict
slowly produced interacting dark radiation (see Fig. 3) with $\Delta N_{\rm
eff}(z_{\rm rec})=0.01-1$, which may potentially alleviate the Hubble
tensions, suppress the formation of DM structures with large densities, and
can be probed in the future CMB and BAO experiments [43, 44, 45]. Above the
lower red dashed line the $Q$ parameter at the recombination era, $z\sim
z_{\rm rec}$, is greater than unity and the dark radiation production is
suppressed. Above the upper red dashed line the system is in the broad
parametric resonance regime in the present Universe.

Above the blue dotted lines, 2-to-2 scattering processes are important and
$T\gg m_{\phi}/2$ in the present Universe. Below the blue dotted line the
naive discussion in Sec. 2.1, which accounts only for 1-to-2 decay/inverse
decays processes (including Pauli blocking) but neglects scattering, applies
in the present Universe (but not necessarily the earlier Universe) since the
$\phi_{c}$ annihilation is not efficient. In this case, the system today has a
Fermi sea component, Eq.(4), due to the Pauli blocked decay of $\phi$.

Finally, in the orange region, the DM has a longer lifetime than the age of
the Universe in vacuum, and thus corresponds to the standard decaying DM
scenario.

##  4 DM stabilized by C$\nu$B

So far, we have investigated DM stabilization due to Pauli-blocking in a
minimal model containing an exotic light BSM fermion. It is natural to
consider the possibility that the light fermion is instead a neutrino, either
a left-handed neutrino or a right-handed neutrino, i.e., $\psi=\nu_{L}$ or
$\nu_{R}$. The $\psi$ in the previous section can be straightforwardly
identified with $\nu_{R}$ without changing our previous conclusions provided
$\nu_{R}$ has a suppressed mass and not too large mixing with the ordinary
neutrinos. The more interesting case is $\psi$ as $\nu_{L}$, which we consider
in the following.

The interaction Lagrangian is given as 171717Neutrino couplings to scalar dark
matter have been investigated in various other contexts, see, e.g., Refs. [68,
69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79].

${\cal
L}\supset-\frac{1}{2}y_{ij}\phi(\bar{\nu}^{c}_{L})_{i}(\nu_{L})_{j}+{\rm
h.c.},$ (89)

where $i={1,2,3}$ is the index in the mass eigenbasis, e.g., in the normal
ordering case, the lightest, heavier, and heaviest neutrinos, respectively. We
assume the normal hierarchy, i.e., normal ordering with a small lightest
neutrino mass,

$m_{\nu_{1}}<m_{\phi}/2,$ (90)

so that DM is stabilized as a consequence of the Pauli-blocking of its decays
to the lightest neutrino. We also assume a generic flavor structure,
$y_{ij}\sim y$, and as such we omit the flavor index in the following.

There are several aspects of this scenario that are different from the minimal
model discussed in Sec. 3. First, there is the presence of the C$\nu$B from
the conventional Big Bang cosmology that contributes $N_{\rm eff}\sim 3$,
which was confirmed in the CMB and BAO observations [1]. Since the C$\nu$B at
early times should not deviate significantly from the standard cosmological
predictions, in the following, we approximate that $T_{\nu}$ of the plasma is
the usual neutrino temperature, $T^{\rm SM}_{\nu}$, i.e., $T_{\nu}[z>z_{\rm
rec}]\sim T_{\nu}^{\rm SM}[z>z_{\rm rec}]$, before the CMB era in most of the
discussion. A slight modification to the standard cosmology will be also
discussed. On the other hand, in the later epochs the C$\nu$B may be altered
significantly without violating the CMB and BAO constraints. As we will see,
the change in the late time C$\nu$B spectrum is a natural prediction of our
scenario.

Second, the extra interaction of the C$\nu$B with $\phi$ may inhibit the free
streaming of neutrinos, in contradiction with inferences from CMB and BAO data
[16]. This leads to a new constraint on the coupling $y$.

Lastly, the heavier neutrinos are expected to decay during the age of the
Universe, even if the coupling of $y$ is not very large. This outcome is
natural since we focus on the parameter region where the $\phi$ proper
lifetime (in vacuum) is shorter than the age of the Universe, and given the
assumption of a generic flavor structure, a coupling of the same order leads
to an even faster decay of the heavier neutrinos ($\nu_{3}~{}{\rm
or}~{}\nu_{2}$). Eventually, the $\phi_{c}$ decay is Pauli-blocked in our
universe, but the decay of $\nu_{3,2}$ is not.

We discuss these three points in more detail in the following.

###  4.1 C$\nu$B free-streaming and DM stabilized by neutrino sea

Under our assumption of a generic flavor-violating structure in the
$\phi$-neutrino couplings, $y_{ij}\sim y$, a heavier neutrino in the rest
frame can decay into $\phi$ and a lighter neutrino with a decay width

$\tau_{\nu_{2,3}}^{-1}\equiv\Gamma_{\nu_{2,3}}\sim\frac{y^{2}m_{{\nu}_{2,3}}}{16\pi}$
(91)

where $m_{\nu_{i}}$ is the neutrino mass of $(\nu_{L})_{i}$, and we neglect
the mass of the decay products. When the lifetime is shorter than the age of
the Universe, the heavier neutrinos are absent in the present C$\nu$B, which
then consists of only the lightest neutrino.

Such decays may be constrained as they can hinder neutrino free-streaming via
the process $\nu_{3}\leftrightarrow\phi+\nu_{1,2}.$ These reactions should not
be too frequent around the recombination era. To study free-streaming it is
important to check the evolution of the anisotropic component of the plasma.
For the decay of a relativistic particle of energy $E_{\nu_{3}}$ and mass
$m_{\nu_{3}}$ in the cosmic frame, daughter particles have momenta in the
collinear angle $m_{\nu_{3}}/E_{\nu_{3}}$. An inverse decay by the
annihilation of a particle with the ambient plasma produces heavier particles
also in a collinear angle $\sim m_{\nu_{3}}/E_{\nu_{3}}$. Thus the initial
momentum diffuses to the whole momentum space after order
$(m_{\nu_{3}}/E_{\nu_{3}})^{-2}$ reactions, similar to a random walk. Thus, it
has been argued that the timescale for momentum transport is further
suppressed by the diffusion factor together with the boost factor
$(m_{\nu_{3}}/T)^{2}\times m_{\nu_{3}}/T$ [38, 39, 40, 41]. Moreover, recently
it was found that there is an additional suppression numerically (and
understood analytically under several assumptions) by $(m_{\nu_{3}}/T)^{2}$
due to a cancellation of the collision terms of the relevant order in a
derived evolution equation for the total anisotropic stress [42]. Here we
adopt the bound [42] (see also [38, 39, 40, 41])181818This is the case with a
generic flavor violation for the $\phi$ coupling. In the case of the Majoron
or $B-L$ Higgs boson, the neutrino decay of $\nu_{2,3}$ is suppressed. But
there is a free-steaming constraint coming from the reaction
$\phi\leftrightarrow\nu_{1}\nu_{1}$. The bound can be similarly obtained [42]
by using the timescale
$\Gamma^{-1}_{\phi\to\nu_{1}\nu_{1}}\left(\frac{m_{\phi}}{T}\right)^{-5}$.

$y\lesssim 10^{-9}\left(\frac{m_{\nu_{3}}}{0.05{\rm\,eV}}\right)^{-3}.$ (92)

If this bound is satisfied, then two-to-two processes involving neutrinos and
thermalized $\phi$ particles will also be suppressed [42].

Let us argue that in our scenario, the condition (92) is sufficient to ensure
neutrino free-streaming, including also the effects of the specific scattering
reactions involving the DM condensate
$\phi\phi_{c}\leftrightarrow\nu_{1}\nu_{1}$ and
$\phi_{c}\nu_{1}\leftrightarrow\phi\nu_{1}$. The reaction rate (56) is not the
appropriate one for damping anisotropic stress since the outgoing particles
are collinear. By including the diffusion factor $(E_{\rm cm}/T)^{2}\sim
m_{\phi}/T$, we estimate the transport rate to be
$\Gamma_{T}\sim\frac{m_{\phi}}{T}\Gamma_{\rm
scat}\left(\frac{m_{\phi}}{T}\right)^{n-1}$, where $n-1\geq 0$ is a number
that represents the possible further cancellation in the evolution equation of
the total anisotropic stress including the potential Pauli blocking effect. We
find that $\Gamma_{T}\sim
y^{4}\frac{T^{2}}{m_{\phi}}\left(\frac{m_{\phi}}{T}\right)^{n-1}\frac{\hat{n}_{\phi}m_{\phi}}{T^{4}}$.
Since $\hat{n}_{\phi}m_{\phi}\sim T^{4}$ around matter-radiation equality and
thus recombination, and $n>0$, the transport rate is always smaller than the
Hubble rate in this parameter region.

It should be noted that if Eq.(92) is satisfied, the usual interaction rate
$\frac{y^{2}}{4\pi}m_{\nu_{3}}^{2}/T_{\nu}$ can only become faster than the
Hubble rate much after the BBN era. Therefore, the entropy of the C$\nu$B
shortly after decoupling around the BBN epoch is the same as that of the SM
case.

###  4.2 Spectrum and composition of C${\nu}$B

Our scenario predicts a non-trivial spectrum and composition of the C$\nu$B as
a consequence of neutrino decays and neutrino production from the $\phi$
condensate after recombination. The spectrum of the C$\nu$B is, in principle,
measurable at future C$\nu$B detection experiments, and thus provides a novel
probe of our model. In the following, we will consider two scenarios. In the
first part (Sec. 4.2.1), we consider the case in which the DM annihilation is
chirality suppressed. In the second part (Sec. 4.2.2), we consider the
chirality unsuppressed case. The second case is similar in most respects to
the first case, but when the coupling is large we may potentially alleviate
the Hubble tension. As we will see, both scenarios typically feature non-
standard components of the C$\nu$B which may be probed in the future.

####  4.2.1 Case of chirality suppression

In the chirality suppressed case, $m_{\nu_{1}}/m_{\phi}\ll 1$, the lifetime of
$\phi$ can be estimated from Eq.(59) to be $\sim(y^{4}T)^{-1}\sim 10^{7}\times
13.8\,{\rm Gyr}\left(\frac{y}{10^{-9}}\right)^{-4}\frac{T^{\rm
SM}_{\nu}[0]}{T}.$

There are two sources of the C$\nu$B in the present Universe. One is the
standard primordial source of neutrinos which originates from the hot Big
Bang. We stress that in our scenario the spectrum, however, may be quite
different from the standard one because the additional interaction involving
$\phi$ leads to the decays of neutrinos. The effect also depends on the size
of the coupling $y$, thus we will examine two regimes for its strength in the
following. The other source of the C$\nu$B is the $\phi$ condensate decays,
which will be discussed lastly, and this contribution should be added to the
primoridal one.

##### C$\nu$B component from primordial neutrinos at weak coupling

If $y\lesssim 10^{-14}$ the heaviest neutrino decay happens at
$T_{\nu}\lesssim m_{\nu_{3}}\approx 0.05{\rm\,eV}.$ Thus $\nu_{3}$ becomes
non-relativistic first and then decays via

$\nu_{3}\to\phi({\rm boosted})+\nu_{1,2}.$ (93)

In this case, the decays happen out of equilibrium. Interestingly, since the
Pauli blocking effect can be neglected for the boosted $\phi$ decay,191919In a
small portion of the parameter space, there is also the possibility that the
further decay of $\nu_{2},\phi(\rm boosted)$ does not happen within the age of
the Universe due to the time dilation effect, which we do not account for in
our order-of-magnitude estimate.

$\nu_{2}\to\phi(\rm boosted)+\nu_{1},~{}~{}\phi(\rm boosted)\to 2\nu_{1}$ (94)

can happen subsequently. Here $\nu_{2}$ may be either the product of a
$\nu_{3}$ decay or a relic that was produced around BBN. Assuming $\nu_{3}$
decays to $\nu_{2}\phi$ and $\nu_{1}\phi$ with the same rate, the total
comoving number density of the neutrinos increases from
$(n_{\nu,1}+n_{\nu,2}+n_{\nu,3})a^{3}\to(8n_{\nu,1})a^{3}$. Some event rates
for the neutrino interaction with electrons are proportional to $|(U_{\rm
PMNS})_{ei}|^{2}=\\{0.7,0.3,0.02\\}$, with $U_{\rm PMNS}$ being the PMNS
matrix. In particular, this is the case for the tritium capture rate in the
PTOLEMY project [19, 20, 21]. When the heavier neutrinos all decay into the
lightest one, the event rate can be enhanced by a factor202020 In contrast,
for the inverted hierarchy case with a negligible lightest neutrino mass, the
event rate is suppressed by $0.02\times 8/(0.7+0.3+0.02)\sim 10\%$.

$\frac{0.7\times(1+3+0.5(3+5))}{(0.7+0.3+0.02)}\sim 6~{}~{}~{}~{}~{}~{}({\rm
typical~{}}E_{\nu}=\mathcal{O}(0.01){\rm\,eV})$ (95)

compared to the SM case. Here we used the fact that the event rate does not
depend much on the energy of the neutrino. Here and hereafter, the bracket
after the formula for the enhancement of the C$\nu$B number density denotes
the typical energy of the component. Thus, the precise measurement of the
C$\nu$B may provide a probe of our scenario. The expected energy resolution of
the PTOLEMY experiment is $0.01-0.1{\rm\,eV}$ [19], which favors energetic
neutrinos. In our case, the part of the numerator $3+0.5(3+5)$ represents the
neutrinos produced in the cascade processes of Eqs.(93) and (94), which are
boosted in comparision to the primordial component. These have energies as
large as $\sim(1/4-1/2)m_{\nu_{2,3}}=\mathcal{O}(0.01){\rm\,eV}$ if
$\nu_{2},\nu_{3}$ decays around the present Universe such that we can neglect
the effects of redshift.

##### C$\nu$B component from primordial neutrinos at moderate coupling

When $y\gtrsim 10^{-14}$ the decay and inverse decay happen thermally at
$T_{\nu}\gtrsim m_{\nu_{3}}$. All neutrinos reach quasi-thermal equilibrium
with $\phi$ well after BBN. The comoving entropy stored in the neutrino-$\phi$
system is conserved and is equal to the SM value carried by the three
neutrinos following their decoupling around the BBN era at which point the
$\phi$ interaction is frozen. Later, after the decay/decoupling of all the
heavier neutrinos, we have the relation

$(1+7/4)\frac{2\pi^{2}}{45}(T_{\nu}[0])^{3}\simeq
3\times(7/4)\frac{2\pi^{2}}{45}(T_{\nu}^{\rm SM}[0])^{3}$ (96)

from comoving entropy conservation. Here we have assumed that the mass of
$\phi$ is smaller than the temperature. We will discuss the case of heavier
$\phi$ later. As a result, the final comoving neutrino number density of the
stable lightest $\nu_{1}$ increases by

$\left(\frac{T_{\nu}[0]}{T^{\rm SM}_{\nu}[0]}\right)^{3}\simeq\frac{21}{11},$
(97)

and thus the detection rate is enhanced by a factor of

$\frac{0.7\times(21/11)}{(0.7+0.3+0.02)}\sim 1.3~{}~{}~{}~{}~{}~{}({{\rm
typical~{}}E_{\nu}=\mathcal{O}(T^{\rm SM}_{\nu}[0])}).$ (98)

In addition, for $10^{-9}\gtrsim y\gtrsim 10^{-13},$ the thermalization of
$\phi-\nu$ system happens before recombination, meaning that the energy
density of the neutrino sector is slightly different from $\Lambda$CDM. Thus,
it may be probed via measurements of $\Delta N_{\rm eff}$. To perform the
estimate we can again make use of the comoving entropy density conservation.
Near the recombination era, all the three neutrinos as well as $\phi$ are
relativistic, implying that the relativistic entropy degrees of freedom in
this sector are $g_{s\nu}\simeq 3\times 2\times\frac{7}{8}+1$ compared to the
SM case $3\times 2\times\frac{7}{8}$. This leads to

$(1+3\times 7/4)(T_{\nu}[z_{\rm rec}])^{3}\frac{2\pi^{2}}{45}\simeq(3\times
7/4)(T^{\rm SM}_{\nu}[z_{\rm rec}])^{3}\frac{2\pi^{2}}{45}.$ (99)

We obtain

$\Delta N_{\rm eff}[z_{\rm rec}]\approx-0.17.$ (100)

A deviation of this size can be tested in future CMB and BAO experiments [43,
44, 45]. We note that the decrease in $\Delta N_{\rm eff}$ during the CMB era
does not change the prediction (98), consistent with comoving entropy
conservation.

##### C$\nu$B components from DM decay

Since the DM stability in the present Universe is a consequence of the Pauli-
blocked decay to the lightest neutrinos, and, at relatively small coupling
strength the $\phi$ annihilation is negligible, we have an extra component of
the C$\nu$B corresponding to the Fermi sea produced by the decay of the
condensate. This component is characterized by the Fermi momentum of
$m_{\phi}/2$ and number density $n_{\nu_{1}}\sim(m_{\phi}/2)^{3}/6\pi^{2}$,
and therefore it dominates over the other components when $m_{\phi}\gg
T_{\nu}^{\rm SM}[0]$. In this case, the interaction rate of an electron with
C$\nu$B is enhanced by a factor

$\sim 10^{3}\left(\frac{m_{\phi}}{5\,\rm
meV}\right)^{3}~{}~{}~{}~{}~{}~{}~{}({\rm
typical~{}}E_{\nu}=\mathcal{O}(m_{\phi})),$ (101)

which is the ratio of the present-day neutrino number density to the standard
one. Although the energy of these neutrinos is relatively small compared to
the component from the heavy neutrino decays due to the bound Eq.(6), there
can be several orders of magnitudes enhancement in the event rate in the
future C$\nu$B experiment. Detecting this component is the smoking gun signal
of our scenario.

####  4.2.2 Case without chirality suppression at large coupling

For the scenario in which $\phi_{c}$ self-annihilation is not chirality
suppressed, differences from Sec. 4.2.1 only appear for relatively large
Yukawa couplings, $y\sim 10^{-9}$. In this case, the energy density of the
C$\nu$B may be slightly enhanced around the recombination era due to the slow
$\phi_{c}$ evaporation. Following the discussion in Sec. 3.2, the C$\nu$B
temperature is close to the temperature (71) at recombination. At later times,
the temperature follows the attractor solution, e.g., at $z=0$ we have
Eq.(74). Thus when the $\nu_{1}$ radiation from $\phi_{c}$ evaporation is
produced in comparable amounts to the primoridal C$\nu$B, the evolution of the
neutrino radiation follows the discussion in Sec. 3.2, resulting in a slight
enhancement of the total C$\nu$B at $z=0$, provided the DM is light enough. In
particular, near recombination we can obtain $N_{\rm eff}[z_{\rm rec}]\sim
3-4$ while today, analogously to Eq.(74), we get $N_{\rm eff}[0]\sim 10N_{\rm
eff}[z_{\rm rec}].$ Thus we expect an enhancement of the C$\nu$B detection
rate by as much as a factor of

$\sim\frac{0.7\times(4\times 10)^{3/4}}{0.7+0.3+0.02}\sim
10~{}~{}~{}~{}~{}~{}~{}({\rm typical~{}}E_{\nu}=\mathcal{O}(T_{\nu}^{\rm
SM}[0])).$ (102)

Such an enhancement of this kind occurs in the narrow blue-shaded region in
Fig. 5. However, when $m_{\phi}>10^{-4}{\rm\,eV}$, the degenerate neutrino sea
from the direct decay of $\phi_{c}$ is more important, as in the chirality
suppressed case.

For $y\ll 10^{-9}$, the additional neutrino radiation from $\phi_{c}$
evaporation is negligible. Still, along with the primordial component, we
predict the novel C$\nu$B components originating from neutrino decay as well
as DM decay, similarly to the chirality-suppressed case discussed in the
previous subsection.

Similarly to Sec. 3.3, there is an interesting possibility that our scenario,
with relatively sizable couplings $y\sim 10^{-9}$ (see Fig. 4), may help to
alleviate the Hubble tension. In this case, additional neutrino radiation is
produced near recombination from $\phi_{c}$ evaporation. The produced
radiation with the standard neutrinos has a shorter free-streaming length than
usual. It was pointed out that a majoron (scalar) coupled to neutrinos may
relax the Hubble tension by introducing extra neutrino radiation $\Delta
N_{\rm eff}[z_{\rm rec}]\sim 0.5$ [80, 81]. In our scenario, extra radiation
corresponding to $\Delta N_{\rm eff}[z_{\rm rec}]\sim 0.5+0.17\sim 0.7$ can be
provided naturally by $\phi_{c}$ evaporation, where the additional amount
$~{}0.17$ is needed to compensate the negative value of Eq.(100) caused by the
late thermalization of $\phi$ particles.

In summary, in the case that the DM condensate annihilation to neutrinos is
not chirality suppressed, we have a special parameter region (corresponding to
the narrow blue shaded range in Fig. 5) where the $\phi_{c}$ evaporation
produces extra neutrino radiation, which may potentially alleviate the Hubble
tension and be probed in C$\nu$B detection experiments. When the DM mass is
larger than about $0.1\,$meV, the Hubble tension may still potentially be
addressed by this additional neutrino radiation, however this component is
less important in the context of the C$\nu$B detection compared to the
degenerate neutrino sea originating from the direct decay of DM in the present
Universe. In both cases, we also predict an extra C$\nu$B component from the
decays of the heavier neutrinos. Below the blue region, we have the
essentially same predictions as in the chirality suppressed case.

####  4.2.3 Discussions on perturbativity

The parameter space of DM stabilized by the C$\nu$B in the $m_{\phi}-y$ plane
is presented in Fig. 5. We take $m_{\nu_{1}}/m_{\phi}=1/5.$ Above the red
solid line, the neutrino is massive during the BBN epoch due to the large
$\phi$ amplitude, which will be discussed at the end of this section. We note
that the boundaries of these regions may vary to some extent depending on the
component values of the coupling matrix Above the red dotted line, the
$Q[z_{\rm rec}]$ parameter at the recombination becomes larger than unity and
the Boltzmann equation at this era becomes non-perturbative. Analyzing this
regime may require a different approach such as e.g., lattice methods.

Fig. 5: The prediction for the present universe C$\nu$B with the DM stabilized
by C$\nu$B in $m_{\phi}-y$ plane. We take $m_{\nu_{1}}=1/5m_{\phi}.$ In the
top blue-shaded region, the extra radiation from the DM evaporation increases
the C$\nu$B-electron interaction rate by at most $\sim 10$ (see the discussion
around Eq.(102)). Below it, the $\rm C\nu B$ is thermalized, but we do not
have the extra radiation but an enhancement of $\sim 1.3$ (see the discussion
around Eq.(98)). In both cases, the C$\nu$B is composed of a plasma of
$\nu_{1}$ and $\phi$ with the chemical potential of $\mu_{1}=m_{\phi}/2.$ In
the right top, the C$\nu$B is composed of the “degenerate” $\nu_{1}$ from the
DM decay. We can have several orders of magnitude of enhancement of the
detection rate (see the discussion around Eq.(101)). In the bottom right
region, the C$\nu$B is composed of boosted $\nu_{1}$ from $\nu_{3},\nu_{2}$
and boosted $\phi$ decays. The enhancement can be at most a factor of $6$.
Above the red region above the red line, the averaged neutrino masses become
heavier than the temperature at the BBN. The other colored regions are the
same as Fig. 4.

###  4.3 Production mechanism of $\phi$ stabilized by C$\nu$B.

Before concluding this section, we discuss a constraint that is relevant if
the DM condensate is produced through the misalignment mechanism. In this
scenario, the onset of the oscillations takes place at the temperature given
by Eq.(85). This is much higher than the BBN temperature and may even be
higher than the electroweak phase transition temperature. The effective
neutrino mass $\bar{M}_{\rm eff}$ during these early times is given by
Eq.(80). In particular, one should avoid non-relativistic neutrinos during the
CMB and BBN epochs, leading to the condition $\bar{M}_{\rm eff}\lesssim
T_{\gamma}$ at $T_{\gamma}=1{\rm\,MeV}$. This translates to the bound

$|y|\lesssim 4\times 10^{-9}\left(\frac{m_{\phi}}{10^{-5}{\rm\,eV}}\right).$
(103)

Furthermore, we note that at such early times and high temperatures that the
free energy density of neutrinos sources a finite temperature contribution to
the $\phi$ potential. Generically, the thermal mass of $\phi$ in this early
epoch will dominate over the bare mass, such that $\phi$ no longer redshifts
as matter during this era. However, once the temperature drops sufficiently so
that the bare scalar mass again dominates, $\phi$ will behave as matter. It is
also possible to imagine that the scalar potential is suitably modified so
that $\phi$ oscillations do not begin until some later epoch following BBN. In
this sense, the bound (103) is model-dependent.

For energies and temperatures larger than the weak scale, we must find a
reasonable UV origin of the Lagrangian (89), consistent with the gauge
symmetries of the SM. The UV model should be consistent with neutrino mass
generation, and furthermore, the thermalization of $\phi$ via the interactions
among $\phi$, neutrinos (or left-handed lepton), and gauge bosons must be
avoided. Consistent renormalizable models utilizing the misalignment mechanism
can be constructed by introducing right-handed neutrinos, and we present one
such example model in Appendix C). An alternative possibility to evade these
problems is to consider other mechanisms for the production of the DM $\phi$,
such as inflaton decays [26, 27].

##  5 Generalizations to other models and signatures

Here we discuss some possible generalizations to other models that may realize
the mechanism of DM stability from Pauli blocking, as well as some additional
potential signatures that may be present depending on how the dark sector
couples to the SM degrees of freedom.

##### Other models with DM stability from Pauli blocking

In the minimal model presented in Sec. 3 we have assumed that $\psi$ is
stable. However, this is mainly for simplicity. Our scenario should work
equally well even if $\psi$ decays to lighter states. This is because the
$\psi$ decay products must include at least one fermion, and thus the decay of
$\psi$ will be Pauli-blocked as well.

Our mechanism for DM stabilization can also be applied to more generic
interactions than the ones considered in this work. For instance, instead of
Eq.(1) one can consider the interaction $\phi\bar{\psi}\Psi$, where $\psi$ is
a light or massless fermion while $\Psi$ is a heavy fermion with mass
$M_{\Psi}$. The DM $\phi$ decays, filling the $\psi$ and $\Psi$ Fermi spheres
such that the decay is eventually Pauli blocked. Interestingly, the
corresponding Fermi momenta can be made arbitrarily small by taking the mass
splitting $m_{\phi}-M_{\Psi}$ to be small. In this limit, the energy
transferred from the DM to the fermions $\psi$ and $\Psi$ also becomes small.
In comparision to the minimal model of Sec. 3, the advantage of this scenario
is that one can have a relatively heavy $\phi$ at the expense of a tuning
between $m_{\phi}$ and $M_{\Psi}.$212121Such tuning can be related with the
radiatively induced neutrino mass in certain UV models [82, 83, 84].  In this
case, however, we must require that the couplings of $\Psi$ to other lighter
particles are small so that $\phi$ decays via $\Psi$-mediated interactions are
ineffective.

Instead of a primary two body decay that stabilizes DM, we may have three body
decay channels involving fermionic decay products. These decays can also
become Pauli blocked and DM can be stabilized. As an interesting example,
consider the Lagrangian $\phi\bar{\psi}\sigma^{\mu\nu}\psi F_{\mu\nu}$, where
$F_{\mu\nu}$ is the photon field-strength. In this scenario, the DM
continuously produces photons as it replenishes the Fermi sea during the
expansion of the Universe. On the other hand, a boosted $\phi$ (e.g., a dark
radiation component) can decay into the fermions and a photon with
substantially higher energies. The two distinct sources of photons can thus
potentially provide an interesting signature of the mechanism. The coupling to
photons may lead to further constraints from early universe cosmology, which
should studied carefully.

In this work we have considered a scalar DM candidate, but our mechanism
should also work for the case of vector boson DM coupled to fermions. Most
discussions do not change qualitatively. A difference is that the annihilation
of the DM is not chirality suppressed even in the massless limit of the decay
product.

##### Indirect signatures from other couplings

So far in this scenario, we have neglected the possible couplings of
$\phi,\psi$ to SM particles, but they can be present on general grounds and
may lead to additional indirect signatures of the dark sector. Here we mention
two possibilities for illustration.

The symmetries of the theory, Eq.(1), allow higher dimensional interactions of
the form $\frac{1}{M}HL\phi\psi\supset\frac{v}{M}\nu\phi\psi,$ with
$H,v,~{}{\rm and}~{}L$ being the SM Higgs doublet, its vacuum expectation
value, and a lepton doublet, respectively. Such couplings cause the C$\nu$B to
decay to the lighter dark sector particles. Interestingly, the lifetime of the
heaviest SM neutrino can be below the age of the Universe when the scale
$M\lesssim M_{\rm pl}$.222222More generally, this fact shows that C$\nu$B
search may provide a nice opportunity for testing a light dark sector with a
fermion and a scalar, even if direct couplings are only generated by
gravitational effects. For instance, a lighter dark SM could be probed if the
ordinary neutrino decays to a dark neutrino and dark Higgs boson. Precisely
measuring the C$\nu$B and finding the absence of the heaviest neutrino may
probe this scenario.

Another interesting possibility arises if $\phi$ couples to a pair of photons,
albeit with a decay rate that is highly suppressed. This scenario is natural
if the $\phi$-photon coupling arises from a loop effect, such that it is
highly suppressed compared to the $\phi$-fermion (Sec. 3) or $\phi$-neutrino
coupling (Sec. 4). In the case of DM stabilized by the C$\nu$B, we have two
sources of photon signals in the present Universe. One is from the non-
relativistic DM rare decay to photons, with photon energy of half the DM mass
(neglecting the redshift). The other is from the decay of $\phi$ produced from
the decay of a neutrino. For the second source, the resulting typical energy
of photon is around half of the mass difference of the neutrinos between the
initial and final states. Interestingly, the latter energy does not depend
much on the mass of $\phi$ if it is small enough. This process can be probed
by C$\nu$B experiments searching for neutrino radiative decay such as COBAND
[85].

##  6 Conclusions and outlook

In this work we have explored the novel implications of the hypothesis that DM
stability is a consequence of the Pauli exclusion principle. We studied a
simple model containing a sub-eV scalar $\phi$ DM candidate with a Yukawa
coupling to a lighter exotic fermion $\psi$. Even though the scalar DM proper
lifetime in vacuum, $1/\Gamma(\phi\rightarrow\psi\psi)$, is much shorter than
the age of the Universe, the available fermionic $\psi$ states produced in the
decay are rapidly filled and the decay is Pauli blocked. We have considered
not only the impact of decay/inverse decay processes on the evolution of the
scalar DM condensate, but also the relevance of scattering processes and
parametric resonance effects, finding a substantial parameter space in which
the effective lifetime of the DM condensate is much longer than the age of the
Universe. In particular, we have demonstrated that for relatively large
couplings, scattering processes can be rapid enough to produce a substantial
quasi-thermal dark sector component behaving as (self-interacting) dark
radiation, which can potentially help to relax the Hubble tensions. Future
precision CMB and BAO observations will provide additional tests of the
scenario. Furthermore, a cutoff in the densities of DM structures is expected,
motivating future work on simulations of structure formation, which may help
furnish an additional robust probe of the scenario.

Perhaps the most natural candidate for the light fermion is the SM left-handed
neutrinos. In addition to the predictions of the exotic fermion scenario
summarized above, this scenario predicts significant modifications to the
C$\nu$B at late times, both in terms of its composition and energy spectrum.
There are several possible effects that may distort the C$\nu$B, including the
cascade decays of the heavier SM neutrinos as well as DM evaporation and/or
decay. In particular, in a large fraction of the viable parameter space, the
interaction rate of the C$\nu$B with an electron is significantly enhanced.
This opens up the possibility of testing our scenario through precise
measurements of the C$\nu$B with experiments such as PTOLEMY.

Looking ahead, there is wide scope for further investigations of the scenarios
proposed here and the general DM stability mechanism. In particular, it would
be very interesting to explore the detailed predictions for precision
cosmological and astrophysical observables through, e.g., the study of
cosmological perturbations and/or structure formation, which will help to more
accurately delineate the constraints and future prospects for the various
probes discussed in this paper. Likewise, for the scenario in which DM is
stabilized by neutrinos, it would be valuable to carry out more detailed
phenomenological investigations of the properties of the present-day C$\nu$B
and its detection prospects, including the exploration of more general flavor
patterns of the $\phi$-neutrino couplings. Finally, it would be worthwhile to
explore other models realizing the mechanism of DM stability by Pauli blocking
and their phenomenology.

## Acknowledgements

W.Y. thanks the Department of Physics and Astronomy of the University of
Pittsburgh for its hospitality when this work was initiated in 2019. B.B. is
supported in part by the U.S. Department of Energy under grant No. DE-
SC0007914 and in part by PITT PACC. W.Y. was supported by JSPS KAKENHI Grant
Nos. 16H06490, 20H05851, 21K20364, 22K14029, 22H01215, 23K22486 and by NRF
Strategic Research Program NRF-2017R1E1A1A01072736.

##  Appendix A Generic formulas for decay and scattering

In the main text, we have focused on real scalar as the DM interacting with a
Majorana fermion. In this appendix, we provide the relevant formulas for decay
and scattering reactions for both scalar and pseudoscalar couplings. To this
end, we consider the Lagrangian interaction

${\cal
L}\supset-\frac{1}{2}\phi\,\overline{\psi}(y_{S}+y_{P}i\gamma^{5})\psi,$ (104)

where $y_{S}$ and $y_{P}$ are both real. Here $\psi$ describes a four-
component Majorana fermion.

###  A.1 Decay $\phi\rightarrow\psi\psi$

The partial decay width for $\phi\rightarrow\psi\psi$ is

$\Gamma_{\phi\to\psi\psi}=\frac{1}{2!}\frac{m_{\phi}}{8\pi}\left(y_{S}^{2}\beta_{\psi}^{3}+y_{P}^{2}\beta_{\psi}\right),~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\beta_{\psi}=\left(1-\frac{4m_{\psi}^{2}}{m_{\phi}^{2}}\right)^{1/2}.$
(105)

We note that for the scalar coupling, the decay width exhibits a $P$-wave
suppression, while for the pseudoscalar coupling, the decay proceeds in the
$S$-wave.

###  A.2 Annihilation $\phi\phi\rightarrow\psi\psi$

Here we provide some details for the annihilation process
$\phi\phi\rightarrow\psi\psi$ for the case of scalar and pseudoscalar
couplings.

####  A.2.1 Scalar coupling

Consider the Lagrangian interaction, i.e. $y_{S}=y,y_{P}=0$ in Eq.(104),

${\cal L}\supset-\frac{y}{2}\phi\overline{\psi}\psi.$ (106)

The amplitude for the annihilation process is

${\cal
M}=-y^{2}\,\overline{u}(p_{3})\left[\frac{\not{\\!p}_{3}-\not{\\!p}_{1}+m_{\psi}}{t-m_{\psi}^{2}}+\frac{\not{\\!p}_{1}-\not{\\!p}_{4}+m_{\psi}}{u-m_{\psi}^{2}}\right]v(p_{4}),$
(107)

where we have introduced the Mandelstam variables with the standard
definitions, $s=(p_{1}+p_{2})^{2}$, $t=(p_{1}-p_{3})^{2}$,
$u=(p_{1}-p_{4})^{2}$. The spin summed squared amplitude is

$|\overline{\cal M}|^{2}=y^{4}\times
2\left[\frac{F_{tt}}{(t-m_{\psi}^{2})^{2}}+\frac{F_{uu}}{(u-m_{\psi}^{2})^{2}}+\frac{F_{tu}+F_{ut}}{(t-m_{\psi}^{2})(u-m_{\psi}^{2})}\right],$
(108)

where

$\displaystyle F_{tt}$
$\displaystyle=t\,u-m_{\psi}^{2}\,(9\,t+u)-7\,m_{\psi}^{4}+8\,m_{\psi}^{2}\,m_{\phi}^{2}-m_{\phi}^{4},$
$\displaystyle F_{uu}$
$\displaystyle=t\,u-m_{\psi}^{2}\,(9\,u+t)-7\,m_{\psi}^{4}+8\,m_{\psi}^{2}\,m_{\phi}^{2}-m_{\phi}^{4},$
(109) $\displaystyle F_{tu}=F_{ut}$
$\displaystyle=-t\,u-3\,m_{\psi}^{2}\,(t+u)-9\,m_{\psi}^{4}+8\,m_{\psi}^{2}\,m_{\phi}^{2}+m_{\phi}^{4}.$

The differential cross-section can be written as

$\displaystyle
d\sigma=\frac{1}{2!}\frac{1}{16\pi\beta_{\phi}^{2}s^{2}}|\overline{\cal
M}|^{2}dt,$ (110)

where $\beta_{i}=(1-4m_{i}^{2}/s)^{1/2}$.

The formula for the integrated cross-section in full generality is complicated
and we will not write it here. However, it is interesting to consider the
nonrelativistic limit and compute the quantity $\sigma v_{\rm rel}$, where
$v_{\rm rel}=|v_{1}-v_{2}|=2\beta_{\phi}$ is the relative velocity of the
$\phi$ particles in the CM frame. In the non-relativistic limit we obtain

$\sigma v_{\rm rel}\simeq\frac{2y^{4}m_{\psi}^{2}}{\pi
m_{\phi}^{4}}\left(1-\frac{m_{\psi}^{2}}{m_{\phi}^{2}}\right)^{3/2}.$ (111)

We see that the annihilation proceeds through the $S$-wave, but exhibits a
chirality suppression. The implications of the two cases are explored further
in Section Sec. 3.

####  A.2.2 Pseudocalar coupling

Consider the Lagrangian interaction, i.e. $y_{S}=0,y_{P}=y$ in Eq.(104),

${\cal L}\supset-i\frac{y}{2}\phi\overline{\psi}\gamma^{5}\psi.$ (112)

The amplitude of the annihilation process is

${\cal
M}=y^{2}\,\overline{u}(p_{3})\,\gamma^{5}\,\left[\frac{\not{\\!p}_{3}-\not{\\!p}_{1}+m_{\psi}}{t-m_{\psi}^{2}}+\frac{\not{\\!p}_{1}-\not{\\!p}_{4}+m_{\psi}}{u-m_{\psi}^{2}}\right]\,\gamma^{5}\,v(p_{4}),$
(113)

We see the appearance of the $\gamma^{5}$ now. The spin-averaged squared
amplitude is

$|\overline{\cal M}|^{2}=y^{4}\times
2\left[(t-m_{\psi}^{2})(u-m_{\psi}^{2})-m_{\phi}^{4}\right]\left[\frac{1}{t-m_{\psi}^{2}}-\frac{1}{u-m_{\psi}^{2}}\right]^{2}.$
(114)

The cross-section is obtained from the general formula listed in the previous
subsection. In the non-relativistic limit we obtain

$\sigma v_{\rm rel}\simeq v_{\rm rel}^{4}\,\frac{y^{4}}{60\pi
m_{\phi}^{2}}\left(1+\frac{3m_{\psi}^{2}}{2m_{\phi}^{2}}\right)\left(1-\frac{m_{\psi}^{2}}{m_{\phi}^{2}}\right)^{3/2}.$
(115)

The annihilation proceeds through the $D$ wave but does not exhibit a
chirality suppression. While we have note studied the pseudoscalar coupling in
detail in this work, we expect this scenario to give qualitatively similar
features to the scalar coupling in the chirality suppressed cased studied in
Sec. 3.1.

##  Appendix B Limits on fermion energy in collision integrals

###  B.1 Decay $\phi(p_{\phi})\rightarrow\psi(p_{1})\psi(p_{2})$

In this collision term (39) we encounter an integral over the fermion energy
$E_{1}$. The limits are determined from the energy conserving delta function
for $\cos\theta=\pm 1$, leading to the condition

$0=E_{\phi}-E_{1}-\sqrt{m_{\psi}^{2}+p_{\phi}^{2}+p_{1}^{2}\pm
2p_{\phi}p_{1}}.$ (116)

Solving this equation, we obtain

$\displaystyle p_{1}^{\pm}$ $\displaystyle=\frac{1}{2}(E_{\phi}\beta_{\psi}\pm
p_{\phi}),$ (117) $\displaystyle E_{1}^{\pm}$
$\displaystyle=\frac{1}{2}(E_{\phi}\pm\beta_{\psi}p_{\phi}),$ (118)

where $\beta_{\psi}$ is given in Eq. (105).

###  B.2 Annihilation
$\phi(p_{\phi})\phi(p_{2})\rightarrow\psi(p_{3})\psi(p_{4})$

In this collision term (2.2.2) we encounter an integral over the fermion
energy $E_{3}$. As above, the limits of integration follow from the energy
conserving delta function when $\cos\theta=\pm 1$, leading to the equation

$0=m_{\phi}+E_{2}-E_{3}-\sqrt{m_{\psi}^{2}+p_{2}^{2}+p_{3}^{2}\pm
2p_{2}p_{3}}.$ (119)

Solving this equation, we obtain

$\displaystyle p_{3}^{\pm}$
$\displaystyle=\frac{1}{2}[(E_{2}+m_{\phi})\kappa_{\psi}\pm p_{2}],$ (120)
$\displaystyle E_{3}^{\pm}$
$\displaystyle=\frac{1}{2}(E_{2}+m_{\phi}\pm\kappa_{\psi}\,p_{2}),$ (121)

where in this case

$\kappa_{\psi}=\left(1-\frac{2m_{\psi}^{2}}{m_{\phi}(E_{2}+m_{\phi})}\right)^{1/2}.$
(122)

##  Appendix C Renormalizable model for DM stabilized by C$\nu$B

In this appendix we present a renormalizable UV model realizing the scenario
discussed in Sec. 4. Along with the scalar $\phi$ which serves as DM, the
theory contains a Dirac fermion mediator $N$. The Lagrangian shares some
features with the Dirac-neutrino portal model [86, 87, 88],

${\cal L}_{\rm
int}=\frac{y_{1}}{2}\phi\bar{{N}}^{c}_{L}N_{L}+\frac{y_{2}}{2}\phi\bar{{N}}^{c}_{R}N_{R}+M_{N}\bar{N}_{R}N_{L}+y_{R}\bar{N}_{R}HL+{\rm
h.c.}~{}$ (123)

The model enjoys an approximate lepton $Z_{4}$ parity where
$(\phi,N_{R},N_{L},L,H)$ are charged as $(2,1,1,1,0)$. By integrating out the
Dirac right-handed neutrinos one obtains the effective Lagrangian,

${\cal L}_{\rm eff}\supset\frac{y}{2}\phi\bar{\nu}_{L}^{c}\nu_{L}+{\rm h.c.}$
(124)

In small $y_{2}$ regime, we obtain $y=y_{1}y_{R}^{2}v^{2}/2M_{N}^{2}$ at
energies much below the electroweak scale $v$. Notice that the neutrino
Majorana mass term $\tfrac{1}{2}m_{\nu}\overline{\nu}^{c}_{L}\nu_{L}$ is
forbidden by the $Z_{4}$ symmetry. To generate neutrino masses, one can
introduce a small lepton parity violating term.

Setting aside possible model-dependent UV contributions, the radiative
corrections to the $\phi$ mass from the new degrees of freedom in this model
are given by

$\delta m_{\phi}^{2}\sim\frac{y_{1}^{2}M_{N}^{2}}{16\pi^{2}}.$ (125)

For example, taking $y_{R}v\sim M_{N}$ such that $y\sim y_{1}$, we obtain

$\delta
m_{\phi}^{2}\sim\left(10^{-4}{\rm\,eV}\right)^{2}\left(\frac{y}{10^{-14}}\right)^{2}\left(\frac{|M_{N}|}{100{\rm\,GeV}}\right)^{2}.$
(126)

In the early Universe, especially in the symmetric phase, the thermalization
rate of $\phi$ particle is estimated as

$\Gamma_{\rm th}\sim\frac{y_{1}^{2}y_{R}^{2}}{4\pi}T.$ (127)

This is much smaller than the expansion rate of the Universe for $T\sim
100{\rm\,GeV}~{}{\rm and}~{}y_{1}\sim y$ in the parameter region in Fig.5. On
the other hand, the induced $\phi$ thermal mass may be larger than the bare
mass depending on the parameter region. In this case, $\phi$ condensate
behaves as radiation at the onset of oscillation, rather than matter because
the effective mass scales with $a^{-1}$. However, our overall conclusions do
not change since we can require that at later times $\phi$ accounts for the DM
once the thermal mass becomes smaller than the bare mass due to the redshift.

## References

  * [1] Planck collaboration, N. Aghanim et al., _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6 [1807.06209].
  * [2] A. G. Riess, L. Macri, S. Casertano, H. Lampeitl, H. C. Ferguson, A. V. Filippenko et al., _A 3% Solution: Determination of the Hubble Constant with the Hubble Space Telescope and Wide Field Camera 3_ , _Astrophys. J._ 730 (2011) 119 [1103.2976].
  * [3] V. Bonvin et al., _H0LiCOW – V. New COSMOGRAIL time delays of HE 0435 $-$1223: $H_{0}$ to 3.8 per cent precision from strong lensing in a flat $\Lambda$CDM model_, _Mon. Not. Roy. Astron. Soc._ 465 (2017) 4914 [1607.01790].
  * [4] A. G. Riess et al., _A 2.4% Determination of the Local Value of the Hubble Constant_ , _Astrophys. J._ 826 (2016) 56 [1604.01424].
  * [5] A. G. Riess et al., _Milky Way Cepheid Standards for Measuring Cosmic Distances and Application to Gaia DR2: Implications for the Hubble Constant_ , _Astrophys. J._ 861 (2018) 126 [1804.10655].
  * [6] S. Birrer et al., _H0LiCOW - IX. Cosmographic analysis of the doubly imaged quasar SDSS 1206+4332 and a new measurement of the Hubble constant_ , _Mon. Not. Roy. Astron. Soc._ 484 (2019) 4726 [1809.01274].
  * [7] A. G. Riess et al., _A Comprehensive Measurement of the Local Value of the Hubble Constant with 1 km/s/Mpc Uncertainty from the Hubble Space Telescope and the SH0ES Team_ , 2112.04510.
  * [8] W. L. Freedman, B. F. Madore, T. Hoyt, I. S. Jang, R. Beaton, M. G. Lee et al., _Calibration of the Tip of the Red Giant Branch (TRGB)_ , 2002.01550.
  * [9] N. Khetan et al., _A new measurement of the Hubble constant using Type Ia supernovae calibrated with surface brightness fluctuations_ , _Astron. Astrophys._ 647 (2021) A72 [2008.07754].
  * [10] W. L. Freedman, _Cosmology at a Crossroads_ , _Nature Astron._ 1 (2017) 0121 [1706.02739].
  * [11] L. Verde, T. Treu and A. G. Riess, _Tensions between the Early and the Late Universe_ , _Nature Astron._ 3 (2019) 891 [1907.10625].
  * [12] E. Di Valentino et al., _Snowmass2021 - Letter of interest cosmology intertwined II: The hubble constant tension_ , _Astropart. Phys._ 131 (2021) 102605 [2008.11284].
  * [13] N. Schöneberg, G. Franco Abellán, A. Pérez Sánchez, S. J. Witte, V. Poulin and J. Lesgourgues, _The $H_{0}$ Olympics: A fair ranking of proposed models_, 2107.10291.
  * [14] E. Mörtsell and S. Dhawan, _Does the Hubble constant tension call for new physics?_ , _JCAP_ 09 (2018) 025 [1801.07260].
  * [15] A. G. Riess, S. Casertano, W. Yuan, L. M. Macri and D. Scolnic, _Large Magellanic Cloud Cepheid Standards Provide a 1% Foundation for the Determination of the Hubble Constant and Stronger Evidence for Physics beyond $\Lambda$CDM_, _Astrophys. J._ 876 (2019) 85 [1903.07603].
  * [16] N. Blinov and G. Marques-Tavares, _Interacting radiation after Planck and its implications for the Hubble Tension_ , _JCAP_ 09 (2020) 029 [2003.08387].
  * [17] LSST Dark Matter Group collaboration, A. Drlica-Wagner et al., _Probing the Fundamental Nature of Dark Matter with the Large Synoptic Survey Telescope_ , 1902.01055.
  * [18] M. Sitwell, A. Mesinger, Y.-Z. Ma and K. Sigurdson, _The Imprint of Warm Dark Matter on the Cosmological 21-cm Signal_ , _Mon. Not. Roy. Astron. Soc._ 438 (2014) 2664 [1310.0029].
  * [19] PTOLEMY collaboration, M. G. Betti et al., _Neutrino physics with the PTOLEMY project: active neutrino properties and the light sterile case_ , _JCAP_ 07 (2019) 047 [1902.05508].
  * [20] D. McKeen, _Cosmic neutrino background search experiments as decaying dark matter detectors_ , _Phys. Rev. D_ 100 (2019) 015028 [1812.08178].
  * [21] A. J. Long, C. Lunardini and E. Sabancilar, _Detecting non-relativistic cosmic neutrinos by capture on tritium: phenomenology and physics potential_ , _JCAP_ 08 (2014) 038 [1405.7654].
  * [22] O. E. Bjaelde and S. Das, _Dark Matter Decaying into a Fermi Sea of Neutrinos_ , _Phys. Rev. D_ 82 (2010) 043504 [1002.1306].
  * [23] P. B. Greene and L. Kofman, _Preheating of fermions_ , _Phys. Lett. B_ 448 (1999) 6 [hep-ph/9807339].
  * [24] P. B. Greene and L. Kofman, _On the theory of fermionic preheating_ , _Phys. Rev. D_ 62 (2000) 123516 [hep-ph/0003018].
  * [25] L. Randall, J. Scholtz and J. Unwin, _Flooded Dark Matter and S Level Rise_ , _JHEP_ 03 (2016) 011 [1509.08477].
  * [26] T. Moroi and W. Yin, _Light Dark Matter from Inflaton Decay_ , _JHEP_ 03 (2021) 301 [2011.09475].
  * [27] T. Moroi and W. Yin, _Particle Production from Oscillating Scalar Field and Consistency of Boltzmann Equation_ , _JHEP_ 03 (2021) 296 [2011.12285].
  * [28] L. Randall, J. Scholtz and J. Unwin, _Cores in Dwarf Galaxies from Fermi Repulsion_ , _Mon. Not. Roy. Astron. Soc._ 467 (2017) 1515 [1611.04590].
  * [29] H. Davoudiasl, P. B. Denton and D. A. McGady, _Ultralight fermionic dark matter_ , _Phys. Rev. D_ 103 (2021) 055014 [2008.06505].
  * [30] M. Carena, N. M. Coyle, Y.-Y. Li, S. D. McDermott and Y. Tsai, _Cosmologically Degenerate Fermions_ , 2108.02785.
  * [31] N. Afshordi, M. Zaldarriaga and K. Kohri, _On the stability of dark energy with mass-varying neutrinos_ , _Phys. Rev. D_ 72 (2005) 065024 [astro-ph/0506663].
  * [32] A. Gogoi, R. K. Sharma, P. Chanda and S. Das, _Early Mass-varying Neutrino Dark Energy: Nugget Formation and Hubble Anomaly_ , _Astrophys. J._ 915 (2021) 132 [2005.11889].
  * [33] A. Y. Smirnov and X.-J. Xu, _Neutrino bound states and bound systems_ , 2201.00939.
  * [34] J. Preskill, M. B. Wise and F. Wilczek, _Cosmology of the Invisible Axion_ , _Phys. Lett. B_ 120 (1983) 127.
  * [35] L. F. Abbott and P. Sikivie, _A Cosmological Bound on the Invisible Axion_ , _Phys. Lett. B_ 120 (1983) 133.
  * [36] M. Dine and W. Fischler, _The Not So Harmless Axion_ , _Phys. Lett. B_ 120 (1983) 137.
  * [37] Planck collaboration, N. Aghanim et al., _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6 [1807.06209].
  * [38] S. Hannestad and G. Raffelt, _Constraining invisible neutrino decays with the cosmic microwave background_ , _Phys. Rev. D_ 72 (2005) 103514 [hep-ph/0509278].
  * [39] M. Archidiacono and S. Hannestad, _Updated constraints on non-standard neutrino interactions from Planck_ , _JCAP_ 07 (2014) 046 [1311.3873].
  * [40] A. Basboll, O. E. Bjaelde, S. Hannestad and G. G. Raffelt, _Are cosmological neutrinos free-streaming?_ , _Phys. Rev. D_ 79 (2009) 043512 [0806.1735].
  * [41] M. Escudero and M. Fairbairn, _Cosmological Constraints on Invisible Neutrino Decays Revisited_ , _Phys. Rev. D_ 100 (2019) 103531 [1907.05425].
  * [42] G. Barenboim, J. Z. Chen, S. Hannestad, I. M. Oldengott, T. Tram and Y. Y. Y. Wong, _Invisible neutrino decay in precision cosmology_ , _JCAP_ 03 (2021) 087 [2011.01502].
  * [43] A. Kogut et al., _The Primordial Inflation Explorer (PIXIE): A Nulling Polarimeter for Cosmic Microwave Background Observations_ , _JCAP_ 07 (2011) 025 [1105.2044].
  * [44] CMB-S4 collaboration, K. N. Abazajian et al., _CMB-S4 Science Book, First Edition_ , 1610.02743.
  * [45] D. Baumann, D. Green and M. Zaldarriaga, _Phases of New Physics in the BAO Spectrum_ , _JCAP_ 11 (2017) 007 [1703.00894].
  * [46] T. Bringmann, F. Kahlhoefer, K. Schmidt-Hoberg and P. Walia, _Converting nonrelativistic dark matter to radiation_ , _Phys. Rev. D_ 98 (2018) 023543 [1803.03644].
  * [47] K. Enqvist, S. Nadathur, T. Sekiguchi and T. Takahashi, _Constraints on decaying dark matter from weak lensing and cluster counts_ , _JCAP_ 04 (2020) 015 [1906.09112].
  * [48] S. Aoyama, T. Sekiguchi, K. Ichiki and N. Sugiyama, _Evolution of perturbations and cosmological constraints in decaying dark matter models with arbitrary decay mass products_ , _JCAP_ 07 (2014) 021 [1402.2972].
  * [49] B. Audren, J. Lesgourgues, G. Mangano, P. D. Serpico and T. Tram, _Strongest model-independent bound on the lifetime of Dark Matter_ , _JCAP_ 12 (2014) 028 [1407.2418].
  * [50] K. Enqvist, S. Nadathur, T. Sekiguchi and T. Takahashi, _Decaying dark matter and the tension in $\sigma_{8}$_, _JCAP_ 09 (2015) 067 [1505.05511].
  * [51] T. Simon, G. Franco Abellán, P. Du, V. Poulin and Y. Tsai, _Constraining decaying dark matter with BOSS data and the effective field theory of large-scale structures_ , _Phys. Rev. D_ 106 (2022) 023516 [2203.07440].
  * [52] S. Alvi, T. Brinckmann, M. Gerbino, M. Lattanzi and L. Pagano, _Do you smell something decaying? Updated linear constraints on decaying dark matter scenarios_ , 2205.05636.
  * [53] M. G. Walker, M. Mateo, E. W. Olszewski, J. Penarrubia, N. W. Evans and G. Gilmore, _A universal mass profile for dwarf spheroidal galaxies?_ , _arXiv preprint arXiv:0906.0341_ (2009) .
  * [54] L. Kofman, A. D. Linde and A. A. Starobinsky, _Reheating after inflation_ , _Phys. Rev. Lett._ 73 (1994) 3195 [hep-th/9405187].
  * [55] L. Kofman, A. D. Linde and A. A. Starobinsky, _Towards the theory of reheating after inflation_ , _Phys. Rev. D_ 56 (1997) 3258 [hep-ph/9704452].
  * [56] J. Baacke, K. Heitmann and C. Patzold, _Nonequilibrium dynamics of fermions in a spatially homogeneous scalar background field_ , _Phys. Rev. D_ 58 (1998) 125013 [hep-ph/9806205].
  * [57] J. Berges, D. Gelfand and J. Pruschke, _Quantum theory of fermion production after inflation_ , _Phys. Rev. Lett._ 107 (2011) 061301 [1012.4632].
  * [58] J. Berges, D. Gelfand and D. Sexty, _Amplified Fermion Production from Overpopulated Bose Fields_ , _Phys. Rev. D_ 89 (2014) 025001 [1308.2180].
  * [59] R. Daido, F. Takahashi and W. Yin, _The ALP miracle: unified inflaton and dark matter_ , _JCAP_ 05 (2017) 044 [1702.03284].
  * [60] F. Takahashi and W. Yin, _QCD axion on hilltop by a phase shift of $\pi$_, _JHEP_ 10 (2019) 120 [1908.06071].
  * [61] D. J. E. Marsh and W. Yin, _Opening the 1 Hz axion window_ , _JHEP_ 01 (2021) 169 [1912.08188].
  * [62] S. Nakagawa, F. Takahashi and W. Yin, _Stochastic Axion Dark Matter in Axion Landscape_ , _JCAP_ 05 (2020) 004 [2002.12195].
  * [63] B. Batell and A. Ghalsasi, _Thermal misalignment of scalar dark matter_ , _Phys. Rev. D_ 107 (2023) L091701 [2109.04476].
  * [64] B. Batell, A. Ghalsasi and M. Rai, _Dynamics of dark matter misalignment through the Higgs portal_ , _JHEP_ 01 (2024) 038 [2211.09132].
  * [65] C. Murgui and R. Plestid, _Coleman-Weinberg dynamics of ultralight scalar dark matter and GeV-scale right-handed neutrinos_ , 2306.13799.