# Freeze-in production of decaying dark matter in five steps

Saniya Heeba  Felix Kahlhoefer  and Patrick Stöcker

###### Abstract

We study the cosmological evolution and phenomenological properties of scalar
bosons in the keV to MeV range that have a tiny mixing with the Standard Model
Higgs boson. The mixing determines both the abundance of light scalars
produced via the freeze-in mechanism and their lifetime. Intriguingly, the
parameters required for such scalars to account for all of the dark matter in
the present Universe generically predict lifetimes comparable to the
sensitivity of present and future indirect detection experiments. In order to
accurately determine the relic abundance of light scalars, we calculate
freeze-in yields including effects from finite temperatures and quantum
statistics and develop a new approach for solving the Boltzmann equation for
number-changing processes in the dark sector. We find that light scalars can
potentially explain the anomalous x-ray emission at 3.5 keV, while evading
constraints from structure formation and predicting potentially observable
self-interaction cross sections.

##  1 Introduction

All known particles of the Standard Model (SM) have sufficiently strong
interactions to enter into thermal equilibrium with each other at high
temperatures. It is hence tempting to assume that also the dark matter (DM)
particle was at some point in thermal equilibrium with the SM and then
obtained its relic abundance via the freeze-out mechanism. While this scenario
is both well-motivated and predictive, there are many interesting alternative
production mechanisms for DM. In particular, it is conceivable that reheating
created only a negligible DM abundance and that the dark sector is only slowly
populated subsequently through tiny interactions with SM particles. This so-
called freeze-in mechanism [1, 2] has received significant interest recently,
illustrated best by the fact that micrOmegas now provides a numerical tool for
the automated calculation of freeze-in yields [3].

While freeze-in is an attractive alternative to the standard paradigm of
thermal freeze-out, an experimental confirmation of the idea is rather
challenging. Colliders such as the LHC are typically unable to probe non-
thermalised hidden sectors [4, 5] and direct detection experiments are only
sensitive if the mediator of the interactions between the DM particle and the
SM is light, such that scattering is strongly enhanced in the non-relativistic
limit [6, 7, 8]. Indirect detection, however, may have a unique opportunity to
probe freeze-in [9], for example if the DM particle is unstable, such that it
can decay via the same interaction responsible for its production [10]. In
this case, the smallness of the coupling implied by the freeze-in mechanism
explains at the same time the long lifetime of DM in the absence of a
stabilising symmetry.

Such a set-up is well-known in the context of keV sterile neutrinos, for which
production and decay both proceed through a tiny mixing with SM neutrinos [11]
(see also Refs. [12, 13, 14] for alternative ways to produce sterile neutrinos
via freeze-in). The case of scalar and vector particles that mix with the SM
Higgs boson or the SM photon, respectively, was first explored in Ref.
[15].111The case of _stable_ scalar singlets produced via the Higgs portal,
was considered previously in Refs. [16, 17, 18]. This work however focused on
the case of a sub-dominant DM component with a lifetime short compared to the
age of the Universe. Moreover, the freeze-in production of light scalars was
calculated in a very approximate way, considering only a limited range of
temperatures and a few production processes.

In the present work we point out that the production in the Early Universe of
light scalars with Higgs mixing is in fact surprisingly complex. First of all,
the freeze-in production proceeds in three stages, corresponding to
temperatures before, during and after electroweak symmetry breaking (EWSB).
Although the dominant contribution to the DM abundance typically arises after
EWSB, it is essential to correctly account for the effect of the electroweak
phase transition (EWPT) in order to avoid unphysical contributions from high
temperatures. Furthermore, since the light scalars are not protected by a
stabilising symmetry, $2\to 3$ and $3\to 2$ processes may play an important
role. This means that the co-moving DM density is not necessarily constant
after the end of freeze-in and additional considerations are needed to
calculate the subsequent evolution of the dark sector.

A number of previous works have studied the freeze-out of number-changing
processes in a dark sector that is initially in kinetic and chemical
equilibrium [19, 20, 21, 22]. In the present work, we extend these studies by
considering the evolution of a dark sector where chemical equilibrium is not
guaranteed.222Refs. [18, 23] provide a similar discussion in the context of
$2\to 4$ processes, while Ref. [24] considers $2\to 3$ processes for the case
of vector DM. For this purpose, we consider $2\to 3$ processes for
relativistic initial states in order to address the question whether or not
number-changing processes are efficient enough to thermalise the dark sector.
Using a combination of analytical approximations and numerical algorithms we
can then calculate the present day abundance of light scalars for arbitrary
model parameters.

Figure 1: Sketch of the evolution of the co-moving number density of light
scalars with decreasing temperature of the SM thermal bath. Note that the
relative size of the different regimes and the magnitude of the different
effects is not to scale.

The different stages of production of light scalars are illustrated in figure
1, which shows schematically the evolution of the co-moving number density $Y$
as a function of the inverse temperature $x$ (precise definitions of these
quantities will be provided below). The first three steps correspond to the
freeze-in production of light scalars before, during and after EWSB. The
fourth and fifth step correspond to thermalisation of the dark sector and dark
sector freeze-out, respectively. While the first three steps always increase
the DM abundance and the final step always leads to a decrease, the fourth
step can either enhance or deplete the co-moving number density, depending on
how the freeze-in yield compares to the equilibrium distribution.

We identify large regions of parameter space where lights scalars can be all
of DM and have a lifetime large compared to the age of the Universe. The
combined constraints from indirect detection experiments and DM self-
interactions, however, force such light scalars to have masses in the keV
range and lifetimes that are only slightly below current experimental
sensitivity, such that the scenario may be testable with future x-ray missions
such as the Hitomi [25] re-flight. Moreover, decaying light scalars from
freeze-in may provide a viable explanation of the claimed observation of an
x-ray line at 3.5 keV in various astrophysical systems [26, 27] while facing
much weaker constraints from structure formation than sterile neutrinos. For
different ways of connecting the 3.5 keV line to the freeze-in mechanism, we
refer to Refs. [28, 29].

This paper is structured as follows. We present the model that we consider in
section 2, derive the mixing with the SM Higgs boson after electroweak
symmetry breaking and discuss the importance of finite-temperature effects.
Section 3 deals with the freeze-in production of light scalars, accounting for
the different production mechanisms before, during and after electroweak
symmetry breaking. The subsequent evolution, in particular the role of $2\to
3$ processes, is discussed in section 4, where we also provides the final
results from the relic density calculation. The phenomenological implications
and existing experimental constraints are discussed in Sec. 5. Appendix A
provides additional details on how to solve the Boltzmann equation for $2\to
3$ processes.

##  2 Light scalars and Higgs mixing

We consider the most general Lagrangian for a real scalar singlet $s_{0}$ with
a $\mathbb{Z}_{2}$ symmetry (see also Ref. [30]):

$\displaystyle\mathcal{L}_{\text{scalar}}$
$\displaystyle=\frac{1}{2}\partial^{\mu}s_{0}\partial_{\mu}s_{0}+\frac{1}{2}\mu_{s}^{2}\,s_{0}^{2}-\frac{1}{4}\lambda_{s}\,s_{0}^{4}-\frac{1}{2}\lambda_{hs}\,|H|^{2}s_{0}^{2}\;.$
(2.1)

The negative mass term for the scalar leads to a spontaneous breaking of the
$\mathbb{Z}_{2}$, such that the scalar acquires a vacuum expectation value
(vev): $s_{0}=s+v_{s}$. We assume for the moment that this happens at
temperatures large compared to the temperature of the EWPT. The Lagrangian
then becomes

$\displaystyle\mathcal{L}_{\text{scalar}}$
$\displaystyle=\frac{1}{2}\partial^{\mu}s\partial_{\mu}s-\frac{1}{2}m_{s}^{2}\,s^{2}-\frac{1}{4}\lambda_{s}\,s^{4}-\lambda_{s}\,v_{s}\,s^{3}-\frac{1}{2}\lambda_{hs}\,|H|^{2}(s^{2}+2\,s\,v_{s})\;.$
(2.2)

Here we have assumed that $\lambda_{hs}\,v_{s}^{2}$ is very small compared to
the bare Higgs mass term and can therefore be neglected, which is a good
approximation for the values of $\lambda_{hs}$ and $v_{s}$ that we will
consider later.

After EWSB we replace $H=(h+v)/\sqrt{2}$. This introduces an additional
contribution to the mass term of the scalar singlet:
$-\mu_{s}^{2}\to-\mu_{s}^{2}+\lambda_{hs}v^{2}/2$. This contribution can
potentially restore the $\mathbb{Z}_{2}$ symmetry, leading to the so-called
vev flip-flop [31, 32]. Here, instead, we will consider the case where the
$\mathbb{Z}_{2}$ remains broken at low temperatures. In this case, the
Lagrangian can be written as

$\displaystyle\mathcal{L}_{\text{scalar}}$
$\displaystyle=\frac{1}{2}\partial^{\mu}s\partial_{\mu}s-\frac{1}{2}m_{s}^{2}\,s^{2}-\frac{1}{4}\lambda_{s}\,s^{4}-\lambda_{s}\,v_{s}\,s^{3}-\frac{1}{4}\lambda_{hs}\,(h^{2}+2\,h\,v)(s^{2}+2\,s\,v_{s})\;.$
(2.3)

We take the free parameters to be $\lambda_{s}$ and $\lambda_{hs}$ as well as
the zero-temperature value of $m_{s}$, such that the vev of the scalar field
is given by

$\displaystyle v_{s}$ $\displaystyle=\frac{m_{s}}{\sqrt{2\,\lambda_{s}}}\;.$
(2.4)

Eq. (2.3) still contains a non-diagonal mass term of the form
$\mathcal{L}_{\text{mixing}}=\lambda_{hs}\,v_{s}v\,\,hs$, leading to the
mixing of the two scalar fields. The mass matrix is given by

$M=\begin{pmatrix}2v^{2}\lambda&\lambda_{hs}v_{s}v\\\
\lambda_{hs}v_{s}v&2v_{s}^{2}\lambda_{s}\end{pmatrix}$ (2.5)

with $\lambda$ denoting the quartic self-coupling of the SM Higgs field. To
find the mass eigenstates we rotate the interaction eigenstates by the mixing
angle $\theta$:

$\begin{pmatrix}h_{\text{SM}}\\\
h_{s}\end{pmatrix}=\begin{pmatrix}\cos\theta&\sin\theta\\\
-\sin\theta&\cos\theta\end{pmatrix}\begin{pmatrix}h\\\ s\end{pmatrix}\;.$
(2.6)

The mass matrix is diagonalised for

$\tan
2\theta\equiv\frac{\lambda_{hs}\,v\,v_{s}}{\lambda\,v^{2}-\lambda_{s}\,v_{s}^{2}}\;,$
(2.7)

which gives the mass eigenvalues

$\displaystyle m^{2}_{h,s}=\lambda v^{2}+\lambda_{s}v_{s}^{2}\pm\sqrt{(\lambda
v^{2}-\lambda_{s}v_{s}^{2})^{2}+\lambda_{hs}^{2}v_{s}^{2}v^{2}}\;.$ (2.8)

We will focus on the case $\theta\ll 1$, such that $h_{\text{SM}}$ is SM-like,
while $h_{s}$ only has tiny couplings to SM states. In this case we can
furthermore neglect the shift in the mass eigenstates.333In particular, the
relative change of $m_{s}$ will be of order $\lambda_{hs}^{2}/\lambda_{s}$,
which is completely negligible. Thus, to simplify notation, we will continue
denoting the mostly SM-like scalar by $h$ and the mostly SM-singlet scalar by
$s$. We also note that in this case eq. (2.7) simplifies to

$\theta=\frac{\lambda_{hs}\,v\,v_{s}}{m^{2}_{h}-m^{2}_{s}}\;.$ (2.9)

While the discussion above has been fully general, we will from now on focus
on the case where the extra scalar is much lighter than the SM Higgs boson,
$m_{s}\ll m_{h}$. The main motivation for this choice is that we will be
interested in the case of long-lived light scalars, which – as we will see
below – requires masses in the keV to MeV range.

We make the crucial observation that the mixing angle $\theta$ depends on
temperature via the temperature dependence of the two vevs. In particular,
$\theta=0$ for $v=0$, i.e. before EWSB. This means that the phenomenology of
our model will be decisively different at high temperatures, when the two
scalar fields cannot mix, and at low temperatures, when mixing becomes
possible. We will treat these two regimes separately below and also discuss in
detail what happens for temperatures close to the EWPT.

The temperature dependence of the electroweak vev can be determined by
considering the effective potential $V$. Assuming the additional scalar field
to be so weakly coupled that its contribution is completely negligible, we
recover the well-known result for the SM [33]:

$V(\phi,T)=D(T^{2}-T_{o}^{2})\phi^{2}-ET\phi^{3}+\frac{\lambda(T)}{4}\phi^{4}\;,$
(2.10)

where $\phi$ is the constant background field and $D,E,T_{o}$ and $\lambda(T)$
can be calculated in terms of SM parameters. Since we are not interested in
the strength of the EWPT, we will ignore the term proportional to $\phi^{3}$,
making use of the fact that $E$ is known to be small in the SM. The
temperature of the EWPT is then simply given by $T_{\mathrm{c}}=T_{o}\approx
164\,\mathrm{GeV}$. For $T>T_{\mathrm{c}}$ the effective potential is
minimized by $\phi=0$, such that the electroweak symmetry is unbroken. While
fermions and gauge bosons are massless, the mass of the (complex) Higgs boson
is given by

$m_{H}(T)^{2}=D(T^{2}-T_{o}^{2})\;.$ (2.11)

For temperatures smaller than but close to the critical temperature, the Higgs
vev and mass are simply given by

$\displaystyle v(T)^{2}$
$\displaystyle=\frac{2D(T_{o}^{2}-T^{2})}{\lambda(T)}\;,$ (2.12)
$\displaystyle m_{h}(T)^{2}$
$\displaystyle=2\lambda(T)v(T)^{2}=4D(T_{o}^{2}-T^{2})\;.$ (2.13)

We will make use of these expressions in the following section.

##  3 Freeze-in production of light scalars

We will be interested in the case where the initial abundance of light scalars
(e.g. at the end of reheating) is completely negligible, such that the
dominant contribution to their relic abundance stems from the “leakage” of
energy from the visible into the dark sector. There are three relevant epochs
to this freeze-in production, corresponding to the temperature-dependent
properties of the SM Higgs boson. As long as the electroweak symmetry is
unbroken, the SM Higgs boson and the light scalar cannot mix and hence only a
small number of production channels are available. During444Strictly speaking,
oscillations happen immendiately _after_ the EWPT, but the relevant
temperatures are very close to $T_{\mathrm{c}}$. the EWPT, the mixing between
the two scalar bosons can be large for a very short period of time,
potentially leading to rapid production of light scalars via oscillations. At
lower temperatures the mixing angle becomes small but non-zero, opening up a
wide range of possible production modes.

###  3.1 Step 1: Production before EWSB

Figure 2:  Processes relevant for the freeze-in production of light scalars
before electroweak symmetry breaking, with $V$ denoting a gauge boson of
$SU(2)$ or $U(1)$. The four-point interaction (a) gives a larger contribution
than diagrams with a $t$-channel Higgs boson (b), (c).

Before EWSB, the SM Higgs field has vanishing vev and therefore cannot mix
with the light scalar. This means that there are only two types of processes
that can contribute to the freeze-in production of light scalars: the four-
point interaction $H^{\dagger}Hs^{2}$ and diagrams with a complex Higgs boson
in the $t$-channel, for example $H^{\dagger}H\to sV$ or $HV\to Hs$, with $V$
being a (massless) gauge boson of either $SU(2)$ or $U(1)$ (see figure 2).
Conversely, diagrams like $t\bar{t}\to gs$ with a fermion in the $t$-channel
are absent, because the light scalar does not at this stage couple to top
quarks. We therefore find that all relevant cross sections vanish in the limit
$\sqrt{S}\to\infty$, such that freeze-in production proceeds dominantly at low
temperatures. This means in particular that in the model we consider the
freeze-in yield is independent of the reheating temperature $T_{\mathrm{r}}$,
as long as $T_{\mathrm{r}}\gg T_{\mathrm{c}}$.

Since the mixing angle $\theta$ vanishes before EWSB and since the mass of the
light scalar can be neglected at high temperatures, the only parameters that
determine the yield from freeze-in before EWSB are the mixed quartic coupling
$\lambda_{hs}$ and the vev $v_{s}$ of the light scalar. We find that processes
involving a $t$-channel Higgs boson are completely negligible, so that the
freeze-in yield is dominated by the process $H^{\dagger}H\to ss$. To leading
order in $\lambda_{hs}$ the corresponding cross section as a function of the
centre-of-mass energy $\sqrt{S}$ is given by

$\displaystyle\sigma_{H^{\dagger}H\to ss}$
$\displaystyle=\frac{\lambda_{hs}^{2}}{32\pi
S}\frac{(S+2m_{s}^{2})^{2}}{(S-m_{s}^{2})^{2}}\sqrt{\frac{S-4m_{s}^{2}}{S-4m_{H}^{2}}}\approx\frac{\lambda_{hs}^{2}}{32\pi\sqrt{S^{2}-4\,S\,m_{H}^{2}}}\;,$
(3.1)

which is in fact independent of $v_{s}$.

The number density of light scalars $n_{s}$ can then be obtained by solving
the Boltzmann equation

$\displaystyle\dot{n}_{s}+3\,H\,n_{s}=C_{H}\,g_{H}\int
C[f_{H}]\frac{\mathrm{d}^{3}p}{(2\pi)^{3}}\;.$ (3.2)

Here $g_{H}=2$ denotes the degrees of freedom of the Higgs field, $C_{H}=2$
reflect the fact that there are two non-identical particles in the initial
state and $H$ is the Hubble expansion rate during radiation domination:

$\displaystyle
H=\sqrt{g_{\mathrm{SM}}^{\ast}}\frac{1.66\,T_{\mathrm{SM}}^{2}}{M_{\mathrm{Pl}}}$
(3.3)

with $M_{\mathrm{Pl}}$ being the Planck mass and $T_{\mathrm{SM}}$ and
$g^{\ast}_{\mathrm{SM}}$ denoting the temperature and number of relativistic
degrees of freedom of the thermal bath of SM particles.

The quantity $C[f_{H}]$ represents the collision term responsible for creating
light scalars from the thermal bath. To obtain a simple expression, we
introduce the dimensionless quantities, $Y=n/s_{\mathrm{SM}}$ and
$x_{\mathrm{SM}}=m_{s}/T_{\mathrm{SM}}$, where $s_{\mathrm{SM}}$ is the
entropy density of the visible sector:

$\displaystyle
s_{\mathrm{SM}}=g_{\mathrm{SM}}^{\ast}\frac{2\pi^{2}\,T_{\mathrm{SM}}^{3}}{45}\;.$
(3.4)

We can then use the methods outlined in Ref. [1] to simplify eq. (3.2) and
obtain

$\displaystyle\frac{\mathrm{d}Y_{s}}{\mathrm{d}x_{\mathrm{SM}}}$
$\displaystyle=\frac{2\,s_{\mathrm{SM}}}{H\,x_{\mathrm{SM}}}\langle\sigma
v\rangle\,Y_{H}^{2}\;.$ (3.5)

Figure 3: Differential yield of light scalars before EWSB as a function of
temperature for the dominant production channel $H^{\dagger}H\to ss$. The blue
solid line corresponds to the result obtained when accounting for quantum
statistics (i.e. using Bose-Einstein and Fermi-Dirac distributions), while the
orange dashed line represents the results obtained when assuming a Maxwell-
Boltzmann distribution. For this plot we have set $\lambda_{hs}=2\times
10^{-10}$.

Calculating the thermally averaged cross-section $\langle\sigma v\rangle$ in a
similar fashion as in Ref. [34], we find the differential production rate of
DM to be given by

$\displaystyle\frac{\mathrm{d}Y_{s}}{\mathrm{d}x_{\mathrm{SM}}}$
$\displaystyle=\frac{90\,M_{\mathrm{Pl}}}{(2\pi)^{6}\,1.66\,(g_{\mathrm{SM}}^{\ast})^{3/2}}\frac{x^{3}}{m_{s}^{4}}\int_{S_{min}}^{\infty}\mathrm{d}S\,(S-4\,m_{H}^{2})\,\sigma\,\sqrt{S}\,K_{1}\left(\sqrt{S}\
x/m_{s}\right)\;,$ (3.6)

where $K_{1}(x)$ is the Modified Bessel function. We show
$\mathrm{d}Y_{s}/\mathrm{d}\log T_{\mathrm{SM}}$ as a function of
$T_{\mathrm{SM}}$ in figure 3. As expected for freeze-in via renormalisable
interactions, we find that the production rate peaks at low temperatures, i.e.
right before the EWPT. Eq. (3.6) can be immediately integrated from $x=0$ to
$x=m_{s}/T_{\mathrm{c}}$ to obtain the total amount of DM production before
EWSB.

We find that for $m_{s}\ll T_{\mathrm{c}}$ the number of DM particles produced
before EWSB is to good approximation independent of the scalar mass and can be
written as

$\displaystyle\Delta Y_{s,1}\approx(3.58\times 10^{9})\ \lambda_{hs}^{2}\;.$
(3.7)

We emphasize that this result consistently takes into account the fact that
the mass of the Higgs bosons in the initial state depends on temperature.
However, since $m_{H}(T)<T$ before EWSB, it is not a good approximation to
treat the Higgs bosons in the initial state as non-relativistic, which we have
done implicitly by using Maxwell-Boltzmann statistics to calculate the
thermally-averaged annihilation cross section.

In order to properly account for the statistics of indistinguishable
relativistic particles, we use micrOmegas [3], which employs Bose-Einstein
statistics to calculate the relic abundance more accurately. It is however not
easily possible to implement a temperature-dependent Higgs mass in micrOmegas.
We therefore determine the value of the Higgs boson mass that gives the same
result as the temperature-dependent Higgs mass for the case of Maxwell-
Boltzmann statistics, and then fix the Higgs boson mass to this value in
micrOmegas also for the case of Bose-Einstein statistics. The suitable value
is found to be $m_{H,\text{fixed}}\approx m_{H}(T=287\,\mathrm{GeV})\approx
136\,\mathrm{GeV}$, which is consistent with the observation that the dominant
contribution to freeze-in production arises shortly before the EWPT.

Following this procedure, we find that the abundance with the more accurate
treatment of statistics is given by

$\displaystyle\Delta Y_{s,1}^{\text{BE}}\approx(6.77\times 10^{9})\
\lambda_{hs}^{2}\;,$ (3.8)

which is larger than the previous result by about a factor of 2 (see figure
3). Such a difference is not unexpected for relativistic particles in the
initial state and is consistent with similar observations made in Ref. [3].

###  3.2 Step 2: Production during EWSB

Electroweak symmetry breaking is triggered by the $\mu^{2}$-term in the
Lagrangian changing sign. At the point when this happens, $T\approx
T_{\mathrm{c}}$, the mass of the SM Higgs boson vanishes and then increases
again with decreasing temperature. This implies in particular that right after
the EWPT the SM Higgs boson will have both a very small (but non-zero) vev
$v(T)$ and a very small mass $m_{h}(T)$. The temperature-dependent mixing
angle is then approximately given by

$\theta(T)=\frac{\lambda_{hs}\,v(T)\,v_{s}}{m_{h}^{2}(T)-m_{s}^{2}}\;,$ (3.9)

where we assume for simplicity that the temperature dependence of $m_{s}$ and
$v_{s}$ is negligible for $T\approx T_{\mathrm{c}}$.

For a short period of time right after the EWPT, the mass of the SM Higgs
boson may be comparable to the mass of the light scalar, $m_{h}(T)\approx
m_{s}$. It then becomes important to include also the imaginary part of the
Higgs boson mass, which depends on the interactions of the Higgs with other
particles in the thermal bath. While this term ensures that $\theta(T)\ll 1$
(assuming sufficiently small $\lambda_{hs}$), the mixing between the two
scalars will be substantially enhanced. At this point, SM Higgs bosons can
directly be converted into light scalars via oscillations.

This production mechanism is analogous to the one of hidden photons via mixing
with the SM photon, which is enhanced at finite temperatures and densities due
to plasma effects. It was shown in Ref. [35] that the resulting abundance of
hidden photons depends essentially only on the mixing parameter and on how
quickly the mass of the visible photon varies with temperature. The
corresponding result for our case can be written as

$\Delta
Y_{s,2}=\frac{\pi\,\zeta(2)\,Y_{h}(T_{\mathrm{osc}})}{\zeta(3)}\frac{\lambda_{hs}^{2}\,v(T_{\mathrm{osc}})^{2}\,v_{s}^{2}}{H(T_{\mathrm{osc}})\,T_{\mathrm{osc}}^{2}}\left|\frac{\mathrm{d}m_{h}^{2}}{\mathrm{d}T}\right|^{-1}_{T=T_{\mathrm{osc}}}\;,$
(3.10)

where $T_{\mathrm{osc}}$ is the temperature for which $m_{h}(T)=m_{s}$. Since
$T_{\mathrm{osc}}\approx T_{\mathrm{c}}$, we define
$\epsilon=(T_{\mathrm{c}}-T_{\mathrm{osc}})/T_{\mathrm{c}}\ll 1$. We then find
$m_{h}^{2}\approx 8\,D\,T_{\mathrm{c}}^{2}\,\epsilon$ and hence

$\epsilon=\frac{m_{s}^{2}}{8\,D\,T_{\mathrm{c}}^{2}}\;.$ (3.11)

For $\epsilon\ll 1$ we find furthermore

$v(\epsilon)^{2}=4\
T_{\mathrm{c}}^{2}\frac{D\,\epsilon}{\lambda(T_{\mathrm{c}})}$ (3.12)

and hence

$\frac{\lambda_{hs}^{2}\,v(T_{\mathrm{osc}})^{2}\,v_{s}^{2}}{H(T_{\mathrm{osc}})\,T_{\mathrm{osc}}^{2}}\left|\frac{\mathrm{d}m_{h}^{2}}{\mathrm{d}T}\right|^{-1}_{T=T_{\mathrm{osc}}}\approx\frac{\lambda_{hs}^{2}\,m_{s}^{2}\,v_{s}^{2}}{16\,D\,\lambda(T_{\mathrm{c}})\,H(T_{\mathrm{c}})\,T_{\mathrm{c}}^{3}}\;.$
(3.13)

We can furthermore make use of the fact that the Higgs boson is highly
relativistic at $T=T_{\mathrm{osc}}$ and hence

$\frac{\pi\,\zeta(2)\,Y_{h}(T_{\mathrm{osc}})}{\zeta(3)}\approx\frac{\zeta(2)\,T_{\mathrm{c}}^{3}}{\pi\,s_{\text{SM}}(T_{\mathrm{c}})}\;.$
(3.14)

Bringing everything together, we obtain for the DM yield from oscillations

$\Delta Y_{s,2}\approx(1.93\times
10^{5}\,\mathrm{GeV^{-4}})\,\lambda_{hs}^{2}\,m_{s}^{2}\,v_{s}^{2}\;.$ (3.15)

For the parameter region that we will be interested in ($m_{s}\ll
v_{s}\sim\mathrm{GeV}$) this yield is small compared to the one from before
EWSB, such that the final DM abundance does not change significantly during
EWSB.

###  3.3 Step 3: Production after EWSB

Figure 4: Processes relevant for freeze-in production after the electroweak
phase transition. For $v_{s}<100\,\mathrm{GeV}$, production is dominated by
Higgs decays (a). Conversely, for $v_{s}>100\,\mathrm{GeV}$, a number of
diagrams like the ones shown in b) and c) give relevant contributions.

As the temperature decreases further, oscillations between SM Higgs bosons and
the light scalars become suppressed, but a number of new production channels
open up (see figure 4 for a few examples). First of all, light scalars can now
be produced in decays of SM Higgs bosons: $h\to ss$. Moreover, processes like
$t\bar{t}\to sg$ now become active and contribute to the abundance of light
scalars. We note that, in contrast to the $t$-channel processes considered in
Sec. 3.1, these processes grow with energy like $\log(S/m_{t}^{2})$, due to
the divergent contribution in the collinear limit $\cos\theta\to 0$. In the
absence of a phase transition at high temperatures, these processes would
therefore dominate the production of light scalars. However, since these
processes become active only after EWSB, their contribution is less important.

Interestingly, the two types of processes depend in different ways on the
fundamental parameters. The production from Higgs boson decays is proportional
to $\lambda_{hs}\,v\propto\sin\theta/v_{s}$, whereas the production from SM
fermions is simply proportional to $\sin\theta$. Hence, by independently
varying $\lambda_{hs}$ and $v_{s}$, it is possible to divide the parameter
space into two regions on the basis of which channel(s) gives the dominant
contribution. We find that for $v_{s}\leq 100\ \mathrm{GeV}$, or equivalently
$\lambda_{s}\geq 5\times 10^{-5}(m_{s}/1\,\mathrm{GeV})^{2}$, the dominant
contribution ($\geq 50\%$) comes from Higgs decays:

$\displaystyle\Gamma_{h\rightarrow ss}=\frac{\sin^{2}\theta}{v_{s}^{2}}\
\frac{\sqrt{m_{h}^{2}-4m_{s}^{2}}\ (m_{h}^{2}+2m_{s}^{2})^{2}}{32\pi\
m_{h}^{2}}\;.$ (3.16)

The corresponding yield can be calculated using a simplified Boltzmann
equation [1]:

$\displaystyle\frac{\mathrm{d}Y^{\mathrm{decay}}_{s,3}}{\mathrm{d}x}=\frac{45\,M_{\mathrm{Pl}}}{2\pi^{4}\,1.66\,(g_{\mathrm{SM}}^{\ast})^{3/2}}\frac{m_{h}^{2}}{m_{s}^{4}}\,\Gamma_{h\rightarrow
ss}\ x^{3}K_{1}\left(m_{h}\,x/m_{s}\right)\;.$ (3.17)

We show $\mathrm{d}Y_{s}/\mathrm{d}\log T_{\mathrm{SM}}$ as a function of
$T_{\mathrm{SM}}$ in figure 5. Integrating this expression from
$x=m_{s}/T_{\mathrm{c}}$ to $x\to\infty$, we obtain

$\displaystyle\Delta Y^{\mathrm{decay}}_{s,3}\sim\left(9\times 10^{15}\
\mathrm{GeV}^{2}\right)\ \frac{\sin^{2}\theta}{v_{s}^{2}}\ \sim 2.2\times
10^{12}\ \lambda_{hs}^{2}\;.$ (3.18)

Figure 5:  Differential yield of light scalars after EWSB as a function of
temperature for a number of different production channels. Blue solid lines
correspond to the results obtained when accounting for quantum statistics
(i.e. using Bose-Einstein and Fermi-Dirac distributions), while orange dashed
lines represent the results obtained when assuming a Maxwell-Boltzmann
distribution. For this plot we have set $\lambda_{hs}=2\times 10^{-10}$ and
$v_{s}=25\,\mathrm{GeV}$. For these parameters, the dominant contribution
arises from the decay $h\to ss$, which peaks around $T_{\text{SM}}\sim
40\,\mathrm{GeV}$.

Conversely, in the opposite regime ($v_{s}\geq 100\ \mathrm{GeV}$), production
is divided among a number of channels, with each channel contributing less
than $10\%$ to the total yield. We perform the relevant calculations using
micrOmegas which identifies more than 100 possible channels. For reference,
one such channel is $W^{+}b\rightarrow ts$ with a cross section proportional
to $\sin^{2}\theta$. In the limit $m_{s}\ll m_{h}$, the total yield from these
processes is calculated to be555Our result differs from the one obtained in
Ref. [15], where only two production channels ($t\bar{t}\to sg$ and $tg\to
ts$) are included. Moreover, Ref. [15] appears not to correctly take into
account the temperature dependence of the production rate, which peaks at
$T\gg m_{s}$ (see figure 5).

$\displaystyle\Delta Y^{2\rightarrow 2}_{s,3}\sim 7.2\times
10^{11}\,\sin^{2}\theta\ \sim\left(1.8\times 10^{8}\
\mathrm{GeV}^{-2}\right)\,\lambda_{hs}^{2}\,v_{s}^{2}\;.$ (3.19)

We also find that including quantum statistics affects the two regimes
differently. In case of decay-dominated production, including quantum
statistics increases the yield by approximately $20\%$. In the second regime,
the resulting change is less than $1\%$. The reason is that a large number of
different channels contribute, which behave in contrary ways on inclusion of
quantum statistics. For example, the yield from channel $W^{+}Z\rightarrow
W^{+}s$ increases whereas that from $W^{+}b\rightarrow ts$ decreases (see
figure 5), so the net effect is washed out.

We can now calculate the total yield from freeze-in by combining the results
obtained above:

$Y_{s,\text{fi}}=\Delta Y_{s,1}+\Delta Y_{s,2}+\Delta
Y_{s,3}^{\text{decay}}+\Delta Y_{s,3}^{2\to 2}$ (3.20)

From eqs. (3.7), (3.15), (3.18) and (3.19), it is easy to see that DM
production is dominated by processes after EWSB, i.e. $Y_{s,\text{fi}}\approx
Y_{s,3}$. For $v_{s}\gg 100\ \mathrm{GeV}$, the contributions from step 1 and
2 could be comparable but remain smaller than that of step 3.

Under the assumption that the initial DM yield after reheating is negligible
compared to yield from freeze-in, we obtain $Y_{s}=Y_{s,\text{fi}}$, which can
be directly converted into the present-day relic abundance $\Omega_{s}$:

$\Omega_{s}=\frac{m_{s}\,s_{\mathrm{SM},0}}{\rho_{\mathrm{c},0}}Y_{s}\;,$
(3.21)

where $s_{\mathrm{SM},0}$ and $\rho_{\mathrm{c},0}$ denote the present-day
entropy density and present-day critical density, respectively.666Here we
assume implicitly that the scalars have a lifetime that is large compared to
the age of the Universe. We will revisit this assumption in Sec. 5. We can
then compare the observed value $\Omega_{\text{DM}}h^{2}=0.12$ [36] to the
prediction in order to determine the value of $\lambda_{hs}$ required for
light scalars to constitute all of DM. For the case of dominant production via
Higgs decays we then obtain

$\lambda_{hs}=4.0\times
10^{-10}\left(\frac{m_{s}}{1\,\mathrm{MeV}}\right)^{-1/2}\;.$ (3.22)

These results can now be used as input for studying the further evolution of
the dark sector.

Before doing so let us briefly revisit our assumption that the
$\mathbb{Z}_{2}$ symmetry of the light scalar is already broken before the
EWPT and that we can neglect the temperature dependence of $v_{s}$ and
$m_{s}$. At first sight, one would expect loop diagrams involving the quartic
self-coupling to have important effects at high temperatures, which would
restore the $\mathbb{Z}_{2}$ symmetry for $T\gtrsim v_{s}$ [37]. However,
since the light scalars are never in thermal equilibrium, these contributions
are in fact absent and the only temperature dependence arises from
interactions between the light scalars and SM Higgs bosons. It was shown in
Ref. [10] that for $\lambda_{hs}>0$ these interactions will restore the
$\mathbb{Z}_{2}$ symmetry at high temperatures. Nevertheless, the dominant
production mode for small $v_{s}$ is $h\to ss$, which does not require mixing
between the two Higgs bosons. This production mechanism therefore works in
exactly the same way also if the $\mathbb{Z}_{2}$ symmetry is still unbroken
after the EWPT. The discussion below is therefore independent of whether the
$\mathbb{Z}_{2}$ symmetry breaks before or after the EWPT.777We note that it
was also shown in Ref. [10] that the energy density of domain walls resulting
from the breaking of the $\mathbb{Z}_{2}$ is sufficiently small to be
consistent with observational constraints.

##  4 Evolution of the dark sector

Once the number density of Higgs bosons becomes strongly Boltzmann suppressed,
the freeze-in production of light scalars terminates and the two sectors of
the theory develop completely independently. This does, however, not
necessarily mean that the co-moving number density of light scalars stays
constant. Since the $\mathbb{Z}_{2}$ symmetry of the light scalar is broken,
$2\to 3$ and $3\to 2$ processes can change the number density and temperature
of the dark sector. In particular, these processes can lead to a period of
chemical equilibrium within the dark sector, which ends when the number-
changing processes freeze out. In the following we will first determine the
regions of parameter space where interactions within the dark sector can be
important and then develop the necessary formalism for calculating the final
abundance of light scalars in these regions.

###  4.1 Step 4: Thermalisation

All the production mechanisms discussed above lead to the production of light
scalars with substantial energy. For example, scalars produced in Higgs decays
will have an energy of $m_{h}/2\gg m_{s}$. Scattering between light scalars
will bring them into kinetic equilibrium with each other, such that we can
define a dark sector temperature $T_{\text{dark}}$, which will in general be
different from the temperature $T_{\text{SM}}$ of the SM particles.888In
principle, scattering between light scalars and SM fermions via Higgs exchange
could lead to the exchange of energy and hence kinetic equilibrium between the
two sectors. The cross section for this process is parametrically given by
$\sigma_{\text{scat}}\propto\lambda_{hs}^{2}\,m_{f}^{2}/m_{h}^{4}$. The
corresponding scattering rate $\Gamma_{\text{scat}}=\sigma_{\text{scat}}n_{f}$
is found to be tiny compared to the Hubble rate for the range of
$\lambda_{hs}$ that we are interested in, so it is safe to neglect these
processes. It is however not clear a priori whether the light scalars will
also achieve chemical equilibrium, such that their distribution would follow
an equilibrium distribution with vanishing chemical potential
$\mu_{\text{dark}}$.999Note that chemical equilibrium does in general not
imply that the chemical potentials vanish, but that they add up to zero. In
our case, however, we only consider number-changing processes with final and
initial states of the same particle species, so chemical equilibrium does
imply vanishing chemical potential. Indeed, we will now show that it is
inconsistent to assume $\mu_{\text{dark}}=0$ at the end of freeze-in.

In analogy to the co-moving number density $Y_{s}$ we define a rescaled energy
density

$Z_{s}=\frac{x_{\mathrm{SM}}\,\rho_{s}}{s_{\mathrm{SM}}\,m_{s}}\;.$ (4.1)

For a completely decoupled relativistic dark sector, $\rho_{s}\propto
T_{\mathrm{SM}}\,n_{s}$ and hence $Z_{s}\propto Y_{s}=\text{const}$. If the
dark sector is populated through Higgs decays, on the other hand, each Higgs
decay changes the total number of particles in the dark sector by 2 and the
total energy by $m_{h}$. Thus
$\mathrm{d}\rho_{s}=\tfrac{m_{h}}{2}\mathrm{d}n_{s}$, and the evolution of
$Z_{s}$ is simply given by

$\frac{\mathrm{d}Z_{s}}{\mathrm{d}x_{\mathrm{SM}}}=\frac{m_{h}}{2}\,x_{\mathrm{SM}}\frac{\mathrm{d}Y_{s}}{\mathrm{d}x_{\mathrm{SM}}}\;,$
(4.2)

where $\mathrm{d}Y_{s}/\mathrm{d}x_{\mathrm{SM}}$ is given in eq. (3.17). For
$m_{s}\ll m_{h}$, we then find $Z_{s}/Y_{s}=1.74$ at the end of freeze-in,
which can be rewritten as

$\frac{\rho_{s}}{n_{s}}=1.74\,T_{\mathrm{SM}}\;.$ (4.3)

In other words, freeze-in production of light scalars imposes a specific ratio
of energy density to number density in the dark sector.

Given kinetic equilibrium, the number density and energy density of light
scalars after freeze-in can be written as

$\displaystyle n_{s,\mathrm{fi}}$
$\displaystyle=\frac{1}{2\pi^{2}}\int\frac{k^{2}\,\mathrm{d}k}{\exp\left[(E-\mu_{\text{dark,fi}})/T_{\text{dark,fi}}\right]-1}$
(4.4) $\displaystyle\rho_{s,\mathrm{fi}}$
$\displaystyle=\frac{1}{2\pi^{2}}\int\frac{E\,k^{2}\,\mathrm{d}k}{\exp\left[(E-\mu_{\text{dark,fi}})/T_{\text{dark,fi}}\right]-1}\;,$
(4.5)

where $T_{\text{dark,fi}}$ and $\mu_{\text{dark,fi}}$ denote the temperature
and chemical potential of the dark sector at the end of freeze-in. The precise
definition of when freeze-in ends is of course arbitrary and does not impact
the subsequent calculation. We can now use eq. (4.3) to eliminate
$T_{\text{dark,fi}}$. For example, if we neglect the $-1$ in the denominator
of the expressions above and assume that the light scalars are highly
relativistic, it can be immediately seen that
$\rho_{s,\mathrm{fi}}/n_{s,\mathrm{fi}}=3\,T_{\text{dark,fi}}$ independent of
$\mu_{\text{dark,fi}}$. In combination with eq. (4.3) we therefore find
$T_{\text{dark,fi}}=0.58\,T_{\text{SM,fi}}$. In this approximation, the
temperature of the dark sector is independent of the _amount_ of freeze-in
production, which can therefore only be captured by choosing an appropriate
chemical potential $\mu_{\text{dark}}$.

It should be clear from the discussion above that it will in general only be
possible to match the freeze-in prediction for both $n_{s}$ and $\rho_{s}$ if
both $T_{\text{dark,fi}}$ and $\mu_{\text{dark,fi}}$ are allowed to
vary.101010Refs. [7, 18] suggest to calculate the dark sector temperature by
using the Stefan-Boltzmann law
$T_{\mathrm{dark}}=\left[\rho_{s}\,g^{\ast}/(\rho\,g_{\mathrm{dark}}^{\ast})\right]^{1/4}T_{\mathrm{SM}}$,
which implicitly assumes vanishing chemical potential. While this formula
gives the correct prediction for the energy density, it fails to predict the
correct number density in the dark sector and therefore cannot be applied to
our case. This implies in particular that the number density of light scalars
may initially be quite different from the equilibrium number density expected
for the temperature $T_{\text{dark}}$. Whether or not it stays that way, i.e.
whether or not the light scalars enter into chemical equilibrium, depends on
the rate of number-changing processes.

Assuming for the moment that these processes are inefficient, the subsequent
evolution of the dark sector will be fully determined by the simultaneous
conservation of the dark sector entropy and co-moving number density. The
initial entropy density of the dark sector is given by

$s_{\text{dark,fi}}=\frac{\rho_{s,\mathrm{fi}}+p_{s,\mathrm{fi}}-\mu_{\text{dark,fi}}\,n_{s,\mathrm{fi}}}{T_{\text{dark,fi}}}\;,$
(4.6)

while the co-moving number density is simply given by
$Y_{s}=n_{s,\mathrm{fi}}/s_{\text{SM,fi}}$, because the contribution of the
dark sector to the Hubble expansion rate is completely negligible. Since
entropy is separately conserved in the two sectors, the ratio of the entropy
densities

$\xi=\frac{s_{\text{SM,fi}}}{s_{\text{dark,fi}}}\;,$ (4.7)

is simply a constant.

We can therefore calculate the temperature and chemical potential of the dark
sector, parametrised by $x_{\text{dark}}=m_{s}/T_{\text{dark}}$ and
$\nu_{\text{dark}}=\mu_{\text{dark}}/T_{\text{dark}}$ as a function of the
visible sector temperature $T_{\text{SM}}$. In the absence of number-changing
processes, the evolution of these quantities is fairly trivial. The only
effect that changes the temperature ratio and the chemical potential is when
particles in either of the two sectors become non-relativistic.

Figure 6: Number-changing processes in the dark sector which result in a
decrease in the dark sector temperature and an increase in the dark matter
number density.

We can use these observations to answer the question whether it is consistent
to ignore the effect of number-changing processes. The most important such
processes are the $2\to 3$ processes shown in figure 6, which convert kinetic
energy into new particles, such that the dark sector temperature decreases
while the co-moving number density increases. For $m_{s}\ll\sqrt{S}$ the cross
section for the simplest of these processes (shown in the leftmost panel) is
given by111111The cross section is in fact proportional to
$\lambda_{s}^{4}\,v_{s}^{2}$, which we write as
$2\,\lambda_{s}^{3}\,m_{s}^{2}$.

$\sigma_{2\to 3}=\frac{27\,\lambda_{s}^{3}m_{s}^{2}}{64\pi^{3}S^{2}}\;.$ (4.8)

We find however that the additional diagrams give a larger contribution to the
total cross section, because the intermediate particles can be nearly on-
shell. For these diagrams it becomes unfeasible to perform the phase-space
integration analytically (in particular in the important limit $\sqrt{S}\to
3\,m_{s}$), so that we use calchep 3.7.1 [38] to obtain the full result for
arbitrary centre-of-mass energies.

The corresponding interaction rate is given by $\Gamma_{2\to 3}=\langle\sigma
v\rangle_{2\to 3}\,n_{s}$, where $\langle\sigma v\rangle_{2\to 3}$ denotes the
thermal average _with respect to the temperature of the dark sector_ , which
can be calculated in terms of the SM temperature. If $\Gamma_{2\to
3}(T_{\text{SM}})\ll H(T_{\text{SM}})$ for all temperatures $T_{\text{SM}}$,
the light scalars will never enter into chemical equilibrium with each other.
In this case, we can indeed treat the co-moving number density $Y_{s}$ as
constant and obtain the final relic abundance of light scalars directly from
the freeze-in yield.

Figure 7: Ratio of the rate of number-changing processes, $\Gamma_{2\to
3}=\langle\sigma_{2\to 3}\rangle\,n_{s}$ and the Hubble expansion rate $H$ as
a function of the inverse SM temperature $x_{\mathrm{SM}}$. We consider fixed
$m_{s}=1\,\mathrm{MeV}$ and several different values of $\lambda_{s}$. A
second $x$-axis indicates the corresponding inverse temperature of the dark
sector $x_{\text{dark}}$ under the assumption that number-changing processes
are negligible. Initially, this temperature is given by
$x_{\text{dark}}\approx 1.72\,x_{\mathrm{SM}}$ but then evolves independently
due to the separate conservation of entropy in the two sectors.

We show in figure 7 the ratio $\Gamma_{2\to 3}/H$ for specific values of
$\lambda_{s}$ and $m_{s}$ under the assumption that the co-moving number
density is given by
$Y_{s}=\Omega_{\text{DM}}\,\rho_{\text{c,0}}/(m_{s}\,s_{\text{SM,0}})=\text{const}$,
such that light scalars would constitute all of the DM in the present
Universe. The two $x$-axis of the figure indicate $x_{\text{SM}}$ and
$x_{\text{dark}}$, respectively. We observe that the interaction rate is
largest relative to the Hubble rate for $x_{\text{dark}}\approx 0.4$, i.e.
shortly before the average centre-of-mass energy becomes insufficient to
produce a new particle. Since the cross section is simply proportional to
$\lambda_{s}^{3}$ for fixed $m_{s}$, we can immediately infer the largest
value of $\lambda_{s}$ for which $\Gamma_{2\to 3}/H<1$ for all temperatures.
For example, for $m_{s}=1\,\text{MeV}$ we find $\lambda_{s}\lesssim 3\cdot
10^{-5}$.

Figure 8: Regions of different dark sector evolution as a function of $m_{s}$
and $\lambda_{s}$. In the red shaded regions the reaction rate $\langle\sigma
v\rangle_{2\to 3}\,n_{s}$ exceeds the Hubble rate $H(T)$ at some point in the
early Universe, leading to a thermalised dark sector with vanishing chemical
potential. In the yellow shaded regions, the ratio of the rates is between 1%
and unity, such that number-changing processes cannot be neglected even though
the dark sector might not reach chemical equilibrium. In the green shaded
regions, the reaction rate never exceeds 1% of the expansion rate, such that
the abundance will not be significantly affected by number-changing processes.
The grey area in the bottom-right corner indicate the parameter range where
$2\to 2$ processes give the dominant contribution to the freeze-in yield.

We illustrate in figure 8 how the maximum value of $\Gamma_{2\to 3}/H$ depends
on the scalar mass $m_{s}$ and its self coupling $\lambda_{s}$ under the
assumption that $Y_{s}$ is constant and corresponds to the observed DM relic
abundance. The five stars correspond to the five curves shown in figure 7. In
the red shaded parameter region the maximum ratio exceeds unity, meaning that
$2\to 3$ processes will be efficient enough to bring the dark sector into
chemical equilibrium. In the green shaded region, on the other hand, the ratio
never exceeds 0.01 and we can safely neglect number-changing processes. In the
intermediate region, shaded in yellow, the reaction rate stays below the
Hubble expansion rate, but is still large enough that we cannot treat $Y_{s}$
as constant. We will deal with this case together with the case of full
chemical equilibrium in the following section.

Finally, we also indicate in figure 8 the parameter region where the freeze-in
abundance is dominantly set by $t$-channel processes (rather than Higgs
decays). We find that in these parameter regions $2\to 3$ processes are
completely negligible, so that the evolution of the dark sector after the end
of freeze-in is trivial. Conversely, in the regions of parameter space where
chemical equilibrium may be established, it is fully justified to estimate the
energy density of the dark sector under the assumption of freeze-in via Higgs
decays, as we have done above.

###  4.2 Step 5: Dark sector freeze-out

Let us now take a closer look at what happens when the dark sector reaches
chemical equilibrium. In this case $\mu_{\text{dark}}=0$ and the number
density is given by the equilibrium distribution

$n^{\text{eq}}_{s}=\frac{1}{2\pi^{2}}\int\frac{k^{2}\,\mathrm{d}k}{\exp\left(E/T_{\text{dark}}\right)-1}\;.$
(4.9)

As the light scalars become non-relativistic, the equilibrium number density
becomes exponentially suppressed. The reason is that the $2\to 3$ processes,
which helped to populate the dark sector, now become inefficient due to the
lack of sufficient kinetic energy to produce a third particle. The inverse
process, on the other hand, remains fully efficient and leads to a depletion
of light scalars in combination with an increase of the dark sector
temperature. This is the well-known cannibalism mechanism [19, 22, 21]. It
ceases to be efficient once the rate for the $3\to 2$ process drops below the
Hubble expansion rate, at which point the interactions freeze out and the co-
moving number density of light scalars becomes constant. The final abundance
of light scalars is then determined by the temperature when the departure from
chemical equilibrium happens.

While it is possible to use entropy conservation to obtain an approximate
estimate of the freeze-out abundance [19], we will calculate the evolution of
$Y_{s}$ numerically by solving the Boltzmann equation

$\frac{\mathrm{d}{Y_{s}}}{\mathrm{d}x_{\text{SM}}}=\frac{s_{\text{SM}}}{H\,x_{\text{SM}}}\langle\sigma
v\rangle_{2\to 3}Y_{s}^{2}\left[1-\frac{Y_{s}}{Y_{s}^{\text{eq}}}\right]\;,$
(4.10)

where $Y^{\text{eq}}_{s}=n^{\text{eq}}_{s}/s_{\text{SM}}$. To obtain this
equation we have made use of the fact that

$\langle\sigma v\rangle_{2\to 3}\,n^{\text{eq}}_{s}=\langle\sigma
v^{2}\rangle_{3\to 2}\,(n^{\text{eq}}_{s})^{2}\;,$ (4.11)

because the rates of the two processes must be equal in thermal equilibrium
(see appendix A of Ref. [20]). The Boltzmann equation can now be solved for
given $Y_{s,\mathrm{fi}}$ and $\langle\sigma v\rangle_{2\to 3}$ to yield the
present-day value of $Y_{s}$, $Y_{s,0}$. To do so, it is convenient to rewrite
eq. (4.11) as

$\frac{\mathrm{d}\log Y_{s}}{\mathrm{d}\log
x_{\text{SM}}}=\frac{s_{\text{SM}}}{H}\langle\sigma v\rangle_{2\to
3}Y_{s}\left[1-\frac{Y_{s}}{Y^{\text{eq}}_{s}}\right]\;.$ (4.12)

Note that our approach differs from the one taken in Refs. [20, 24] in that we
base our calculation on $\langle\sigma v\rangle_{2\to 3}$ rather than
$\langle\sigma v^{2}\rangle_{3\to 2}$, which we can calculate for arbitrary
centre-of-mass energies (whereas Refs. [20, 24] only consider the non-
relativistic limit). This makes it possible to treat thermalisation and
freeze-out in a unified framework, i.e. we can use the Boltzmann equation also
in the case where the dark sector never reaches chemical equilibrium.

The disadvantage of our approach is that we need to track the dark sector
temperature throughout the evolution of the dark sector (whereas
$\langle\sigma v^{2}\rangle_{3\to 2}$ becomes independent of the dark sector
temperature in the non-relativistic limit). To do so, we again make use of the
fact that the entropy ratio $\xi$ remains constant, such that we can write

$\frac{1}{Y_{s}}=\frac{\xi\,s_{\text{dark}}}{n_{s}}=\xi\left(x_{\text{dark}}\frac{\rho_{s}+p_{s}}{m_{s}\,n_{s}}-\nu_{\text{dark}}\right)\;,$
(4.13)

where as before $x_{\text{dark}}=m_{s}/T_{\text{dark}}$ and
$\nu_{\text{dark}}=\mu_{\text{dark}}/T_{\text{dark}}$. If we neglect the $-1$
in the denominator of the Bose-Einstein distribution, which is a good
approximation as long as the dark sector is sufficiently sparsely populated,
we find that

$\zeta\equiv\frac{\rho_{s}+p_{s}}{m_{s}\,n_{s}}$ (4.14)

depends only on $x_{\text{dark}}$, but not on $\nu_{\text{dark}}$, which in
turn is given by

$\nu_{\text{dark}}=\log\left(\frac{Y_{s}}{Y^{\text{eq}}_{s}}\right)\;.$ (4.15)

Combining these equations gives

$\frac{1}{Y_{s}}=\xi\left[x_{\text{dark}}\,\zeta(x_{\text{dark}})-\log\left(\frac{Y_{s}}{Y^{\text{eq}}_{s}(x_{\text{dark}})}\right)\right]\;,$
(4.16)

which is now independent of $\nu_{\mathrm{dark}}$. Eq. (4.16) can be solved
numerically for $x_{\text{dark}}$ for given $Y_{s}$, $\xi$ and $s_{\text{SM}}$
(the latter entering in the calculation of $Y^{\text{eq}}_{s}$).

Figure 9:  Evolution of $Y_{s}$ as a function of $x_{\mathrm{SM}}$ for fixed
scalar mass $m_{s}=1\,\text{MeV}$ and different values of the coupling
$\lambda_{s}$. The initial value for $Y_{s}$ is chosen such that in the
absence of number-changing processes the light scalars would constitute all of
DM in the present Universe.

We are now in the position to solve eq. (4.11) numerically (see appendix A for
details). Figure 9 shows a number of examples for how $Y_{s}$ evolves as
function of $x_{\text{SM}}$. All examples correspond to the same scalar mass,
$m_{s}=1\,\mathrm{MeV}$, but different values of $\lambda_{s}$ (as indicated
in figure 8). Note that $x_{\text{dark}}$ depends on $\lambda_{s}$ and hence
we do not show a second $x$-axis in this plot. We identify the following cases
(cf. figure 8):

  * •

For $\lambda_{s}\lesssim 5\cdot 10^{-6}$ the reaction rate always stays well
below the Hubble rate and the co-moving number density does not change.

  * •

For $\lambda_{s}\sim 10^{-5}$ the reaction rate stays only slightly below the
Hubble rate. The $2\to 3$ processes can therefore not be entirely neglected
and lead to a slight increase in $Y_{s}$. Nevertheless, the dark sector does
not reach chemical equilibrium in this case.

  * •

For $\lambda_{s}\gtrsim 3\cdot 10^{-5}$ the dark sector reaches chemical
equilibrium and subsequently freezes out. With increasing $\lambda_{s}$
equilibrium is reached more quickly and is maintained longer before freeze-
out.

As for standard freeze-out we therefore find that larger couplings lead to a
smaller relic abundance. For the case of freeze-out via $3\to 2$ processes,
however, the dependence on the couplings is much milder. Indeed, as shown in
appendix A, the co-moving number density after dark sector freeze-out can
approximately be written as

$Y_{s}=\frac{4}{3\,\xi\left[\log\left(\frac{M_{\mathrm{p}}\,\lambda_{s}^{3}}{m_{s}\,\xi^{2/3}}\right)-a\right]}$
(4.17)

with $a\approx 4\text{--}5$. In particular, $Y_{s}$ depends only
logarithmically on $\lambda_{s}$ and therefore changes only slightly even if
$\lambda_{s}$ varies substantially.

###  4.3 Relic abundance of light scalars

Figure 10:  Left: Value of $\lambda_{hs}$ required to obtain the correct
abundance from freeze-in, such that at the end of dark sector freeze-out the
light scalars account for all of DM. Right: Ratio of the scalar abundance at
the end of freeze-in, $Y_{s,\text{fi}}$, and the scalar abundance today,
$Y_{s,0}$, under the assumption that light scalars constitute all of DM. One
can clearly identify parameter regions where $2\to 3$ processes lead to an
increase in the abundance ($Y_{s,0}>Y_{s,\text{fi}}$) and regions where dark
sector freeze-out leads to a depletion of light scalars
($Y_{s,0}<Y_{s,\text{fi}}$).

We now have all the necessary ingredients to combine the calculations outlined
in steps 1–5 above and calculate the final relic abundance of light scalars as
a function of $m_{s}$, $\lambda_{s}$ and $\lambda_{hs}$. We show in the left
panel of figure 10 the value of $\lambda_{hs}$ needed to obtain
$\Omega_{s}h^{2}=0.12$ as a function of $m_{s}$ and $\lambda_{s}$. Following
eq. (3.22), one would expect $\Omega_{s}$ to be independent of $\lambda_{s}$
in the absence of $2\to 3$ processes, such that lines of constant
$\lambda_{hs}$ would simply be vertical. This is indeed what we find for small
values of $\lambda_{s}$, when number-changing processes are inefficient.

For larger self-couplings, however, the behaviour becomes more complicated and
lines of constant $\lambda_{hs}$ are no longer vertical. This effect can be
most easily understood by consulting the right panel of figure 10, which shows
the ratio of the scalar abundance after freeze-in, $Y_{s,\text{fi}}$, and the
present-day abundance after dark sector freeze-out, $Y_{s,0}$, for parameter
points that satisfy the relic density requirement. In the bottom-right corner
we see once again that number-changing processes have a negligible effect on
the abundance of light scalars, such that $Y_{s,\text{fi}}/Y_{s,0}=1$. For
somewhat larger values of $\lambda_{s}$ (indicated by dark blue shading) the
dark sector still does not reach chemical equilibrium, but the effect of $2\to
3$ processes cannot be completely neglected, leading to a net increase in the
co-moving number density (and a decrease the dark sector temperature relative
to the SM temperature). This means that $Y_{s,\text{fi}}/Y_{s,0}<1$ and hence
smaller freeze-in yields (and smaller values of $\lambda_{hs}$) are sufficient
to reproduce the observed DM abundance in the present Universe. For even
larger values of $\lambda_{s}$ (indicated by orange and red shading) the dark
sector enters into chemical equilibrium and dark sector freeze-out leads to a
decrease in the co-moving number density, as well as to heating of the dark
sector. Since in this parameter region $Y_{s,\text{fi}}/Y_{s,0}>1$, larger
freeze-in yields (and hence larger values of $\lambda_{hs}$) are required to
reproduce the observed DM relic abundance.

Finally, we note that in the bottom-right corner of the left panel of figure
10 the freeze-in contribution of additional channels beyond Higgs decays
becomes important (see section 3.3), leading to an additional dependence of
$Y_{s}$ on $\lambda_{s}$, such that lines of constant $\lambda_{hs}$ are no
longer vertical, even though the effect of number-changing processes is
negligible.

##  5 Phenomenological consequences

Let us now answer the question whether the light scalars considered above can
constitute all of the DM in the present Universe. For this purpose we need to
consider their lifetime (which must not only exceed the age of the Universe
but also satisfy various constraints from indirect detection) as well as the
self-interaction cross section (which is constrained by a number of
astrophysical observations). A similar discussion can be found in Ref. [10].

###  5.1 Light scalar decays

In the present work we focus on scalar DM particles in the mass range between
1 keV and 100 MeV, such that there are only two possible decay modes: $s\to
e^{+}e^{-}$ and $s\to\gamma\gamma$ (see also Ref. [39] for a similar
discussion in the context of relaxion-Higgs mixing). The corresponding decay
widths can be calculated in complete analogy to the case of a light Higgs
boson [40, 41]. The electronic partial decay width is given by

$\Gamma(s\to
e^{+}e^{-})=\sin^{2}\theta\,\frac{m_{s}\,m_{e}^{2}}{8\pi\,v^{2}}(1-z_{e})^{3/2}$
(5.1)

with $z_{e}=4m_{e}^{2}/m_{s}^{2}$ and $\Gamma(s\to e^{+}e^{-})=0$ for
$z_{e}\geq 1$. Decays into photons are absent at tree-level, but are induced
via loops of SM fermions and bosons. The resulting partial width can be
written as

$\Gamma(s\to\gamma\gamma)=\sin^{2}\theta\,\frac{\alpha^{2}\,m_{s}^{3}}{256\pi^{3}\,v^{2}}\left|f(z_{e})+\frac{7}{3}\right|^{2}$
(5.2)

where $\alpha$ is the fine-structure constant and

$f(\tau)=2\tau\left[1+(1-\tau)\,\text{arctan}^{2}\left(\frac{1}{\sqrt{\tau-1}}\right)\right]$
(5.3)

is the form factor for fermionic loops. We include this form factor only for
electrons and evaluate all other loops in the limit $z_{f}\to\infty$, which is
a good approximation even for up and down quarks, because their effective mass
in the loop is of the order of the pion mass [42].

As expected, we find that decays into electrons completely dominate once they
become kinematically allowed. In fact, for the typical coupling strengths
required for freeze-in production ($\sin\theta\sim
10^{-12}\text{--}10^{-14}$), the resulting lifetime is well below
$10^{20}\,\mathrm{s}$. Such lifetimes are strongly excluded by observations of
the Cosmic Microwave Background [36, 43, 44] and, more recently, measurements
of the 21cm radiation temperature in the re-ionisation epoch [45, 46]. In
other words, only a small fraction of DM could be in the form of unstable MeV
scalars produced via the freeze-in mechanism.

For $m_{s}<1\,\mathrm{MeV}$ only the photonic decay mode remains kinematically
allowed and the scalar lifetimes become much larger. At the same time,
however, the observational signature is much more striking: a mono-energetic
$\gamma$-ray or x-ray line with energy $E_{\gamma}=m_{s}/2$. Such lines have
been searched for in the context of keV sterile neutrinos [47]. To convert a
bound on the active-sterile mixing angle $\sin^{2}2\theta$ into a bound on the
scalar lifetime, we note that

$\Gamma_{N}=\frac{9\,\alpha\,G_{\mathrm{F}}^{2}}{2048\pi^{4}}\sin^{2}2\theta\,m_{s}^{5}$
(5.4)

and that the bound on the scalar lifetime is stronger by a factor of 2,
because each scalar decay produces two photons rather than one. The strongest
bounds come from INTEGRAL [48] for $m_{s}\sim 100\text{--}1000\,\mathrm{keV}$,
from NuSTAR [49] for $m_{s}\sim 10\text{--}100\,\mathrm{keV}$ and from a
combination of x-ray observations of M31 [50] for even smaller masses. These
constraints typically require $\tau\gtrsim
10^{27}\text{--}10^{29}\,\mathrm{s}$.

###  5.2 Self-interactions and structure formation

The quartic self-coupling between the light scalars gives rise to velocity-
independent self-interactions. The corresponding cross section is given by
[18]

$\frac{\sigma_{s}}{m_{s}}=\frac{9\,\lambda_{s}^{2}}{32\pi\,m_{s}^{3}}\;.$
(5.5)

Observations from the Bullet Cluster constrain this cross section to be
smaller than approximately $1\,\mathrm{cm^{2}\,g^{-1}}$ [51, 52, 53], which
translates to

$\lambda_{s}\lesssim 0.007\left(\frac{m_{s}}{1\,\mathrm{MeV}}\right)^{3/2}\;.$
(5.6)

Since we are considering DM particles with mass in the keV range, another
potential concern are constraints from structure formation, in particular from
the Lyman-$\alpha$ forest [54, 55]. The impact of our model on structure
formation is however very different from the case of sterile neutrinos. The
reason is that self-interactions between DM particles prevent free-streaming,
and structures can therefore only be washed out by diffusion.121212Different
ways to evade structure formation constraints for keV-scale DM are discussed
in Refs. [29, 56]. The diffusion length $l_{s}$ in the presence of self-
interactions is given by [57]

$l_{s}^{2}=\int_{0}^{t_{\mathrm{dec}}}\frac{\mathrm{d}t\,\langle
v_{s}\rangle^{2}}{a^{2}\,n_{s}\,\langle\sigma_{s}v_{s}\rangle}\;.$ (5.7)

The integration range in principle extends to the time $t_{\mathrm{dec}}$ when
self-interactions become inefficient, i.e. when
$n_{s}\langle\sigma_{s}v_{s}\rangle<H$, but the main contribution to the
diffusion stems from relativistic particles. We can therefore approximate

$l_{s}^{2}\approx\int_{0}^{a_{\mathrm{NR}}}\frac{\mathrm{d}a}{H\,a^{3}\,n_{s}\,\langle\sigma_{s}v_{s}\rangle}\;,$
(5.8)

where $a_{\mathrm{NR}}$ denotes the scale factor when $x_{\mathrm{dark}}\sim
1$.

If self-interactions (and hence $3\to 2$ processes) are extremely weak, we can
treat the temperature ratio $\alpha\equiv T_{\mathrm{dark}}/T_{\mathrm{SM}}$
as constant, with values typically in the range $\alpha\approx
0.1\text{--}0.3$. Using $n_{s}=Y_{0}\,s_{\mathrm{SM}}$ and

$\langle\sigma_{s}v_{s}\rangle=\frac{9\,\lambda_{s}^{2}}{64\pi\,T_{\mathrm{dark}}^{2}}$
(5.9)

for $T_{\mathrm{dark}}\gg m_{s}$, it is easy to see that the integrand of Eq.
(5.8) is independent of $a$ during radiation domination (ignoring the slight
temperature dependence of $g_{\mathrm{SM}}^{\ast}$) and hence

$l_{s}^{2}\propto\frac{a_{\mathrm{NR}}\,\alpha^{2}}{Y_{0}\,\lambda_{s}^{2}}\;.$
(5.10)

Using the same approximation as above, $a_{\mathrm{NR}}$ is defined by $\alpha
T_{\mathrm{SM}}=m_{s}$ and hence $a_{\mathrm{NR}}\propto\alpha/m_{s}$. Since
$Y_{0}\propto 1/m_{s}$, we therefore find that $l_{s}$ is largely independent
of $m_{s}$. This leads to the extremely simple result

$l_{s}\sim 10^{-11}\,\mathrm{Mpc}\frac{\alpha^{3/2}}{\lambda_{s}}\;.$ (5.11)

Hence, for $\lambda_{s}>10^{-11}$ the matter power spectrum remains unaffected
on observable scales.

We note however that $a_{\text{dec}}$ can potentially be significantly larger
than $a_{\text{NR}}$, such that dark matter behaves like a non-relativistic
collisional fluid until rather late times. For example, for
$m_{s}=10\,\mathrm{keV}$ and $\lambda_{s}=10^{-6}$ we find
$a_{\text{dec}}\approx 2\times 10^{-5}$ (corresponding to
$T_{\text{SM,dec}}\approx 10\,\mathrm{eV}$). Such large values of
$a_{\text{dec}}$ may affect the growth of density perturbations and hence
potentially lead to modifications of the CMB [58]. We leave a study of these
effects, which may strengthen the bound on the self-interaction cross section,
to future work.

###  5.3 Results

Figure 11:  Constraints on decaying lights scalars from a combination x-ray
observations (orange) and from the bound on self-interactions obtained from
the Bullet Cluster (blue). The dashed blue line indicates a self-interaction
cross section one order of magnitude below current bounds, while the region
shaded in red corresponds to the parameter region preferred by the claimed
excesses of x-ray emission around 3.5 keV.

The results of our analysis are summarised in figure 11. For each value of
$\lambda_{s}$ and $m_{s}$ we fix $\lambda_{hs}$ such that the light scalars
constitute all of the DM in the present Universe (see figure 10). Across the
entire plot the lifetime of the scalars exceeds the age of the Universe by
many orders of magnitude, but for large values of $m_{s}$, as well as for
small values of $\lambda_{s}$ (corresponding to large values of $v_{s}$ and
hence to larger mixing with the SM Higgs boson) there are strong constraints
from searches for x-ray lines. For small masses and large self-couplings, on
the other hand, constraints from DM self-interactions become strong and
exclude a large part of the parameter space. The two constraints are highly
complementary and in combination rule out light scalars with masses above
about 10 keV. For smaller scalar masses, however, viable parameter regions
remain.

This mass range is of particular interest in the context of interpreting a
number of claimed observations of an anomalous x-ray emission around 3.5 keV
in several galaxy clusters [26, 27]. We indicate the parameter region where
these excesses can be accommodated by a red box, noting the well-known tension
between these claims and the constraints obtained from M31 [59]. Intriguingly,
this parameter region corresponds to self-interaction cross sections of the
order of $0.1\,\mathrm{cm^{2}\,g^{-1}}$ (indicated by the dashed blue line),
which are only an order of magnitude below current bounds and may be testable
with novel approaches, such as measurements of core sizes in galaxy clusters
[60, 61].

In conclusion, should the x-ray line emission around 3.5 keV be confirmed by
future observations, light scalars produced via the freeze-in mechanism may
provide an attractive explanation. In contrast to sterile neutrinos, these
particles do not free-stream in the early Universe, but the effects of self-
interactions may play an important role. Light scalars therefore have a
distinct phenomenology that can be explored with an improved understanding of
structure formation at small scales.

##  6 Conclusions

In this paper we have studied in detail the cosmological evolution and
phenomenological properties of light scalar bosons in the keV to MeV range
that mix with the SM Higgs boson. While the interactions with the Higgs render
the light scalars unstable, for sufficiently small couplings their lifetime
can be large compared to the age of the Universe. In these parameter regions
light scalars never enter into thermal equilibrium with SM particles, so that
their abundance is set via the freeze-in mechanism. Nevertheless, self-
interactions and number-changing processes may be sufficiently large that the
light scalars thermalise amongst themselves, leading to a freeze-out in the
dark sector.

A close investigation of the freeze-in production of light scalars reveals a
number of subtleties not fully appreciated previously. First of all, we point
out the importance of the electroweak phase transition, which dramatically
changes the production channels and – crucially – renders the freeze-in yield
insensitive to the reheating temperature. Furthermore, we discuss the effects
of a temperature-dependent Higgs mass and vev, which can potentially lead to
the production of light scalars through oscillations (although the resulting
yield is found to be sub-dominant in practice). Finally, we have studied in
detail the impact of a proper treatment of quantum statistics for the
particles in the thermal bath and find substantial differences to the case of
classical statistics.

Due to the existence of number-changing processes, the evolution of the dark
sector does not terminate at the end of freeze-in. We calculate the rate of
$2\to 3$ processes in order to determine whether the dark sector reaches
chemical equilibrium. In large regions of parameter space these processes play
an important role and we solve the corresponding Boltzmann equation in order
to calculate the abundance of light scalars after dark sector freeze-out. For
the first time we consistently include the evolution of the dark sector
temperature and chemical potential in the Boltzmann equation in order to
obtain accurate results even when the dark sector does not fully thermalise or
when particles are still semi-relativistic during dark sector freeze-out.

We find that it is possible for light scalars in the keV to MeV range to
account for all of DM, provided the Higgs portal coupling $\lambda_{hs}$ is of
the order of $10^{-10}\text{--}10^{-8}$. However, the combined constraints
from the CMB, 21cm astronomy, x-ray observations and bounds on the self-
scattering cross section require such light scalar DM particles to have a mass
of less than about 10 keV (see figure 11). Excitingly, this is precisely the
mass range where light scalars can provide an explanation of the anomalous
x-ray emission at 3.5 keV observed in various astrophysical systems. In
contrast to sterile neutrinos, light scalars do not free-stream and therefore
do not suppress structures on small scales, but their self-interactions may
lead to other observable effects such as the formation of cores in galaxy
clusters.

Finally, we note that large regions of parameter space opens up if we do not
require light scalars to constitute all of DM. In particular, there are
effectively no bounds from DM self-interactions on a sub-dominant fraction of
light scalars. Moreover, it would be interesting to explore the case where the
light scalars have a lifetime that is small compared to the present Universe,
such that they can modify the reionisation history. Such light scalars could
for example act as the mediator between the SM and another stable DM particle,
which also obtains its abundance via the freeze-in mechanism. A detailed study
of these possibilities will be left for future work.

## Acknowledgments

We thank Camilo Garcia-Cely, Sebastian Wild and Bryan Zaldivar for valuable
comments on the manuscript, Maria Archidiacono, Basudeb Dasgupta, Andreas
Goudelis, Jan Heisig, Julien Lesgourgues, Kai Schmidt-Hoberg and Tommi
Tenkanen for enlightening discussions and the authors of micrOmegas for
technical support. This work is funded by the DFG Emmy Noether Grant No. KA
4662/1-1. Numerical calculations were performed with computing resources
granted by RWTH Aachen University under project rwth0280.

##  Appendix A Solving the Boltzmann equation

In this appendix we discuss how to solve the Boltzmann equation for $2\to 3$
processes both numerically and analytically, and provide details on the
implementation that we have employed to obtain our results.

Unlike in the case of standard thermal production of dark matter, where we can
identify $x_{\mathrm{SM}}=x_{\mathrm{dark}}$, the inverse temperature
parameter of a decoupled dark sector $x_{\mathrm{dark}}$ has a non-trivial
dependence on $x_{\mathrm{SM}}$, which is given by eq. (4.16). Since in our
set-up the light scalars are relativistic after freeze-in, $x_{\mathrm{dark}}$
can only be determined by solving this equation numerically.131313We note in
passing that in the non-relativistic limit eq. (4.16) can be solved for
$x_{\mathrm{dark}}$, giving
$x_{\mathrm{dark}}=\frac{1}{2\pi}\,\exp{\left(\frac{5}{3}-\frac{2}{\xi
Y_{s}}\right)}\,m_{s}^{2}\,\left(Y_{s}s_{\text{SM}}\right)^{-\frac{2}{3}}.$ We
find, however, that this approximation is not sufficiently accurate for our
purposes and therefore always use the exact numerical solution. In our
numerical implementation to solve the Boltzmann equation for $2\to 3$
processes, this step turns out to be the bottleneck.

Furthermore the structure of eq. (4.10) corresponds to a so-called “stiff”
differential equation. For this class of differential equations ODE-solving-
algorithms without adaptive step size will fail to obtain the correct
solution. It is therefore preferable to use algorithms that adapt their step
size depending on how quickly the differential equation changes. For points
where the right-hand side of the differential equation changes rapidly if the
input parameters are varied, the chosen step size will be smaller than for
points where the right-hand side of the differential equation is nearly
invariant. To choose the appropriate step size, such adaptive algorithms need
to evaluate the Jacobian of the right hand side of the differential equation,
which in our case requires substantial computational expense.141414Note that
both $\langle\sigma v\rangle_{2\to 3}$ and $Y_{s}^{\text{eq}}$ depend
implicitly on $x_{\text{dark}}$, which in turn depends on $Y_{s}$ via the
numerical solution of eq. (4.16). It is hence necessary to determine the
Jacobian numerically, which requires repeated calls of the function to
determine $x_{\mathrm{dark}}$ for given inputs $\log{x_{\mathrm{SM}}}$ and
$\log{Y_{s}}$. We therefore need to find a good trade-off between stability of
the solution and the accuracy of the step size.

In our analysis we solve eq. (4.10) using the implementation of the multistep-
backward-differentiation method of the GSL-v1.15-library [62] which we found
to give robust results with relatively few computational expensive steps.

An important simplification is possible whenever we are not interested in
tracking the full solution of the Boltzmann equation, but only in finding the
final value of $Y_{s}$. Indeed, if the coupling is large enough to thermalise
the dark sector, it is usually not important when exactly chemical equilibrium
is reached but it suffices to know that $Y_{s}$ follows $Y_{s}^{\mathrm{eq}}$
until the number-changing processes decouple. Indeed, such large couplings are
computationally particularly expensive, so it is desireable to have an
improved treatment. For this we adapt the method implemented in MadDM [63] for
the case of thermal freeze-out. Whenever it is safe to assume that chemical
equilibrium is reached at some point, we make an initial guess of the value
for $x_{\mathrm{SM,start}}$ when number-changing processes decouple
($\Gamma_{2\to 3}\lesssim H$). We then start solving the Boltzmann equation at
this point assuming
$Y_{s}(x_{\mathrm{SM,start}})=Y_{s}^{\text{eq}}(x_{\mathrm{SM,start}})$ and
evolve to the regime where $Y_{s}$ approaches a constant value. To determine
the quality of our initial guess of $x_{\mathrm{SM,start}}$, we then move
$x_{\mathrm{SM,start}}$ to a smaller value (i.e. to higher temperatures) and
repeat the evolution. This procedure is repeated until the final value of
$Y_{s}$ converges. This way we avoid solving the Boltzmann equation in the
regime where $\Gamma_{2\to 3}\gg H$ and hence tiny step sizes would be needed
for an accurate solution.

Figure 12:  Comparison of the final result $Y_{s}^{0}$ of both evolution
methods for fixed mass as a function of the dark scalar self-coupling. The
colour code of the background is the same as in figure 8. For large couplings
(when $\max\left(\Gamma/H\right)\gtrsim 1$) both methods agree, whereas for
small couplings the method based on the assumption of chemical equilibrium
leads results that are too large.

In figure 12 we compare the results obtained from both methods for a fixed
mass of $m_{s}=1\,\mathrm{MeV}$ as a function of the coupling $\lambda_{s}$.
We use the same background colours as in figure 8 to indicate the naive ratio
of $\Gamma_{2\to 3}/H$. We see that both methods lead to the same results when
the coupling is high enough (in the red regime) such that chemical equilibrium
is maintained. For small couplings (green and yellow regimes) the results
differ, as the assumption of chemical equilibrium, which is essential for the
iterative method, is no longer satisfied. In this figure we can also see that
the transition from the yellow to the red regime serves as a smooth transition
between the two methods. In all our calcualtions we therefore use the
iterative method for parameter point in the red regime and the full evolution
for points in the yellow and green regime.

Finally, we note that under the assumption that DM is sufficiently non-
relativistic during dark sector freeze-out, it is also possible to derive an
approximate analytical estimate of the relic abundance after dark sector
freeze-out. As always the condition for freeze-out is that

$H=n_{s,\mathrm{fo}}\,\langle\sigma v\rangle_{2\to
3}=Y_{s,\mathrm{fo}}\,s_{\mathrm{SM}}\,\langle\sigma v\rangle_{2\to 3}\;.$
(A.1)

In addition to the quantity of interest, i.e. the co-moving number density
$Y_{s,\mathrm{fo}}=Y_{s,0}$ at freeze-out, this equation also depends on two
further unknowns. The temperature of the visible sector at freeze-out,
$T_{\mathrm{SM,fo}}$, enters via the Hubble rate and the SM entropy density,
and the freeze-out temperature of the dark sector, $T_{\mathrm{dark,fo}}$,
enters through the thermally averaged cross section. Indeed, in the non-
relativistic regime, this cross section be approximately written as

$\langle\sigma v\rangle_{2\to 3}\propto
e^{-x_{\mathrm{dark}}}\frac{\lambda_{s}^{3}}{m_{s}^{2}}\;,$ (A.2)

where the exponential reflects the fraction of DM particles with enough energy
to induce a $2\to 3$ process and the constant of proportionality needs to be
determined numerically.

To eliminate $T_{\mathrm{dark,fo}}$, we can make use of the fact that in
chemical equilibrium and for $x_{\mathrm{dark}}\gg 1$ eq. (4.13) simplifies to
[19]

$Y_{s}=\frac{1}{\xi\,x_{\mathrm{dark}}}\;.$ (A.3)

This leads to the somewhat surprising result that the right-hand side of eq.
(A.1) depends not linearly but exponentially on $Y_{s}$, in striking contrast
to the case of standard freeze-out. As a result, we expect the DM abundance
when freeze-out happens to depend only logarithmically on the fundamental
model parameters. To make this explicit, we need to eliminate
$T_{\mathrm{SM,fo}}$ by making use of the fact that $\xi=\mathrm{const}$,
which implies

$x_{\mathrm{SM}}^{3}\propto
e^{x_{\mathrm{dark}}}\frac{g^{\ast}_{\mathrm{SM}}\,\sqrt{x_{\mathrm{dark}}}}{\xi}$
(A.4)

during radiation domination.

Combining the above equations, one finds the approximate solution

$Y_{s}=\frac{4}{3\,\xi\left[\log\left(\frac{M_{\mathrm{p}}\,\lambda_{s}^{3}}{m_{s}\,\xi^{2/3}}\right)-a\right]}\;,$
(A.5)

where the parameter $a$ includes all the numerical factors not explictly
included in the first term. By fitting to our results, we determine $a\approx
4\text{--}5$ with a slight dependence on both $x_{\text{dark}}$ and
$x_{\text{SM}}$. This expression is found to give a good fit to the numerical
solution of the Boltzmann equation provided that $\lambda_{s}$ is large enough
for the non-relativistic approximation to be justified.151515By comparing eq.
(A.3) to eq. (A.5) one can see that the non-relativistic limit is a good
approximation if
$\frac{3}{4}\left[\log\left(\frac{M_{\mathrm{p}}\,\lambda_{s}^{3}}{m_{s}\,\xi^{2/3}}\right)-a\right]\gg
1\;.$ As anticipated, $Y_{s}$ depends only mildly on $\lambda_{s}$, such that
varying $\lambda_{s}$ by orders of magnitude typically only changes $Y_{s}$ by
a factor of a few. Finally, note that the freeze-in yield enters only via the
entropy ratio $\xi$. Since larger freeze-in yields correspond to smaller
$\xi$, we find that increasing the freeze-in yield also increases the final DM
abundance after dark sector freeze-out.

## References

  * [1] L. J. Hall, K. Jedamzik, J. March-Russell, and S. M. West, Freeze-In Production of FIMP Dark Matter,  JHEP 03 (2010) 080, [0911.1120].
  * [2] N. Bernal, M. Heikinheimo, T. Tenkanen, K. Tuominen, and V. Vaskonen, The Dawn of FIMP Dark Matter: A Review of Models and Constraints,  Int. J. Mod. Phys. A32 (2017), no. 27 1730023, [1706.07442].
  * [3] G. Bélanger, F. Boudjema, A. Goudelis, A. Pukhov, and B. Zaldivar, micrOMEGAs5.0 : freeze-in, 1801.03509.
  * [4] R. T. Co, F. D’Eramo, L. J. Hall, and D. Pappadopulo, Freeze-In Dark Matter with Displaced Signatures at Colliders,  JCAP 1512 (2015), no. 12 024, [1506.07532].
  * [5] F. Kahlhoefer, On the LHC sensitivity for non-thermalised hidden sectors,  Phys. Lett. B779 (2018) 388–392, [1801.07621].
  * [6] R. Essig, J. Mardon, and T. Volansky, Direct Detection of Sub-GeV Dark Matter,  Phys. Rev. D85 (2012) 076007, [1108.5383].
  * [7] X. Chu, T. Hambye, and M. H. G. Tytgat, The Four Basic Ways of Creating Dark Matter Through a Portal,  JCAP 1205 (2012) 034, [1112.0493].
  * [8] T. Hambye, M. H. G. Tytgat, J. Vandecasteele, and L. Vanderheyden, Direct Detection is testing Freeze-in, 1807.05022.
  * [9] M. Heikinheimo, T. Tenkanen, and K. Tuominen, Prospects for indirect detection of frozen-in dark matter,  Phys. Rev. D97 (2018), no. 6 063002, [1801.03089].
  * [10] K. S. Babu and R. N. Mohapatra, 7 keV Scalar Dark Matter and the Anomalous Galactic X-ray Spectrum,  Phys. Rev. D89 (2014) 115011, [1404.2220].
  * [11] S. Dodelson and L. M. Widrow, Sterile-neutrinos as dark matter,  Phys. Rev. Lett. 72 (1994) 17–20, [hep-ph/9303287].
  * [12] S. K. Kang and A. Patra, keV Sterile Neutrino Dark Matter and Low Scale Leptogenesis,  J. Korean Phys. Soc. 69 (2016), no. 8 1375–1382, [1412.4899].
  * [13] B. Shakya, Sterile Neutrino Dark Matter from Freeze-In,  Mod. Phys. Lett. A31 (2016), no. 06 1630005, [1512.02751].
  * [14] J. König, A. Merle, and M. Totzauer, keV Sterile Neutrino Dark Matter from Singlet Scalar Decays: The Most General Case,  JCAP 1611 (2016), no. 11 038, [1609.01289].
  * [15] J. Berger, K. Jedamzik, and D. G. E. Walker, Cosmological Constraints on Decoupled Dark Photons and Dark Higgs,  JCAP 1611 (2016) 032, [1605.07195].
  * [16] C. E. Yaguna, The Singlet Scalar as FIMP Dark Matter,  JHEP 08 (2011) 060, [1105.1654].
  * [17] M. Blennow, E. Fernandez-Martinez, and B. Zaldivar, Freeze-in through portals,  JCAP 1401 (2014) 003, [1309.7348].
  * [18] M. Heikinheimo, T. Tenkanen, K. Tuominen, and V. Vaskonen, Observational Constraints on Decoupled Hidden Sectors,  Phys. Rev. D94 (2016), no. 6 063506, [1604.02401]. [Erratum: Phys. Rev.D96,no.10,109902(2017)].
  * [19] E. D. Carlson, M. E. Machacek, and L. J. Hall, Self-interacting dark matter,  Astrophys. J. 398 (1992) 43–52.
  * [20] N. Bernal, C. Garcia-Cely, and R. Rosenfeld, WIMP and SIMP Dark Matter from the Spontaneous Breaking of a Global Group,  JCAP 1504 (2015), no. 04 012, [1501.01973].
  * [21] D. Pappadopulo, J. T. Ruderman, and G. Trevisan, Dark matter freeze-out in a nonrelativistic sector,  Phys. Rev. D94 (2016), no. 3 035005, [1602.04219].
  * [22] M. Farina, D. Pappadopulo, J. T. Ruderman, and G. Trevisan, Phases of Cannibal Dark Matter,  JHEP 12 (2016) 039, [1607.03108].
  * [23] M. Heikinheimo, T. Tenkanen, and K. Tuominen, WIMP miracle of the second kind,  Phys. Rev. D96 (2017), no. 2 023001, [1704.05359].
  * [24] N. Bernal, X. Chu, C. Garcia-Cely, T. Hambye, and B. Zaldivar, Production Regimes for Self-Interacting Dark Matter,  JCAP 1603 (2016), no. 03 018, [1510.08063].
  * [25] Hitomi, F. A. Aharonian et al., $Hitomi$ constraints on the 3.5 keV line in the Perseus galaxy cluster,  Astrophys. J. 837 (2017), no. 1 L15, [1607.07420].
  * [26] E. Bulbul, M. Markevitch, A. Foster, R. K. Smith, M. Loewenstein, et al., Detection of An Unidentified Emission Line in the Stacked X-ray spectrum of Galaxy Clusters,  Astrophys. J. 789 (2014) 13, [1402.2301].
  * [27] A. Boyarsky, O. Ruchayskiy, D. Iakubovskyi, and J. Franse, Unidentified Line in X-Ray Spectra of the Andromeda Galaxy and Perseus Galaxy Cluster,  Phys. Rev. Lett. 113 (2014) 251301, [1402.4119].
  * [28] V. Brdar, J. Kopp, J. Liu, and X.-P. Wang, X-Ray Lines from Dark Matter Annihilation at the keV Scale,  Phys. Rev. Lett. 120 (2018), no. 6 061301, [1710.02146].
  * [29] J. Heeck and D. Teresi, Cold keV dark matter from decays and scatterings,  Phys. Rev. D96 (2017), no. 3 035018, [1706.09909].
  * [30] J. A. Evans, S. Gori, and J. Shelton, Looking for the WIMP Next Door,  JHEP 02 (2018) 100, [1712.03974].
  * [31] M. J. Baker and J. Kopp, Dark Matter Decay between Phase Transitions at the Weak Scale,  Phys. Rev. Lett. 119 (2017), no. 6 061801, [1608.07578].
  * [32] M. J. Baker, M. Breitbach, J. Kopp, and L. Mittnacht, Dynamic Freeze-In: Impact of Thermal Masses and Cosmological Phase Transitions on Dark Matter Production,  JHEP 03 (2018) 114, [1712.03962].
  * [33] M. Quiros, Finite temperature field theory and phase transitions, in Proceedings, Summer School in High-energy physics and cosmology: Trieste, Italy, June 29-July 17, 1998, pp. 187–259, 1999.  hep-ph/9901312.
  * [34] P. Gondolo and G. Gelmini, Cosmic abundances of stable particles: Improved analysis,  Nuclear Physics B 360 (1991), no. 1 145 – 179.
  * [35] J. Redondo and M. Postma, Massive hidden photons as lukewarm dark matter,  JCAP 0902 (2009) 005, [0811.0326].
  * [36] Planck, N. Aghanim et al., Planck 2018 results. VI. Cosmological parameters, 1807.06209.
  * [37] V. Vaskonen, Electroweak baryogenesis and gravitational waves from a real scalar singlet,  Phys. Rev. D95 (2017), no. 12 123515, [1611.02073].
  * [38] A. Belyaev, N. D. Christensen, and A. Pukhov, CalcHEP 3.4 for collider physics within and beyond the Standard Model,  Comput. Phys. Commun. 184 (2013) 1729–1769, [1207.6082].
  * [39] T. Flacke, C. Frugiuele, E. Fuchs, R. S. Gupta, and G. Perez, Phenomenology of relaxion-Higgs mixing,  JHEP 06 (2017) 050, [1610.02025].
  * [40] J. F. Gunion, H. E. Haber, G. L. Kane, and S. Dawson, The Higgs Hunter’s Guide, Front. Phys. 80 (2000) 1–404.
  * [41] A. Djouadi, The Anatomy of electro-weak symmetry breaking. I: The Higgs boson in the standard model,  Phys. Rept. 457 (2008) 1–216, [hep-ph/0503172].
  * [42] M. J. Dolan, F. Kahlhoefer, C. McCabe, and K. Schmidt-Hoberg, A taste of dark matter: Flavour constraints on pseudoscalar mediators,  JHEP 03 (2015) 171, [1412.5174]. [Erratum: JHEP07,103(2015)].
  * [43] V. Poulin, J. Lesgourgues, and P. D. Serpico, Cosmological constraints on exotic injection of electromagnetic energy,  JCAP 1703 (2017), no. 03 043, [1610.10051].
  * [44] P. Stöcker, M. Krämer, J. Lesgourgues, and V. Poulin, Exotic energy injection with ExoCLASS: Application to the Higgs portal model and evaporating black holes,  JCAP 1803 (2018), no. 03 018, [1801.01871].
  * [45] J. D. Bowman, A. E. E. Rogers, R. A. Monsalve, T. J. Mozdzen, and N. Mahesh, An absorption profile centred at 78 megahertz in the sky-averaged spectrum,  Nature 555 (2018), no. 7694 67–70.
  * [46] S. Clark, B. Dutta, Y. Gao, Y.-Z. Ma, and L. E. Strigari, 21 cm limits on decaying dark matter and primordial black holes,  Phys. Rev. D98 (2018), no. 4 043006, [1803.09390].
  * [47] M. Drewes et al., A White Paper on keV Sterile Neutrino Dark Matter,  JCAP 1701 (2017), no. 01 025, [1602.04816].
  * [48] H. Yuksel, J. F. Beacom, and C. R. Watson, Strong Upper Limits on Sterile Neutrino Warm Dark Matter,  Phys. Rev. Lett. 101 (2008) 121301, [0706.4084].
  * [49] K. Perez, K. C. Y. Ng, J. F. Beacom, C. Hersh, S. Horiuchi, et al., Almost closing the νMSM sterile neutrino dark matter window with NuSTAR,  Phys. Rev. D95 (2017), no. 12 123002, [1609.00667].
  * [50] S. Horiuchi, P. J. Humphrey, J. Onorbe, K. N. Abazajian, M. Kaplinghat, et al., Sterile neutrino dark matter bounds from galaxies of the Local Group,  Phys. Rev. D89 (2014), no. 2 025017, [1311.0282].
  * [51] M. Markevitch, A. H. Gonzalez, D. Clowe, A. Vikhlinin, L. David, et al., Direct constraints on the dark matter self-interaction cross-section from the merging galaxy cluster 1E0657-56,  Astrophys. J. 606 (2004) 819–824, [astro-ph/0309303].
  * [52] S. W. Randall, M. Markevitch, D. Clowe, A. H. Gonzalez, and M. Bradac, Constraints on the Self-Interaction Cross-Section of Dark Matter from Numerical Simulations of the Merging Galaxy Cluster 1E 0657-56,  Astrophys. J. 679 (2008) 1173–1180, [0704.0261].
  * [53] F. Kahlhoefer, K. Schmidt-Hoberg, M. T. Frandsen, and S. Sarkar, Colliding clusters and dark matter self-interactions,  Mon. Not. Roy. Astron. Soc. 437 (2014), no. 3 2865–2881, [1308.3419].
  * [54] R. Murgia, A. Merle, M. Viel, M. Totzauer, and A. Schneider, ”Non-cold” dark matter at small scales: a general approach,  JCAP 1711 (2017) 046, [1704.07838].
  * [55] T. Kobayashi, R. Murgia, A. De Simone, V. Iršič, and M. Viel, Lyman-$\alpha$ constraints on ultralight scalar dark matter: Implications for the early and late universe,  Phys. Rev. D96 (2017), no. 12 123514, [1708.00015].
  * [56] R. S. L. Hansen and S. Vogl, Thermalizing sterile neutrino dark matter,  Phys. Rev. Lett. 119 (2017), no. 25 251305, [1706.02707].
  * [57] X. Chu, B. Dasgupta, and J. Kopp, Sterile neutrinos with secret interactions—lasting friendship with cosmology,  JCAP 1510 (2015), no. 10 011, [1505.02795].
  * [58] A. Das, B. Dasgupta and R. Khatri, Ballistic Dark Matter oscillates above $\Lambda$CDM, 1811.00028.
  * [59] K. N. Abazajian, Sterile neutrinos in cosmology,  Phys. Rept. 711-712 (2017) 1–28, [1705.01837].
  * [60] M. Kaplinghat, S. Tulin, and H.-B. Yu, Dark Matter Halos as Particle Colliders: Unified Solution to Small-Scale Structure Puzzles from Dwarfs to Clusters,  Phys. Rev. Lett. 116 (2016), no. 4 041302, [1508.03339].
  * [61] S. Tulin and H.-B. Yu, Dark Matter Self-interactions and Small Scale Structure,  Phys. Rept. 730 (2018) 1–57, [1705.02358].
  * [62] M. Galassi, J. Davies, J. Theiler, B. Cough, G. Jungman, et al., GNU Scientific Library Reference Manual for GSL Version 1.15, 3rd ed., 2011.
  * [63] M. Backovic, K. Kong, and M. McCaskey, MadDM v.1.0: Computation of Dark Matter Relic Abundance Using MadGraph5,  Physics of the Dark Universe 5-6 (2014) 18–28, [1308.4955].