# Structure Formation Limits on Axion-Like Dark Matter

Sven Baumholzer  Vedran Brdar  Enrico Morgante

###### Abstract

We derive structure formation limits on dark matter (DM) composed of keV-scale
axion-like particles (ALPs), produced via freeze-in through the interactions
with photons and Standard Model (SM) fermions. We employ Lyman-alpha
(Ly-$\alpha$) forest data sets as well as the observed number of Milky Way
(MW) subhalos. We compare results obtained using Maxwell-Boltzmann and quantum
statistics for describing the SM bath. It should be emphasized that the
presence of logarithmic divergences complicates the calculation of the
production rate, which can not be parameterized with a simple power law
behaviour. The obtained results, in combination with X-ray bounds, exclude the
possibility for a photophilic “frozen-in” ALP DM with mass below $\sim
19\,\mathrm{keV}$. For the photophobic ALP scenario, in which DM couples
primarily to SM fermions, the ALP DM distribution function is peaked at
somewhat lower momentum and hence for such realization we find weaker limits
on DM mass. Future facilities, such as the upcoming Vera C. Rubin observatory,
will provide measurements with which the current bounds can be significantly
improved to $\sim 80\,\mathrm{keV}$.

MITP-20-073

nuhep-th/20-14

##  1 Introduction

The waning of the WIMP [1] is redirecting (astro)particle DM research towards
alternative DM realizations. For instance, non-thermally produced light DM in
form of sterile neutrinos [2, 3, 4, 5, 6], fuzzy DM [7, 8], hidden photons [9,
10] and ALPs [11, 12, 13] is not only receiving a significant attention in
light of existing and forthcoming terrestrial experiments [14], but also has
the potential to resolve discrepancies between DM observations and simulations
[15]. Additionally, non-thermal DM candidates can be related to new solutions
of the old particle physics problems, such as the hierarchy problem in the
case of relaxion DM [16, 17, 18].

In this work, we focus on ALP DM with $\mathcal{O}(1-100)\,\mathrm{keV}$ mass,
which is the scale motivated by hints in X-ray data [19, 20] as well as the
recent measurement of an excess in electron recoil spectrum performed by the
XENON1T collaboration [21]. In fact, the former measurement can be explained
by ALP coupling to photons [12], while for the latter, ALP coupling to
fermions (namely electrons) suffices [17]. In light of these statements, we
study ALP production via freeze-in through feeble interactions with photons
and SM fermions. Our primary mission is to compute the structure formation
limits that have not been derived to date for keV-scale ALPs; for this purpose
we employ recent Ly-$\alpha$ data as well as the observed number of MW
subhalos. Such structure formation limits are widely scrutinized for light
sterile neutrino DM [22, 23]. While using the results from these studies would
allow one to get a very rough estimate on the structure formation limits for
ALP DM, we find it valuable, especially in light of aforementioned experiments
that may start to detect DM, to perform a dedicated study and determine with
large precision a viable parameter space for the considered model. We also
point out that for bosonic DM such bounds are important irrespective of their
strength; in contrast, structure formation limits on fermionic DM, namely
sterile neutrinos, could be overridden by those stemming from phase space,
_i.e._ Pauli blocking [24].

The paper is organized as follows. In Section 2 we introduce the model and
discuss channels for ALP DM production. We then detail computations for the
approaches in which Maxwell-Boltzmann (2.1) and quantum (2.2) statistics are
employed to describe particles in the SM bath. Then, in Section 3 we discuss
experimental observations that allow us to constrain the ALP DM parameter
space. In Section 4 we discuss theoretical aspects of ALP DM in the early
Universe, focusing on the interplay between misalignment and freeze-in
production. In Section 5 our main results are presented. Finally, we summarize
in Section 6.

##  2 The Model and DM production.

The relevant part of the Lagrangian reads

$\displaystyle\mathcal{L}=\,$
$\displaystyle\frac{1}{2}\partial_{\mu}a\,\partial^{\mu}a+\frac{1}{2}m_{a}^{2}a^{2}+\bar{f}(i\not{\partial}-m_{f})f-\frac{1}{4}F_{\mu\nu}F^{\mu\nu}$
$\displaystyle-q\,e\,\bar{f}\not{A}\,f+\frac{c_{a\gamma\gamma}}{4f_{a}}aF_{\mu\nu}\widetilde{F}^{\mu\nu}+\frac{c_{aff}}{f_{a}}\partial_{\mu}a\,\bar{f}\gamma^{\mu}\gamma_{5}f\,.$
(2.1)

Here, we first write down the kinetic and mass term for the pseudoscalar ALP
field $a$; the following three terms describe quantum electrodynamics, namely
a SM fermion $f$ with an electric charge $q$ interacting with photon field
$A^{\mu}$ whose corresponding (dual) field strength tensor is denoted by
$F^{\mu\nu}$ ($\tilde{F}^{\mu\nu}$); the last two terms describe the
interaction of ALPs with $A^{\mu}$ and $f$, respectively. We will chiefly
analyze ALP coupling to photons and fermions separately, _i.e._ working under
the assumption that either $c_{aff}$ or $c_{a\gamma\gamma}$ vanishes. These
scenarios are dubbed as _photophilic_ and _photophobic_ , respectively.

Regarding the photophilic case, at temperatures above $\sim 160\,\mathrm{GeV}$
the EW symmetry is restored, and one should work in the
$\\{B_{\mu},W^{i}_{\mu}\\}$ basis. We assume for simplicity that ALP couples
only to the $U(1)$ gauge fields through

$\frac{c_{aBB}}{4f_{a}}aB_{\mu\nu}\widetilde{B}^{\mu\nu}\,.$ (2.2)

The inclusion of the $W_{\mu}^{i}$ and gluon fields would increase the
production rate of ALPs, but we chose to limit ourselves to the coupling to
photons, which is the scenario typically considered in the literature for
heavy ALP searches [25, 26, 27, 28].

In the photophobic case, using the equations of motion, one can replace the
last term in Eq. 2.1 with the term proportional to the fermion mass [29],
$(2m_{f}/f_{a})\,a\,\bar{f}\gamma_{5}f$. This does not imply that such
coupling is only operative once electroweak symmetry is broken, for $T\lesssim
160\,\mathrm{GeV}$. On the contrary, after a phase redefinition of the quark
and Higgs fields, a term $iy_{t}(c_{t}a/f_{a})\bar{Q}_{3}Ht_{R}$ is generated
[30], where $Q_{3}$ is the third-generation left-handed doublet, $t_{R}$ is
the right handed top quark, and $c_{t}$ is a free dimensionless parameter.
Such a term, thanks to the large top Yukawa, contributes very efficiently to
the ALP production, and it can be the dominant source when all ALP couplings
to SM particles have the same order of magnitude. In this work, however, we
limit ourselves to the sub-TeV values of the reheating temperature in the
photophobic scenario, namely $T_{\mathrm{RH}}<160\,\mathrm{GeV}$. In that
regime, the top quark is Boltzmann suppressed and does not contribute to ALP
production. Due to the $m_{f}$ suppression of the ALP-fermion coupling, the
relevant fermions for ALP production are bottom and charm quarks as well as
tau leptons.

The kinetic and interaction terms in Eq. 2.1 facilitate ALP production from
thermalized fermions and vector bosons $V$ (the photon, a gluon, or the
$U(1)_{Y}$ gauge field) through $Vf\to af$ and $\bar{f}f\to Va$ channels. For
both of these processes there are one photon- and two $f$-mediated Feynman
diagrams contributing at the tree-level111We have also considered ALP
production from Higgs bosons that couple to $U(1)_{Y}$ gauge field. In the
photophilic case we found such contribution to be negligible, namely at the
percent level with respect to the production through diagrams shown in Fig. 1
(upon summing over all fermion species). Note that in the photophobic
scenario, given the low $T_{\text{RH}}$ under our consideration, the Higgs
boson abundance in the bath is Boltzmann suppressed.(see Fig. 1).

Figure 1: Tree level Feynman diagrams for the two processes of interest,
$fV\to fa$ (top) and $\bar{f}f\to Va$ (bottom), in the photophilic (left) and
photophobic (right) cases.

We note that for intermediate scenarios where $c_{aff}\sim c_{a\gamma\gamma}$,
our results turn out to be driven by couplings to photons, namely they match
the photophilic scenario. This is because, in the photophobic case, the ALP
interaction operator is proportional to the fermion mass, and thus the
production rate is suppressed with respect to the photophilic case by a factor
$m_{f}/T_{\mathrm{RH}}$.

Independently of the scenario under consideration, the production of keV-scale
ALP DM through the processes depicted in Fig. 1 should occur via “freeze-in”
[31], because for $m_{a}\gtrsim 100\,\mathrm{eV}$ the freeze-out would lead to
DM relic abundance that greatly exceeds present measurements [12]. The absence
of ALP DM thermalization simply implies that the reheating temperature,
$T_{\text{RH}}$, should not exceed the temperature at which ALP decouples from
the thermal plasma. The decoupling temperature can be estimated by calculating
the scattering rate, $\Gamma$, from the momentum integrated collision term
given in Eq. 2.20 and the number density of SM fermions,
$n_{f}=72\zeta(3)T^{3}/\pi^{2}$. In the photophilic case it yields

$\displaystyle\Gamma=\frac{\int\mathrm{d}^{3}p\,C(x,r)}{n_{f}}=\frac{5g^{\prime
2}T^{3}}{864\pi^{3}\zeta(3)\cos^{2}\theta_{W}}\frac{c_{a\gamma\gamma}^{2}}{f_{a}^{2}}\left(23-24\gamma_{E}-12\log\left(\frac{11g^{\prime
2}}{48}\right)\right)\,,$ (2.3)

where $\theta_{W}$ ($\cos\theta_{W}\sim 0.88$) is the weak mixing angle and
$\gamma_{E}$ is the Euler-Mascheroni constant. For the numerical values of QED
and $U(1)_{Y}$ coupling constants, $e$ and $g^{\prime}$, we use
$\sqrt{4\pi/137}$ and $0.35$, respectively. We neglect the running of these
parameters.

Fixing the coupling $c_{a\gamma\gamma}/f_{a}$ (using the distribution of Eq.
2.17 and the expression for $\Gamma$ in above equation) in order to match the
observed relic abundance of $\Omega_{\text{DM}}h^{2}\approx 0.12$ we obtain

$\displaystyle\frac{T_{\text{RH}}}{T_{\text{dec}}}\simeq 2.3\times
10^{-4}\left(\frac{g_{*}(T)}{106.75}\right)\left(\frac{10\,\mathrm{keV}}{m_{a}}\right).$
(2.4)

Hence, for $m_{a}\sim\mathcal{O}$(keV) it is guaranteed that the ALP DM never
thermalizes.

For the photophobic scenario we note that $\Gamma/H$ is largest for
temperatures $T\simeq m_{f}$. However, as we will show, demanding that
$\Omega_{\text{DM}}h^{2}$ is matched, one requires $c_{aff}/f_{a}\lesssim
10^{-9}\,\,\mathrm{GeV}^{-1}$ for $m_{a}=\mathcal{O}(10\,\mathrm{keV})$ and
with such feeble couplings it is guaranteed that photophobic ALP DM does not
thermalize with SM species.

Using the distribution functions presented later in the text, for the
contribution from $\bar{f}f\to Ba$ in photophilic scenario we find

$\displaystyle\Omega_{\text{DM}}^{\bar{f}f\to Ba}h^{2}\approx
0.12\left(\frac{106.75}{g_{*}}\right)^{3/2}\left(\frac{c_{a\gamma\gamma}/f_{a}}{10^{-17}\,\mathrm{GeV}^{-1}}\right)^{2}\left(\frac{m_{\text{DM}}}{10\,\text{keV}}\right)\left(\frac{T_{\text{RH}}}{6.7\cdot
10^{16}\,\text{GeV}}\right)\,,$ (2.5)

and the $fB\to fa$ process yields

$\displaystyle\Omega_{\text{DM}}^{fB\to
fa}h^{2}\approx\,0.12\left(\frac{106.75}{g_{*}}\right)^{3/2}\left(\frac{c_{a\gamma\gamma}/f_{a}}{10^{-17}\,\mathrm{GeV}^{-1}}\right)^{2}\left(\frac{m_{\text{DM}}}{10\,\text{keV}}\right)\left(\frac{T_{\text{RH}}}{2.7\cdot
10^{15}\,\text{GeV}}\,\right)\,.$ (2.6)

Including both contributions we get

$\displaystyle\Omega_{\text{DM}}h^{2}\approx\,0.12\left(\frac{106.75}{g_{*}}\right)^{3/2}\left(\frac{c_{a\gamma\gamma}/f_{a}}{10^{-17}\,\mathrm{GeV}^{-1}}\right)^{2}\left(\frac{m_{\text{DM}}}{10\,\text{keV}}\right)\left(\frac{T_{\text{RH}}}{2.6\cdot
10^{15}\,\text{GeV}}\,\right)\,.$ (2.7)

Clearly, the scattering process $fB\to fa$ represents the dominant production
channel in the photophilic scenario and hence will also give the dominant
impact in the determination of structure formation limits. The reason for that
lies in the logarithmic enhancement arising due to regularization of t-channel
$B$-mediated process with a thermal gauge boson mass.

The equivalent equations to Eqs. 2.5 and 2.6 for the photophobic scenario read

$\Omega_{\text{DM}}^{\bar{f}f\to Va}h^{2}\approx
0.12\left(\frac{80}{g_{*}}\right)^{3/2}\left(\frac{m_{\text{DM}}}{10\,\text{keV}}\right)\left(\frac{\sum_{f}\kappa_{f}}{38.9\,\mathrm{GeV}}\right)\left(\frac{c_{aff}/f_{a}}{9.6\cdot
10^{-11}\,\mathrm{GeV}^{-1}}\right)^{2}\,,$ (2.8)

and

$\Omega_{\text{DM}}^{fV\to fa}h^{2}\approx
0.12\left(\frac{80}{g_{*}}\right)^{3/2}\left(\frac{m_{\text{DM}}}{10\,\text{keV}}\right)\left(\frac{\sum_{f}\kappa_{f}}{38.9\,\mathrm{GeV}}\right)\left(\frac{c_{aff}/f_{a}}{1.2\cdot
10^{-10}\,\mathrm{GeV}^{-1}}\right)^{2}\,,$ (2.9)

while the sum gives

$\Omega_{\text{DM}}h^{2}\approx
0.12\left(\frac{80}{g_{*}}\right)^{3/2}\left(\frac{m_{\text{DM}}}{10\,\text{keV}}\right)\left(\frac{\sum_{f}\kappa_{f}}{38.9\,\mathrm{GeV}}\right)\left(\frac{c_{aff}/f_{a}}{7.6\cdot
10^{-11}\,\mathrm{GeV}^{-1}}\right)^{2}\,.$ (2.10)

The factor $\kappa_{f}$ is equal to $m_{f}\,n_{c}\,q^{2}e^{2}$ for leptons and
to $m_{f}\,n_{c}\,q^{2}e^{2}+4\,m_{f}\,g_{s}^{2}$ for quarks, where $n_{c}$ is
the number of colors of the fermion $f$ and $q$ its electric charge. For
quarks, we took into account the diagrams with an external gluon. In this
case, $g_{s}$ is the strong gauge coupling which we fixed to $1.31$ (value
obtained by averaging $\alpha_{s}$ over the relevant energy range between the
b-quark and $Z$ mass), and the factor $4$ comes from the trace over the
$SU(3)$ generators. The parameter $\kappa_{f}$ takes the approximate values
$\\{0.16,8.97,29.0,38.1,38.9\\}\,\mathrm{GeV}$ for tau, charm, bottom, the sum
over these three, and the sum over all SM fermions (excluding the top),
respectively.

One can infer from Eqs. 2.8 and 2.9 that both $\bar{f}f\to Va$ and $fV\to fa$
are of similar strength in the photophobic scenario.

The key ingredient to address structure formation is the ALP DM distribution
function. It is obtained by solving the Boltzmann equation, which, for the
process $1+2\to 3+a$, reads

$\displaystyle\left[\frac{\partial}{\partial t}-Hp_{a}\frac{\partial}{\partial
p_{a}}\right]f(p_{a},t)=\mathcal{C}(p_{a})\,,$ (2.11)

with

$\displaystyle\mathcal{C}(p_{a})=\frac{1}{2E}\int$
$\displaystyle\frac{d^{3}p_{1}}{(2\pi)^{3}2E_{1}}\frac{d^{3}p_{2}}{(2\pi)^{3}2E_{2}}\frac{d^{3}p_{3}}{(2\pi)^{3}2E_{3}}(2\pi)^{4}\delta^{4}(P_{1}+P_{2}-P_{3}-P_{a})\times$
$\displaystyle\times|\mathcal{M}|^{2}f_{1}(E_{1},T)f_{2}(E_{2},T)(1-f_{3}(E_{3},T))\,.$
(2.12)

Here, $|\mathcal{M}|^{2}$ is the squared amplitude for the considered process.
In what follows, we present the computation of the ALP DM distribution
function both for the case of Maxwell-Boltzmann and quantum statistics
describing the SM bath.

###  2.1 ALP distribution function computed using Maxwell-Boltzmann statistics

Using the Maxwell-Boltzmann statistics and assuming $f_{i}\ll 1$ for all
particles involved, simplifies the Boltzmann equation considerably [32].
Indeed, one can approximate
$f_{1}(E_{1})f_{2}(E_{2})(1-f_{3}(E_{3}))\approx\exp[(E_{1}+E_{2})/T]=\exp(P_{0}/T)$,
where $P=P_{1}+P_{2}=P_{3}+P_{4}$ and $P_{0}$ is the first component (energy)
of the 4-vector $P$. The collision term can then be factorized as

$\displaystyle\mathcal{C}(p_{a})$
$\displaystyle=\frac{1}{2E}\int\frac{d^{4}P}{(2\pi)^{4}}\frac{e^{-P_{0}/T}}{2E_{3}}(2\pi)\,\delta(E_{3}+E_{a}-P_{0})\times$
$\displaystyle\times\int\frac{d^{3}p_{2}}{(2\pi)^{3}2E_{2}}\frac{d^{3}p_{3}}{(2\pi)^{3}2E_{3}}(2\pi)^{4}\,\delta^{4}(P_{1}+P_{2}-P)\,|\mathcal{M}|^{2}\,.$
(2.13)

Neglecting any CP violation, the second line is nothing but the reduced cross
section, $\hat{\sigma}$, of the inverse process $3+a\to 1+2$, multiplied by a
phase space factor $\lambda(s,m_{1}^{2},m_{2}^{2})^{1/2}/s$. It is related to
the usual cross section through
$\hat{\sigma}=2\,[\lambda(s,m_{1}^{2},m_{2}^{2})/s]\,\sigma$, where $s$ is a
Mandelstam variable and $\lambda$ is the Källén function. This function
approximately reads $\lambda\approx s^{2}$ when masses
$m_{a},m_{f}\lesssim\sqrt{s}\sim T$ which is the case for the relevant epoch
in the early Universe associated to ALP DM production. The second line of Eq.
2.1 is invariant under longitudinal boosts, which makes it simple to compute
it in the center of momentum frame where $\vec{p}=0$. We computed
$\hat{\sigma}$ analytically at tree level considering the diagrams shown in
Fig. 1. For QED, with one fermion of charge $\pm 1$, we obtained

$\displaystyle\hat{\sigma}_{\bar{f}f\to\gamma
a}=\frac{e^{2}}{12\pi}\frac{c_{a\gamma\gamma}^{2}}{f_{a}^{2}}s\left(1+2\frac{m_{f}^{2}}{s}\right)\sqrt{1-\frac{4m_{f}^{2}}{s}}$
$\displaystyle-\frac{2}{\pi}e^{2}m_{f}^{2}\left(2\frac{c_{aff}^{2}}{f_{a}^{2}}+\frac{c_{aff}c_{a\gamma\gamma}}{f_{a}^{2}}\right)\log\left(\frac{s-2m_{f}^{2}-s\sqrt{1-4m_{f}^{2}/s}}{2m_{f}^{2}}\right)\,,$
(2.14)

and

$\displaystyle\hat{\sigma}_{f\gamma\to
fa}=\frac{e^{2}}{16\pi}\frac{c_{a\gamma\gamma}^{2}}{f_{a}^{2}}s\left(1-\frac{m_{f}^{2}}{s}\right)^{2}\left[4\log\left(\frac{(s-m_{f}^{2})^{2}}{s\,m_{\gamma}^{2}}\right)-3-2\frac{m_{f}^{2}}{s}+\frac{m_{f}^{4}}{s^{2}}\right]$
$\displaystyle-\frac{e^{2}}{2\pi}\left(2\frac{c_{aff}^{2}}{f_{a}^{2}}+\frac{c_{aff}c_{a\gamma\gamma}}{f_{a}^{2}}\right)m_{f}^{2}\left(1-\frac{m_{f}^{2}}{s}\right)\left(2\log\left(\frac{s}{m_{f}^{2}}\right)-3+4\frac{m_{f}^{2}}{s}-\frac{m_{f}^{4}}{s^{2}}\right)\,.$
(2.15)

Here, $m_{\gamma}\approx eT/3$ is the plasmon mass that serves as a regulator
for the diagram including t-channel photon exchange. The expressions in Eqs.
2.14 and 2.15 are obtained for general cases of both $c_{a\gamma\gamma}$ and
$c_{aff}\neq 0$ and have well defined limits if any of these two couplings
vanishes. We note that Eq. 2.15 holds for both $f$ and $\bar{f}$ scattering
off photons, and thus it must be summed over twice.

The expressions in Eqs. 2.14 and 2.15 must be corrected to account for the
presence of SM fermions. In the photophilic scenario, in which the ALP
population is generated at $T\gg 160\,\mathrm{GeV}$, $e^{2}c_{a\gamma\gamma}$
must be replaced by $g^{\prime 2}c_{aBB}=g^{\prime
2}c_{a\gamma\gamma}/\cos^{2}\theta_{W}$, and the result must be multiplied by
the sum over SM fermions $g_{Y}\equiv\sum_{f}Y^{2}n_{c}=10$. The plasmon mass
in this case is $(11/12)^{1/2}g^{\prime}T$ [30]. In the photophobic scenario,
instead, the electric charge $e^{2}$ should be multiplied by $q^{2}n_{c}$ for
leptons, and replaced with $q^{2}n_{c}e^{2}+4g_{s}^{2}$ for quarks, and
finally summed over the SM fermions involved (typically $b,c,\tau$).

It is convenient to rewrite the Boltzmann equation (2.1) by defining
dimensionless quantities $r=m_{H}/T$ and $x=p/T$, where $m_{H}$ is a reference
mass which we fix to be the Higgs mass; its actual value is irrelevant as
$m_{H}$ cancels in the calculation of physical quantities. With this
redefinition one obtains the master formula in Eq. (11) of [32] for ALP DM
distribution function

$\displaystyle f(x,r)$
$\displaystyle=\frac{M_{0}}{16\pi^{2}\,m_{H}\,x^{2}}\int_{r_{i}}^{r_{f}}\,dr\int_{y^{*}}^{\infty}dy\,\hat{\sigma}\left(\frac{m_{H}^{2}\,y}{r^{2}}\right)\,\text{Exp}\left[-x-\frac{y}{4x}\right]\,,$
(2.16)

where $M_{0}=M_{\mathrm{Pl}}\sqrt{45/(4\pi^{3}g_{*})}$ and
$M_{\mathrm{Pl}}=1/G_{N}\approx 1.22\times 10^{19}\,\mathrm{GeV}$ is the
Planck mass; $r_{i,f}$ are the lower and upper integration boundaries
$m_{H}/T_{\mathrm{RH}}$ and $m_{H}/m_{f}$, respectively;
$y^{*}=4rm_{f}^{2}/m_{H}^{2}$ for the fermion annihilation process and
$y^{*}=rm_{f}^{2}/m_{H}^{2}$ for fermion scattering. Using this approach, the
change of the number of relativistic degrees of freedom, $g_{*}$, is not taken
into account but in the relevant temperature range (that is above QCD phase
transition) such approximation is justified.

In Eq. 2.16 and in the following, we stop at the first non-zero order in
$m_{f}$, i.e. we set $m_{f}=0$ everywhere apart from the cross section where
we keep the factor $m_{f}^{2}$ multiplying the $c_{aff}^{2}$ terms. We also
keep only the first non-zero $m_{f}$-dependent term inside the logarithms (see
Eqs. 2.14 and 2.15). The same holds for the integration limits $r_{f}$ and
$y^{*}$, which can be sent, respectively, to infinity and to zero in the
photophilic ALP case. Finally, we drop the interference terms $\propto
c_{a\gamma\gamma}c_{aff}$, as we are going to consider separately the cases
$c_{a\gamma\gamma}=0$ and $c_{aff}=0$. Note that the interference terms in
Eqs. 2.14 and 2.15 are proportional to the corresponding $c_{aff}^{2}$ terms,
and hence, even if considered, they would not lead to any major effect, in
particular the energy dependence of the cross section would be unaltered.

The distribution function can thus be integrated analytically. We obtain

$f(x)=\frac{g^{\prime
2}}{12\pi^{3}\cos^{2}\theta_{W}}g_{Y}M_{0}T_{\mathrm{RH}}\frac{c_{a\gamma\gamma}^{2}}{f_{a}^{2}}e^{-x}\left\\{1+\frac{3}{2}\left[1-4\gamma_{E}+4\log\left(\frac{48x}{11g^{\prime
2}}\right)\right]\right\\}\,,$ (2.17)

in the photophilic and

$\displaystyle f(x)$
$\displaystyle=\sum_{f}\frac{1}{2\pi^{3}}M_{0}\,\kappa_{f}\,\frac{c_{aff}^{2}}{f_{a}^{2}}\frac{e^{-x}}{x}\bigg{\\{}\left[2\sqrt{\pi
x}\,(1+\log
2)\,\mathrm{erf}\left(\frac{1}{\sqrt{x}}\right)+2\Gamma\left(0,\frac{1}{x}\right)\right]$
$\displaystyle+\left[\sqrt{\pi
x}\,\mathrm{erf}\left(\frac{1}{2\sqrt{x}}\right)+2\Gamma\left(0,\frac{1}{4x}\right)\right]\bigg{\\}}\,,$
(2.18)

in the photophobic scenario. In these expressions, _erf_ is the error and
$\Gamma(a,x)$ is the incomplete Gamma function. In both (2.17) and (2.1), the
first term in the curly brackets comes from the fermion annihilation process
$f\bar{f}\to Va$, while the second arises from scattering $fV\to fa$.

Figure 2: ALP momentum distribution. _Left panel:_ In _blue_ we show the
photophilic case for Maxwell-Boltzmann statistics; _black_ lines represent the
photophilic case with quantum statistics and in _green_ we show the
photophobic case for Maxwell-Boltzmann statistics. For each, the dashed (dot-
dashed) line represents $\bar{f}f\to Va$ ($fV\to fa$) and the solid line is
their sum. The vertical dashed line marks $p=g^{\prime}T$, left of which
thermal corrections should be added [33]. _Right panel:_ Photophilic scenario
for Maxwell-Boltzmann and quantum statistics in QED, namely, with only a
single fermion with unit electric charge included.

The left panel of Figure 2 shows the momentum distribution obtained for the
photophilic (blue) and photophobic (green) ALP DM, in units of
$M_{0}T_{\mathrm{RH}}c_{a}^{2}/f_{a}^{2}$, where $c_{a}$ denotes
$c_{a\gamma\gamma}$ or $c_{aff}$, respectively. The black lines correspond to
the photophilic case with quantum statistics, and will be discussed later (see
Sec. 2.2 below). The dashed lines represent the contribution from fermion
annihilation, the dot-dashed ones correspond to fermion scattering, and solid
lines are the sum. We rescaled the photophobic ALP lines by
$T_{\mathrm{RH}}/\kappa_{f}$, as otherwise, for $c_{a\gamma\gamma}\sim
c_{aff}$ choice, green lines would not be visible on the same panel next to
those corresponding to photophilic ALP. Moreover, since for the photophobic
ALP $f(x)\propto m_{f}\,q^{2}n_{c}$, with such rescaling the computed green
curve is independent of the fermion choice. It can be seen from the plot that,
if both $c_{a\gamma\gamma}$ and $c_{aff}$ are present, the contribution from
$c_{aff}$ is negligible as long as $c_{aff}\lesssim
c_{a\gamma\gamma}T_{\mathrm{RH}}/m_{f}$.

Importantly, at $p/T\simeq 4\times 10^{-2}$ the dot-dashed blue line turns
negative and hence the momentum distribution becomes unphysical for smaller
values of $p/T$. The reason for that is a simplified treatment of the IR
divergence in the t-channel photon exchange diagram. For the purpose of our
computations, we make a sharp cut at the location where the function turns
negative. In the case of a Maxwell-Boltzmann distribution, we inferred this
condition to be

$\displaystyle
p/T=\text{Exp}\,[-1/4+\gamma_{E}]\,\frac{m_{\gamma}^{2}}{4T^{2}}\,.$ (2.19)

The same condition holds also for QED, with appropriately adjusted plasmon
mass $m_{\gamma}$. Such single fermion scenario is shown in the right panel of
Fig. 2 and one can observe the shift of the cutoff toward an order of
magnitude smaller values of $p/T$; this effect only arises due to different
values for the plasmon mass in the full model and QED, respectively. In any
case, for both panels, the cutoff occurs at a value of $p/T$ at least two
orders of magnitude smaller than the expected mean, $\langle p\rangle/T$. If
the cutoff was at $p/T\simeq\mathcal{O}(1)$, the calculation would not be
justified.

In passing, we mention that this effect has, to the best of our knowledge, not
been discussed in the literature before. Typically, for UV freeze-in
realizations where the production occurs through higher dimensional operators,
the amplitude for the process is expressed only as a function of Mandelstam
$s$ variable [34, 35, 36] and in such cases there are clearly no problems with
IR divergences. Here, we have, however, demonstrated that in a specific model,
where t-channel diagram exists, full calculation should be performed and the
impact of IR divergences should be assessed in order to have a control over
consistency of the calculation. We expect the unphysical low energy behaviour
to disappear if a fully thermal treatment is adopted and postpone this to a
future work.

Coming back to the left panel of Fig. 2, one can infer that in the photophilic
case, the annihilation process is subdominant, and as such it is often
neglected in the literature. On the other hand, in the photophobic case it
yields stronger contribution with respect to scattering (around a factor of 2,
see also Eqs. 2.8 and 2.9). For the average momentum we obtained $\langle
p\rangle/T=3.24$ for a photophilic ALP and $\langle p\rangle/T=2.36$ for a
photophobic one. Up to higher order corrections, this result is independent of
$m_{f}$.

###  2.2 ALP distribution function computed using quantum statistics

The procedure described in Sec. 2.1 is only approximate, and it fails to
capture two features that may a priori be important. First, it cannot account
correctly for the phase space regions where $p_{i}\lesssim T$, as in these
regions the classical statistics deviates from the quantum one. Second,
thermal field theory effects are important for $p_{i}\lesssim eT$, and must be
correctly resummed. This is especially true for the process $fB\to fa$, which
has a logarithmic IR divergence when $|\vec{p}_{\gamma}-\vec{p}_{a}|\to 0$. In
that context, we are curious whether the usage of quantum statistics could in
principle cure the occurrence of negative momentum distribution function at
low $p/T$, as discussed in Section 2.1. While we found there that the cutoff
to be imposed is at values of $p/T$ that do not effectively impact $\langle
p/T\rangle$, it would still be welcome to have physical values of momentum
distribution function across all available momenta. In order to properly take
into account aforementioned features, we follow the procedure discussed in
Ref. [37].

In the photophilic scenario, and for the process $fB\to fa$, the collision
term of Eq. (2) (the “hard” process) is evaluated by isolating the momentum
flowing in the $t$-channel propagator $\vec{k}=\vec{p}_{\gamma}-\vec{p}_{a}$
and imposing an IR cutoff in the integration, $|\vec{k}|>k_{\mathrm{cut}}$.
For $k_{\mathrm{cut}}\to 0$ the integral diverges as
$\log(T/k_{\mathrm{cut}})$. This divergence is cured by adding a “soft” term
which is extracted from the ALP self-energy evaluated with a resummed photon
propagator at finite temperature. It diverges as
$\log(k_{\mathrm{cut}}/m_{\gamma})$; therefore the sum is finite:
$\log(T/m_{\gamma})=-\log(\sqrt{11/12}g^{\prime})$. The resulting collision
term cannot be integrated fully analytically. The relevant expressions can be
found in Eq. (29) and Eq. (21) of Ref. [37], respectively for the hard and
soft contribution (where the collision term, $\mathcal{C}$, is obtained by
multiplying the interaction rate, $\Gamma$, by the Bose-Einstein distribution
$1/(\exp(E_{a}/T)-1)$. We notice that these expressions are only accurate for
$p_{a}>g^{\prime}T$ (but they correctly account $p_{1,2,3}<g^{\prime}T$).
Below this value, the ALP-photon vertex should be thermally corrected as well
[33].

The collision term for the annihilation process $f\bar{f}\to\gamma a$ in the
photophilic case can also be integrated with the full quantum statistics. By
isolating the $s$-channel momentum $\vec{q}=\vec{p}_{f}+\vec{p}_{\bar{f}}$ and
following a procedure similar to that of Ref. [37], the integrals over $q$ and
the angles can be computed analytically, while a numerical integration must be
performed over the energies of the incoming particles. We did not include
thermal field theory corrections in this case. We checked that, even with the
full quantum statistics, the fermion annihilation process is subdominant
compared to the scattering one, as we already found for the Maxwell-Boltzmann
case (see again Section 2.1).

The momentum distribution is obtained from the collision term by solving the
Boltzmann equation, which is most conveniently written in terms of already
introduced dimensionless quantities

$\displaystyle\frac{m_{H}^{2}}{M_{0}r}\frac{\partial}{\partial
r}f(x,r)=\mathcal{C}(x,r)\,.$ (2.20)

The result is shown in black in the left panel of Fig. 2. The distribution is
similar in shape compared to the one obtained with Maxwell-Boltzmann
statistics. The overall normalization differs by a factor $\int
d^{3}p_{a}f_{\mathrm{MB}}/\int d^{3}p_{a}f_{\mathrm{QS}}\sim 5$; note that
this is not relevant for the purpose of deriving structure formation limits.
The average momentum is, instead, very similar between two cases; we obtained
$\langle p/T\rangle=3.18$ for quantum statistics that is to be compared with
$3.24$ from Maxwell-Boltzmann (see Section 2.1), yielding a $\approx 2\%$
difference. Due to practically identical results on $\langle p/T\rangle$ that
stem from the two approaches, for the calculation of structure formation
limits in Section 5, we choose to proceed with the classical statistics
description that allows us to work with fully analytical expressions, and
leads to physical results even for lower momenta. We note that, as apparent
from Fig. 2 (black lines), the change of the sign of the distribution function
is not resolved with this procedure, and the dot-dashed lines come with a
cutoff (see also right panel for QED). This is, however, not surprising
because as stated above, the calculation is not valid for $p_{a}<g^{\prime}T$.
Interestingly, when quantum statistics is employed, the cutoff appears at
somewhat larger values of momenta in comparison to Maxwell-Boltzmann case.

In passing, we briefly mention that the calculation of the collision term
becomes more involved for a photophobic ALP. For each process, two Feynman
diagrams have to be summed, in different channels (see Fig. 1). The
interference between $s-$ and $t-$channel and between $t-$ and $u-$channel
makes it impossible to define the exchanged momentum univocally, preventing an
analytical integration of the collision term. For this scenario we content
ourselves with the classical statistics results assuming that, as for the case
of photophilic ALP, the quantum effects would affect very mildly the shape of
the distribution and the average momentum $\langle p\rangle/T$.

##  3 Structure formation probes: Lyman-$\alpha$ and satellite counts

Generically, a DM candidate which features a non vanishing distribution
function for rather large values of $p/T$ can be considered to be warm. As
such, it starts to wash out structures at small scales and this can be
quantified by the suppression of the matter power spectrum at such scales.
The respective DM momentum distribution function $f(x,r)$ is essential in the
computation of the matter power spectrum, $P(k)$ (the value of $k$
characterizes inverse scale): It explicitly enters in the calculation of the
DM energy density fluctuation $\delta\rho$, which is needed to construct
$P(k)$. Hence, computing the matter power spectrum, is a starting point and
the essence for assessing the structure formation limits222A more precise
estimation calls for $N$-body simulations approach which goes beyond the scope
of this work..

Qualitatively, this suppression can be quantified by the average momentum
$\langle p/T\rangle$ of the respective ALP DM momentum distribution function
corresponding to both production channels for which we find

$\displaystyle\langle p/T\rangle_{\bar{f}f\to
Va}=\frac{\int_{0}^{\infty}\mathrm{d}x\,x^{3}f_{\bar{f}f\to\gamma
a}(x,r_{\text{RH}})}{\int_{0}^{\infty}\mathrm{d}x\,x^{2}f_{\bar{f}f\to\gamma
a}(x,r_{\text{RH}})}=3.0\,\,(2.31)\,,$ $\displaystyle\langle p/T\rangle_{fV\to
fa}=\frac{\int_{0}^{\infty}\mathrm{d}x\,x^{3}f_{f\gamma\to
fa}(x,r_{\text{RH}})}{\int_{0}^{\infty}\mathrm{d}x\,x^{2}f_{f\gamma\to
fa}(x,r_{\text{RH}})}\approx 3.24\,\,(2.44)\,,$ $\displaystyle\langle
p/T\rangle=\frac{\int_{0}^{\infty}\mathrm{d}x\,x^{3}f(x,r_{\text{RH}})}{\int_{0}^{\infty}\mathrm{d}x\,x^{2}f(x,r_{\text{RH}})}\approx
3.24\,\,(2.36)\,,$ (3.1)

where $r_{\text{RH}}=m_{H}/T_{\text{RH}}$; numbers in each line correspond to
the photophilic and photophobic scnenario, the latter given in brackets. The
last line refers to the sum of the two process, and it is the one which is
physically relevant.

There are two complementary probes that we employ in our analysis: _Ly-
$\alpha$ forest_ and the _number of MW subhalos_. The former stands for a
number of absorption lines occurring in the spectra of quasars and galaxies at
higher redshift, stemming from the hydrogen in the inter-galactic medium [38,
39]. From such “forests” one can obtain the matter power spectrum in one
dimension, along the line of sight.

In the literature, it is customary to define a transfer function, $T(k)$,

$\displaystyle T(k)\equiv\sqrt{\frac{P(k)}{P(k)_{\Lambda\text{CDM}}}}\,,$
(3.2)

where subscript $\Lambda$CDM in the denominator denotes cold dark matter and
the expression in the numerator is the power spectrum for ALP DM. In the
following we calculate it  numerically by feeding $f(x,r)$ into CLASS [40,
41].  Furthermore, we have to take into account entropy dilution effects due
to a change in $g_{*}(T)$ between the time of production and today. This
effect is quantified via an effective DM temperature $T_{\mathrm{DM}}$ which
is related to the photon temperature $T_{\gamma}$ by

$T_{\mathrm{DM}}=\left(\frac{g_{*}(T_{0})}{g_{*}(T_{\mathrm{prod}})}\right)^{1/3}T_{\gamma}=\left(\frac{3.94}{g_{*}(T_{\mathrm{prod}})}\right)^{1/3}T_{\gamma}.$
(3.3)

As can be seen this leads to a “cooling” of the ALP DM and the strength of
this effect depends on the time of production $T_{\mathrm{prod}}$. Above the
electroweak phase transition we have $g_{*}(T)=106.75$ for the SM.

Then, the computed transfer function is compared to an analytical fit of the
transfer function of a warm thermal relic with mass $m_{\mathrm{TR}}$ and
abundance $\Omega_{\mathrm{TR}}$ [38]:

$\displaystyle T(k)=(1+(\alpha k)^{2\nu})^{-5/\nu}\,,$ (3.4)

where $\nu=1.12$ and

$\displaystyle\alpha=0.049\left(\frac{m_{\text{TR}}}{1\,\text{keV}}\right)^{-1.11}\left(\frac{\Omega_{\text{TR}}}{0.25}\right)^{0.11}\left(\frac{h}{0.7}\right)^{1.22}\,h^{-1}\,\text{Mpc}\,.$
(3.5)

Here, TR denotes a thermal relic which constitutes all of DM and for which
structure formation limits are typically derived. The values for the Hubble
constant and $\Omega_{\text{TR}}$ are taken from [42]. The comparison is done
by utilizing the half mode analyis [43] where we define a scale, $k_{1/2}$, at
which the squared transfer function has dropped to $0.5$, i.e.
$T^{2}(k_{1/2})=0.5$. Then, for a given model and the corresponding transfer
function, we check whether $T(k)^{2}\geq
T_{\text{lim}}(k)^{2},\;\forall\,k\leq k_{1/2}$. If not, the considered
parameter choice can be discarded. $T_{\text{lim}}$ is based on Eq. 3.4 and
uses $m_{\text{TR}}$ adopted from dedicated searches using Ly-$\alpha$ forest
observations. The lower bounds on the mass of thermal relic, $m_{\text{TR}}$,
span the range between $m_{\text{TR}}=1.9$ and $5.58\,\mathrm{keV}$ (see [44,
38, 45, 46, 47]). When showing our limits on ALP DM parameter space, we will
present both conservative and aggressive bounds, corresponding to this mass
range of $m_{\text{TR}}$. In Fig. 3 we show a comparison between Eq. 3.2 and
Eq. 3.4 for several values of ALP DM mass. One can infer that the more
stringent limit ($m_{\text{TR}}=5.58\,\mathrm{keV}$) is essentially excluding
ALP masses smaller than $\lesssim 15\,\mathrm{keV}$, whereas by taking
conservative limit of $m_{\text{TR}}=1.9\,\mathrm{keV}$, only $m_{a}\lesssim
5\,\mathrm{keV}$ are disfavored.

While we make use of the half mode analysis to infer limits on $m_{a}$ for a
given $m_{\text{TR}}$, one can relate these two masses also by equalizing the
free-streaming length, $\lambda_{\text{FS}}$, for both ALP and TR models. In
the thermal relic case, $\lambda_{FS}\sim
0.22\,\text{Mpc}\left(\text{keV}/m_{\text{TR}}\right)^{4/3}$, whereas for the
ALP,
$\lambda_{FS}\sim\text{Mpc}\left(\text{keV}/m_{a}\right)\left(g_{*}(T_{\nu})/g_{*}(T_{\mathrm{prod}})\right)^{1/3}$
[32]. Since we deal with  large production temperatures, $T_{\mathrm{prod}}$,
we have  $g{(T_{\mathrm{prod}})}=106.75$, while $g_{*}(T_{\nu})=10.75$. This
allows one to derive the following relation between the two masses,
$(m_{a})_{\mathrm{min}}\simeq
2.1\,\mathrm{keV}\left(\frac{m_{\text{TR}}}{1\,\mathrm{keV}}\right)^{4/3}$; a
similar result was quoted in [48]. A more appropriate derivation which takes
the different momentum distributions into account is to use the half mode
scale $k_{1/2}$. Equalizing it for both models we found the following
numerical relation

$\displaystyle(m_{a})_{\mathrm{min}}\simeq
2.1\,\mathrm{keV}\left(\frac{m_{\text{TR}}}{1\,\mathrm{keV}}\right)^{1.27}.$
(3.6)

Assuming a different value for $g_{*}(T\mathrm{prod})$ one has to modify this
relation by taking the entropy dilution (see Eq. 3.3) into account. As such,
the previous equation is modified by an additional factor

$\displaystyle(m_{a})_{\mathrm{min}}\simeq
2.1\,\mathrm{keV}\left(\frac{m_{\text{TR}}}{1\,\mathrm{keV}}\right)^{1.27}\left(\frac{106.75}{g_{*}(T_{\mathrm{prod}})}\right)^{1/3}.$
(3.7)

Therefore, smaller values of $g_{*}(T_{\mathrm{prod}})$ have to be compensated
by larger ALP DM masses. We will come back to this when discussing limits for
the photophobic scenario.

Figure 3: Transfer function (see Eq. 3.2) shown for photophilic ALP DM with
different masses compared to thermal relic limits presented as an analytic fit
given in Eq. 3.4.

As far as the number of MW satellites is concerned, ALP DM with momentum
distribution as in Eq. 3.1 would cause the occurrence of smaller number of
subhalos, compared to $\Lambda$CDM. The predicted number of satellites has an
intriguing connection with the type of DM and its production; we are able to
set the limit on $m_{a}$ by requiring that the number of satellites is not
smaller that what was observed. For counting subhalos we follow the approach
outlined in [23, 49]: 11 “classical” satellites are combined with 15 ultra-
faint satellites found by SDSS. Latter number is multiplied by a factor of
$7/2$ because of the limited sky coverage of SDSS; this yields in total
$N_{\text{sub}}=64$. One should note that, in addition to SDSS, several more
ultra-faints satellites or candidates have at least been reported (see for
instance [50]) by other surveys, so we regard $N_{\text{sub}}=64$ as a
conservative estimate for the number of MW companions.

We follow the approach from [51, 52, 23] where the authors present a formula
to derive $N_{\text{sub}}$ for a given matter power spectrum, based on an
extended Press-Schechter approach

$\displaystyle\frac{\mathrm{d}N_{sub}}{\mathrm{d}M_{sub}}=\frac{1}{C}\frac{1}{6\pi^{2}}\frac{M_{MW}}{M_{sub}^{2}}\frac{P(1/R_{sub})}{R_{sub}^{3}\sqrt{2\pi(S_{sub}-S_{MW})}}\,.$
(3.8)

Here, $R_{sub}$ is the radius of a subhalo that has a mass $M_{sub}$; these
two quantities are related through
$M_{sub}=\frac{4\pi}{3}\Omega_{m}\rho_{c}(2.5R_{sub})^{3}$. The power spectrum
$P(k)$ of the warm ALP DM has to be integrated for scales smaller compared to
the respective subhalo or MW scale

$\displaystyle
S_{i}=\frac{1}{2\pi^{2}}\int_{0}^{R_{i}^{-1}}\mathrm{d}k\,k^{2}P(k)\,.$ (3.9)

Decreasing $m_{a}$ leads to smaller variances $S_{i}$ (smaller denominator in
Eq. 3.8) as well as a suppression of $P(k)$ which then features an earlier
drop (smaller numerator in Eq. 3.8); the latter effect is, however, stronger
and therefore less numbers of subhalos $N_{sub}$ are predicted.

We further note here that the limits clearly depend on the MW mass, $M_{MW}$.
The precise value of this quantity is still under investigation and, depending
on the analysis, the reported values range roughly between $1\times
10^{12}\,M_{\odot}<M_{MW}<1.5\times 10^{12}\,M_{\odot}$ (see [53, 54, 55, 56,
57] and references therein), using recent data from the GAIA survey. In the
following, we will refer to the lower mass as aggressive, while the higher
mass bound is dubbed as conservative. We want to stress the following two
subtleties: First, the quoted MW masses are defined with respect to densities
$200$ times larger than the critical density, the so called virial mass of the
MW. As such, the constant $C$ of Eq. 3.8 is given by $C=34$. Second, Eq. 3.8
uses MW masses in units of $M_{\odot}/h$ as an input and one has to take this
additional factor into account.

##  4 Suppressing ALP production from misalignment and topological defects

An important assumption of our discussion is that ALP is predominantly
produced through a freeze-in process. This means that the misalignment
contribution, as well as any ALP population produced from the decay of
topological defects such as cosmic strings, must be suppressed.  As we are
going to show in what follows, this is not trivial for the photophilic ALP, as
the requirement of a small misalignment energy density implies bounds on
$T_{\mathrm{RH}}$ and $f_{a}$ which are in contrast with what is needed in
order to match the observed DM abundance with the freeze-in population. For
this reason, in order to justify the assumption that the misalignment energy
density is suppressed, we will need some specific assumption about the axion
potential and the specific symmetry breaking pattern through which it develops
in the early universe. This statement holds only for the photophilic axion, as
we will detail below.

The energy density from ALP misalignment, assuming a constant ALP mass, is
[11]

$\displaystyle\rho_{a,0}=0.17\,\frac{\mathrm{keV}}{\mathrm{cm}^{3}}\left(\frac{m_{a}}{\mathrm{eV}}\right)^{1/2}\left(\frac{a_{i}}{10^{11}\,\mathrm{GeV}}\right)^{2}\,.$
(4.1)

Requiring that this value does not exceed the observed DM abundance yields

$a_{i}<2.6\times
10^{10}\,\mathrm{GeV}\left(\frac{10\,\mathrm{keV}}{m_{a}}\right)^{1/4}\,,$
(4.2)

where $a_{i}$ is the displacement from the minimum when ALP starts to
oscillate. If the ALP is present during inflation333If the ALP is the pseudo-
Goldstone boson of a spontaneously broken global symmetry, this corresponds to
the scenario where the symmetry is broken before or during inflation , with
$H_{I}<f_{a}$. In any case, this assumption is not crucial for our
discussion., the minimal misalignment is set by quantum fluctuations during
inflation, $a_{i}>\sqrt{N_{e}}\frac{H_{I}}{2\pi}$. Here, $N_{e}$ is the number
of efolds of inflation and $H_{I}$ is the Hubble rate during inflation. This,
combined with Eq. (4.2), results in a rather stringent upper limit for
$H_{I}$. Recalling that the maximal reheating temperature is obtained under
the assumption of instantaneous reheating (i.e. the case in which the entire
energy density of inflation is converted into radiation), we obtain a bound on
the reheating temperature

$\displaystyle T_{\mathrm{RH}}<1.2\times
10^{14}\,\mathrm{GeV}\left(\frac{106.75}{g_{*}}\right)^{1/4}\left(\frac{10\,\mathrm{keV}}{m_{a}}\right)^{1/8}\left(\frac{60}{N_{e}}\right)^{1/4}\,.$
(4.3)

With such values of $T_{\text{RH}}$ in the photophilic scenario, matching the
DM abundance would typically be achievable for large $c_{a\gamma\gamma}/f_{a}$
couplings, disfavored by X-ray constraints (see
$T_{\text{RH}}=10^{14}\,\text{GeV}$ line in Fig. 4).  For the photophobic
scenario, instead, this bound is always satisfied in our parameter space.

Let us now present two alternative ways in which the misalignment contribution
can be suppressed without implying a low reheating temperature as in Eq.
(4.3).

If the Peccei-Quinn symmetry breaks after inflation, the misalignment angle,
neglecting anharmonicities in the potential, is averaged to $\langle
a_{i}^{2}\rangle/f_{a}^{2}=\pi^{2}/3$. After inserting this value into Eq.
(4.1) and imposing $\rho_{a,0}<\rho_{\mathrm{CDM},0}$ we obtain

$\displaystyle f_{a}<1.4\times
10^{10}\,\mathrm{GeV}\left(\frac{10\,\mathrm{keV}}{m_{a}}\right)^{1/4}\,.$
(4.4)

In this case we have no upper bound on the reheating temperature. We only need
to assume that the ALP-SM couplings are suppressed compared to $f_{a}$, i.e.
$c_{a\gamma\gamma},c_{aff}\ll\alpha/(2\pi)$, the RHS being the reference value
of this coupling in typical axion models. Further, one needs to make sure that
the contribution from the decay of cosmic string is negligible. For the
typical example of a pNGB from the breaking of a global $U(1)$ symmetry, such
as the standard QCD axion, the axion density from the decay of topological
defects is similar in magnitude to the one from misalignment. Assuming this is
the case, we expect both contributions to be subdominant when Eq. (4.4) is
satisfied.

The second way of suppressing the misalignment energy density is to have a
very small angle from the start. This is achieved if Peccei-Quinn symmetry is
broken already during inflation and the axion is heavy, $m_{a}>H_{I}$. If the
axion mass is constant, this would imply a too strong upper bound on
$T_{\mathrm{RH}}$. Still, the mass of the axion could have been much larger
during inflation than today (see e.g. [58, 59] and references therein for some
explicit realizations). In particular, for $m_{a}\approx
10^{16}\,\mathrm{GeV}$, the reheating temperature can be as large as
$10^{17}\,\mathrm{GeV}$ during inflation. If, after inflation, the ALP is
light again, thermal fluctuations could in principle increase the energy
density from the misalignment mechanism. This is, however, not the case for
the photophilic scenario under our consideration as ALPs do not thermalize
with the SM plasma.

We again stress that, for the photophobic ALP, the requirement of a vanishing
misalignment energy density is easy to satisfy. If the axion is present during
inflation, Eq. (4.3) applies, which is always satisfied for our choice
$T_{\mathrm{RH}}\lesssim\mathcal{O}(100\mathrm{GeV})$. If, instead, the PQ
symmetry is broken after inflation, $f_{a}$ should satisfy Eq. (4.4).
Comparing this requirement with the preferred value of $c_{aff}/f_{a}$ shown
in Fig. 5, obtained by imposing that the freeze-in abundance matches the
observed DM one, results in $c_{aff}\lesssim\mathcal{O}(1)$.

##  5 Results

We use the distribution functions in Eqs. 2.17 and 2.1 to calculate the matter
power spectrum using CLASS, derive the corresponding transfer function and the
number of MW subhalos and finally compare against results from observations
following the strategy outlined in Section 3. As long as the observed DM
abundance is matched, the result does not depend on the couplings
$c_{a\gamma\gamma}$ and $c_{aff}$; hence, we performed a scan over $m_{a}$ to
derive limits.

The results for the Ly-$\alpha$ forest and the MW subhalo counts for the
photophilic scenario are shown in Fig. 4. For both probes, we show
conservative and aggressive bounds based on the arguments presented in Section
3. The explicit values for the lower limits on the DM mass $m_{a}$ are also
quoted in Table 1.  For this scenario we set $g_{*}(T_{\mathrm{prod}})=106.75$
for all parameter choices because of the high reheating scale involved.

| _cons._ | _agg._
---|---|---
Lyman-$\alpha$ | $4.9\,\mathrm{keV}$ | $19.1\,\mathrm{keV}$
MW subhalo | $10.3\,\mathrm{keV}$ | $17.4\,\mathrm{keV}$
Table 1: Structure formation limits in the _photophilic_ scenario.

For this photophilic scenario we also superimpose limits from X-ray searches
which severely narrow the viable parameter space in the $m_{a}\gtrsim
17\,\mathrm{keV}$ region where structure formation constraints fade away. We
have calculated X-ray limits on ALP DM by utilizing existing ones on keV-scale
sterile neutrino DM [60] as well as the expression for photon fluxes in both
models.

The diagonal lines indicate the parameter space for which reheating
temperature equals $T_{\text{RH}}=10^{14}\,\mathrm{GeV}$,
$10^{16}\,\mathrm{GeV}$ and the reduced Planck mass $M_{\text{Pl}}$,
respectively. The middle one corresponds to the upper bound on the energy
scale during inflation set by Planck through the tensor-to-scalar ratio $r$
[61]. If there are scenarios in which the limit can be avoided, the maximal
allowed reheating temperature would correspond to the reduced Planck mass
(lowest line in the figure). As indicated by the plot, most of the parameter
space for a photophilic ALP which is not constrained by structure formation is
already probed by X-ray searches. More available parameter space can be
achieved by considering additional ALP couplings to $W^{a}$ bosons, which
allow for reheating temperatures smaller by a factor $\approx 1/20$ due to the
larger gauge coupling and to the different trace factors. Namely, ALP coupling
$\sim 5$ times smaller with respect to the axion-photon interaction scenario,
suffices for producing the observed amount of dark matter. Finally, we want to
stress that structure formation limits on the photophilic ALP are exceeding
X-ray bounds for $m_{a}\lesssim$ 2 keV independently of the reheating
temperature.

We should briefly discuss the previously reported unidentified line at $\sim
3.5\,\mathrm{keV}$ in a X-ray data spectrum [19, 20]. While such discovery has
chiefly received explanations in terms of decaying keV-scale sterile neutrino
DM, $7\,\mathrm{keV}$ ALP DM was also discussed [12]. Our aggressive structure
formation limits are, however, clearly disfavoring such an interpretation,
while conservative ones are marginally consistent with it.  We should still
stress that our findings are not in general disfavoring a DM interpretation of
a $3.5\,\mathrm{keV}$ line and hold only for the freeze-in production of ALPs.
For instance, scenarios where ALPs are dominantly produced via the
misalignment mechanism are still viable in this regard.

Figure 4: Structure formation limits on the _photophilic_ ALP parameter space
derived using the half mode analysis technique as well as MW subhalo count.
The former (latter) are denoted as $Ly-\alpha$ (MW). For both, conservative
and aggressive limits are shown (see Section 3 for details). What is also
shown are X-ray limits from INTEGRAL [62], NuSTAR [63] and M31 [64]. The
diagonal dashed lines indicate two specific values of $T_{\text{RH}}$,
calculated by the requirement of producing the observed amount of DM.
Sensitivity projection at $m_{a}\sim 80$ keV from forthcoming Vera Rubin
observatory is shown with red line. Figure 5: Structure formation limits on
the _photophobic_ ALP parameter space derived using the half mode analysis
technique as well as MW subhalo count. The former (latter) are denoted as
$Ly-\alpha$ (MW). For both, conservative and aggressive limits are shown (see
Section 3 for details). Theoretically viable parameter space is shown as a
green diagonal band, spanning cases in which ALP couples only to $\tau$-lepton
(upper edge) and where it interacts with $\tau$-lepton, c- and b-quark (lower
edge). Sensitivity projection at $m_{a}\sim 60$ keV from forthcoming Vera
Rubin observatory is shown with red line.

Next, we apply the same techniques to the photophobic scenario; in that case
we derived Section 2.1 by expanding Eq. 2.16 for small $m_{f}$ assuming that
$T_{\text{RH}}\gg m_{f}$. By doing so, the dependence on the reheating
temperature actually drops out444We implicitly assumed that
$T_{\text{RH}}\lesssim 160\,\mathrm{GeV}$ and this approximation holds up to
$T_{\text{RH}}\geq 10\,\mathrm{GeV}$ or so. For smaller reheating temperatures
one has to keep the full expression and an explicit $T_{\text{RH}}$ dependence
is recovered.. The accessible parameter space is then represented by a linear
relation between $m_{a}$ and $c_{aff}$. The results are shown in Fig. 5 where
the green diagonal band denotes the viable parameter space of the model. The
upper edge of the band corresponds to a $\tau$-lepton contribution only, while
the lower edge shows the sum of $\tau$-lepton, c- and b-quark contributions,
with equal coefficients.

The regions shown in blue are disfavored using Ly-$\alpha$ forest and the MW
subhalo counts. The most aggressive limit (see also Table 2) reads
$m_{a}\approx 15\,\mathrm{keV}$ and we hence predict that only the parameter
space above such values could be realizable in Nature. The derived constraints
are slightly weaker compared to the photophilic case, because ALP DM produced
in the photophobic scenario is “cooler” compared to the other case555A recent
study on updated numbers of MW satellites quotes results in terms of
$m_{\text{TR}}>6.5\,\mathrm{keV}$ [65]. As such we can use the ALP DM transfer
function directly to derive mass bounds; after doing so we found
$m_{a}>23.0\,\,\mathrm{keV}$ for photophilic and $m_{a}>17.5\,\,\mathrm{keV}$
for the photophobic scenario, which extends even beyond _aggressive_ bounds
derived in our own analysis. (see Section 2).  However, in the photophobic
scenario we assumed $T_{\mathrm{RH}}\lesssim 160\,\mathrm{GeV}$ and so one has
to take a change in $g_{*}(T_{\mathrm{prod}})$ into account. The results were
derived with $g_{*}(T_{\mathrm{prod}})=106.75$ and so the limits on the ALP
mass have to be rescaled according to Eq. 3.7. For instance, if
$g_{*}(T_{\mathrm{prod}})=80$, all mass limits are increased by a factor
$\simeq 1.1$ compared to the results shown in Table 2.

| _cons._ | _agg._
---|---|---
Lyman-$\alpha$ | $3.7\,\mathrm{keV}$ | $15.5\,\mathrm{keV}$
MW subhalo | $7.8\,\mathrm{keV}$ | $13.3\,\mathrm{keV}$
Table 2: Structure formation limits in the _photophobic_ scenario.

If the ALP has flavor universal couplings or if the strength of ALP couplings
to fermions is at least comparable across several flavors including electrons,
then the red giants [66, 67] and XENON1T [21] limits, shown on Fig. 5, are
applicable since they constrain the ALP-electron coupling. The latter is
clearly more stringent and relevant for larger span of ALP masses, excluding
$c_{aee}/f_{a}\gtrsim 10^{-10}\,\mathrm{GeV}^{-1}$. Clearly, the XENON1T line
is chiefly below the green band in the region $m_{a}\lesssim 10$ keV,
disfavoring the flavor universal ALP coupling scenario for such parameter
choice.

We have also evaluated constraints from X-ray searches using expressions from
[68]. Namely, if ALP couples to fermions, interaction with photons will be
induced at the quantum level. We have, however, found that such X-ray limits
are only relevant if ALP couples to electrons, because the loop induced
coupling to photons is suppressed by the mass ratio $m_{a}^{2}/m_{f}^{2}$. In
such case, even parameter space above $m_{a}\sim 10$ keV, unexcluded by
XENON1T, would be disfavored in flavor universal scenario. We note, however,
that it was recently shown that by adding more new physics, destructive
interference between loop-level diagrams contributing to decays into SM
photons can be achieved; hence the limit can be relaxed substantially [69],
particularly for such higher values of $m_{a}$.

We note here that our results constrain the relaxion DM model of Ref. [16],
in which relaxion production occurs via freeze-in. The relaxion can constitute
100% of DM only if a large mass is allowed by a large hierarchy in the
relaxion couplings. Even assuming the most conservative bounds from the number
of MW subhalos, the model provides an explanation of the XENON1T anomaly only
if the relaxion is responsible for a fraction of the DM abundance [17].

Finally, let us address future sensitivity projections. Utilizing data from
the upcoming Vera C. Rubin observatory, designed to measure MW mass halo
function down to $10^{6}\,M_{\odot}$, it will be possible to improve limits
dramatically. Using its projected thermal mass limit
$m_{\text{TR}}>18\,\mathrm{keV}$ [70, 71] we find that, for the photophilic
case, the lower limit on the ALP DM mass would be pushed to
$m_{a}>83\,\mathrm{keV}$. Similarly, we find $m_{a}>62\,\mathrm{keV}$ for the
photophobic case (both values are illustrated in Figs. 4 and 5 with red
lines). It is rather intriguing that for the case of photophilic ALP, the
structure formation limits are expected to start competing with the long-
standing X-ray limits in $m_{a}\sim\mathcal{O}(100)$ keV region.

##  6 Summary and Conclusion

ALPs are currently one of the most popular beyond the SM extensions, being
studied and tested across several mass scales. In this paper we studied keV-
scale ALP DM produced via freeze-in through feeble interactions with photons
and SM fermions. The respective momentum distribution has been calculated
assuming Maxwell-Boltzmann and quantum statistics. Although these approaches
do not capture the full picture at small values of $p/T$, they turn out to be
robust for constraining the parameter space using the matter power spectrum
derived from the respective momentum distribution function. Using
Lyman-$\alpha$ forest data as well as the observed number of MW companions we
derived structure formation limits that were missing to date. For the
photophilic ALP DM, we found the most aggressive limits to exclude ALP DM
masses below $\sim 19\,\mathrm{keV}$, complementing constraints from X-ray
data. For photophobic ALP, the obtained limits are somewhat milder, because in
that case we found DM to be “cooler”. Utilizing measurements from the upcoming
Vera Rubin observatory, such bounds will be strongly improved to $m_{a}\gtrsim
60\,\mathrm{keV}$ for both scenarios.

## Acknowledgements

We would like to thank Joerg Jaeckel for discussions. The work of SB and EM is
supported by the Cluster of Excellence “Precision Physics, Fundamental
Interactions, and Structure of Matter” (PRISMA+ EXC 2118/1) funded by the
German Research Foundation (DFG) within the German Excellence Strategy
(Project ID 39083149). Fermilab is operated by the Fermi Research Alliance,
LLC under contract No. DE-AC02-07CH11359 with the United States Department of
Energy.

## References

  * [1] G. Arcadi, M. Dutra, P. Ghosh, M. Lindner, Y. Mambrini, M. Pierre, S. Profumo, and F. S. Queiroz, The waning of the WIMP? A review of models, searches, and constraints, Eur. Phys. J. C 78 (2018), no. 3 203, [arXiv:1703.07364].
  * [2] S. Dodelson and L. M. Widrow, Sterile-neutrinos as dark matter, Phys. Rev. Lett. 72 (1994) 17–20, [hep-ph/9303287].
  * [3] X.-D. Shi and G. M. Fuller, A New dark matter candidate: Nonthermal sterile neutrinos, Phys. Rev. Lett. 82 (1999) 2832–2835, [astro-ph/9810076].
  * [4] A. Merle, V. Niro, and D. Schmidt, New Production Mechanism for keV Sterile Neutrino Dark Matter by Decays of Frozen-In Scalars, JCAP 03 (2014) 028, [arXiv:1306.3996].
  * [5] 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, [arXiv:1710.02146].
  * [6] A. De Gouvêa, M. Sen, W. Tangarife, and Y. Zhang, Dodelson-Widrow Mechanism in the Presence of Self-Interacting Neutrinos, Phys. Rev. Lett. 124 (2020), no. 8 081802, [arXiv:1910.04901].
  * [7] J.-W. Lee, Brief History of Ultra-light Scalar Dark Matter Models, EPJ Web Conf. 168 (2018) 06005, [arXiv:1704.05057].
  * [8] L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten, Ultralight scalars as cosmological dark matter, Phys. Rev. D 95 (2017), no. 4 043541, [arXiv:1610.08297].
  * [9] J. Redondo and M. Postma, Massive hidden photons as lukewarm dark matter, JCAP 02 (2009) 005, [arXiv:0811.0326].
  * [10] T. Hambye, M. H. Tytgat, J. Vandecasteele, and L. Vanderheyden, Dark matter from dark photons: a taxonomy of dark matter production, Phys. Rev. D 100 (2019), no. 9 095018, [arXiv:1908.09864].
  * [11] P. Arias, D. Cadamuro, M. Goodsell, J. Jaeckel, J. Redondo, and A. Ringwald, WISPy Cold Dark Matter, JCAP 06 (2012) 013, [arXiv:1201.5902].
  * [12] J. Jaeckel, J. Redondo, and A. Ringwald, 3.55 keV hint for decaying axionlike particle dark matter, Phys. Rev. D 89 (2014) 103511, [arXiv:1402.7335].
  * [13] S. H. Im and K. S. Jeong, Freeze-in Axion-like Dark Matter, Phys. Lett. B 799 (2019) 135044, [arXiv:1907.07383].
  * [14] J. L. Ouellet et al., First Results from ABRACADABRA-10 cm: A Search for Sub-$\mu$eV Axion Dark Matter, Phys. Rev. Lett. 122 (2019), no. 12 121802, [arXiv:1810.12257].
  * [15] D. J. E. Marsh and J. Silk, A Model For Halo Formation With Axion Mixed Dark Matter, Mon. Not. Roy. Astron. Soc. 437 (2014), no. 3 2652–2663, [arXiv:1307.1705].
  * [16] N. Fonseca and E. Morgante, Relaxion Dark Matter, Phys. Rev. D 100 (2019), no. 5 055010, [arXiv:1809.04534].
  * [17] N. Fonseca and E. Morgante, Probing photophobic (rel)axion dark matter, arXiv:2009.10974.
  * [18] A. Banerjee, H. Kim, and G. Perez, Coherent relaxion dark matter, Phys. Rev. D 100 (2019), no. 11 115026, [arXiv:1810.01889].
  * [19] E. Bulbul, M. Markevitch, A. Foster, R. K. Smith, M. Loewenstein, and S. W. Randall, Detection of An Unidentified Emission Line in the Stacked X-ray spectrum of Galaxy Clusters, Astrophys. J. 789 (2014) 13, [arXiv:1402.2301].
  * [20] 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, [arXiv:1402.4119].
  * [21] XENON Collaboration, E. Aprile et al., Excess electronic recoil events in XENON1T, Phys. Rev. D 102 (2020), no. 7 072004, [arXiv:2006.09721].
  * [22] A. Merle, A. Schneider, and M. Totzauer, Dodelson-Widrow Production of Sterile Neutrino Dark Matter with Non-Trivial Initial Abundance, JCAP 04 (2016) 003, [arXiv:1512.05369].
  * [23] A. Schneider, Astrophysical constraints on resonantly produced sterile neutrino dark matter, JCAP 04 (2016) 059, [arXiv:1601.07553].
  * [24] S. Tremaine and J. Gunn, Dynamical Role of Light Neutral Leptons in Cosmology, Phys. Rev. Lett. 42 (1979) 407–410.
  * [25] J. L. Feng, I. Galon, F. Kling, and S. Trojanowski, Axionlike particles at FASER: The LHC as a photon beam dump, Phys. Rev. D 98 (2018), no. 5 055021, [arXiv:1806.02348].
  * [26] M. J. Dolan, T. Ferber, C. Hearty, F. Kahlhoefer, and K. Schmidt-Hoberg, Revised constraints and Belle II sensitivity for visible and invisible axion-like particles, JHEP 12 (2017) 094, [arXiv:1709.00009].
  * [27] S. Alekhin et al., A facility to Search for Hidden Particles at the CERN SPS: the SHiP physics case, Rept. Prog. Phys. 79 (2016), no. 12 124201, [arXiv:1504.04855].
  * [28] V. Brdar, B. Dutta, W. Jang, D. Kim, I. M. Shoemaker, Z. Tabrizi, A. Thompson, and J. Yu, Axion-like Particles at Future Neutrino Experiments: Closing the ”Cosmological Triangle”, arXiv:2011.07054.
  * [29] M. Bauer, M. Neubert, and A. Thamm, Collider Probes of Axion-Like Particles, JHEP 12 (2017) 044, [arXiv:1708.00443].
  * [30] A. Salvio, A. Strumia, and W. Xue, Thermal axion production, JCAP 01 (2014) 011, [arXiv:1310.6982].
  * [31] L. J. Hall, K. Jedamzik, J. March-Russell, and S. M. West, Freeze-In Production of FIMP Dark Matter, JHEP 03 (2010) 080, [arXiv:0911.1120].
  * [32] J. Heeck and D. Teresi, Cold keV dark matter from decays and scatterings, Phys. Rev. D 96 (2017), no. 3 035018, [arXiv:1706.09909].
  * [33] E. Braaten and T. C. Yuan, Calculation of screening in a hot plasma, Phys. Rev. Lett. 66 (1991) 2183–2186.
  * [34] F. Elahi, C. Kolda, and J. Unwin, UltraViolet Freeze-in, JHEP 03 (2015) 048, [arXiv:1410.6157].
  * [35] G. Ballesteros, M. A. Garcia, and M. Pierre, How warm are non-thermal relics? Lyman-$\alpha$ bounds on out-of-equilibrium dark matter, arXiv:2011.13458.
  * [36] F. D’Eramo and A. Lenoci, Lower Mass Bounds on FIMPs, arXiv:2012.01446.
  * [37] M. Bolz, A. Brandenburg, and W. Buchmuller, Thermal production of gravitinos, Nucl. Phys. B 606 (2001) 518–544, [hep-ph/0012052]. [Erratum: Nucl.Phys.B 790, 336–337 (2008)].
  * [38] M. Viel, J. Lesgourgues, M. G. Haehnelt, S. Matarrese, and A. Riotto, Constraining warm dark matter candidates including sterile neutrinos and light gravitinos with WMAP and the Lyman-alpha forest, Phys. Rev. D 71 (2005) 063534, [astro-ph/0501562].
  * [39] V. K. Narayanan, D. N. Spergel, R. Dave, and C.-P. Ma, Constraints on the mass of warm dark matter particles and the shape of the linear power spectrum from the Ly$\alpha$ forest, Astrophys. J. Lett. 543 (2000) L103–L106, [astro-ph/0005095].
  * [40] J. Lesgourgues, The Cosmic Linear Anisotropy Solving System (CLASS) I: Overview, arXiv:1104.2932.
  * [41] J. Lesgourgues and T. Tram, The Cosmic Linear Anisotropy Solving System (CLASS) IV: efficient implementation of non-cold relics, JCAP 09 (2011) 032, [arXiv:1104.2935].
  * [42] Planck Collaboration, N. Aghanim et al., Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641 (2020) A6, [arXiv:1807.06209].
  * [43] J. König, A. Merle, and M. Totzauer, keV Sterile Neutrino Dark Matter from Singlet Scalar Decays: The Most General Case, JCAP 11 (2016) 038, [arXiv:1609.01289].
  * [44] A. Garzilli, O. Ruchayskiy, A. Magalich, and A. Boyarsky, How warm is too warm? Towards robust Lyman-$\alpha$ forest bounds on warm dark matter, arXiv:1912.09397.
  * [45] C. Yèche, N. Palanque-Delabrouille, J. Baur, and H. du Mas des Bourboux, Constraints on neutrino masses from Lyman-alpha forest power spectrum with BOSS and XQ-100, JCAP 06 (2017) 047, [arXiv:1702.03314].
  * [46] V. Iršič et al., New Constraints on the free-streaming of warm dark matter from intermediate and small scale Lyman-$\alpha$ forest data, Phys. Rev. D 96 (2017), no. 2 023522, [arXiv:1702.01764].
  * [47] J.-W. Hsueh, W. Enzi, S. Vegetti, M. Auger, C. D. Fassnacht, G. Despali, L. V. Koopmans, and J. P. McKean, SHARP – VII. New constraints on the dark matter free-streaming properties and substructure abundance from gravitationally lensed quasars, Mon. Not. Roy. Astron. Soc. 492 (2020), no. 2 3047–3059, [arXiv:1905.04182].
  * [48] K. J. Bae, A. Kamada, S. P. Liew, and K. Yanagi, Light axinos from freeze-in: production processes, phase space distributions, and Ly-$\alpha$ forest constraints, JCAP 01 (2018) 054, [arXiv:1707.06418].
  * [49] R. Murgia, A. Merle, M. Viel, M. Totzauer, and A. Schneider, ”Non-cold” dark matter at small scales: a general approach, JCAP 11 (2017) 046, [arXiv:1704.07838].
  * [50] DES Collaboration, A. Drlica-Wagner et al., Milky Way Satellite Census. I. The Observational Selection Function for Milky Way Satellites in DES Y3 and Pan-STARRS DR1, Astrophys. J. 893 (2020) 1, [arXiv:1912.03302].
  * [51] A. Schneider, R. E. Smith, and D. Reed, Halo Mass Function and the Free Streaming Scale, Mon. Not. Roy. Astron. Soc. 433 (2013) 1573, [arXiv:1303.0839].
  * [52] A. Schneider, Structure formation with suppressed small-scale perturbations, Mon. Not. Roy. Astron. Soc. 451 (2015), no. 3 3117–3130, [arXiv:1412.2133].
  * [53] W. Wang, J. Han, A. P. Cooper, S. Cole, C. Frenk, and B. Lowing, Estimating the dark matter halo mass of our Milky Way using dynamical tracers, Mon. Not. Roy. Astron. Soc. 453 (2015), no. 1 377–400, [arXiv:1502.03477].
  * [54] T. Callingham, M. Cautun, A. J. Deason, C. S. Frenk, W. Wang, F. A. Gómez, R. J. Grand, F. Marinacci, and R. Pakmor, The mass of the Milky Way from satellite dynamics, arXiv:1808.10456.
  * [55] M. Cautun, A. Benitez-Llambay, A. J. Deason, C. S. Frenk, A. Fattahi, F. A. Gómez, R. J. Grand, K. A. Oman, J. F. Navarro, and C. M. Simpson, The Milky Way total mass profile as inferred from Gaia DR2, Mon. Not. Roy. Astron. Soc. 494 (2020), no. 3 4291–4313, [arXiv:1911.04557].
  * [56] E. Karukes, M. Benito, F. Iocco, R. Trotta, and A. Geringer-Sameth, A robust estimate of the Milky Way mass from rotation curve data, JCAP 05 (2020) 033, [arXiv:1912.04296].
  * [57] W. Wang, J. Han, M. Cautun, Z. Li, and M. N. Ishigaki, The mass of our Milky Way, Sci. China Phys. Mech. Astron. 63 (2020), no. 10 109801, [arXiv:1912.02599].
  * [58] G. Dvali, Removing the cosmological bound on the axion scale, hep-ph/9505253.
  * [59] R. T. Co, E. Gonzalez, and K. Harigaya, Axion Misalignment Driven to the Bottom, JHEP 05 (2019) 162, [arXiv:1812.11186].
  * [60] M. Drewes et al., A White Paper on keV Sterile Neutrino Dark Matter, JCAP 01 (2017) 025, [arXiv:1602.04816].
  * [61] Planck Collaboration, Y. Akrami et al., Planck 2018 results. X. Constraints on inflation, Astron. Astrophys. 641 (2020) A10, [arXiv:1807.06211].
  * [62] A. Boyarsky, D. Malyshev, A. Neronov, and O. Ruchayskiy, Constraining DM properties with SPI, Mon. Not. Roy. Astron. Soc. 387 (2008) 1345, [arXiv:0710.4922].
  * [63] B. M. Roach, K. C. Ng, K. Perez, J. F. Beacom, S. Horiuchi, R. Krivonos, and D. R. Wik, NuSTAR Tests of Sterile-Neutrino Dark Matter: New Galactic Bulge Observations and Combined Impact, Phys. Rev. D 101 (2020), no. 10 103011, [arXiv:1908.09037].
  * [64] S. Horiuchi, P. J. Humphrey, J. Onorbe, K. N. Abazajian, M. Kaplinghat, and S. Garrison-Kimmel, Sterile neutrino dark matter bounds from galaxies of the Local Group, Phys. Rev. D 89 (2014), no. 2 025017, [arXiv:1311.0282].
  * [65] DES Collaboration, E. Nadler et al., Milky Way Satellite Census. III. Constraints on Dark Matter Properties from Observations of Milky Way Satellite Galaxies, arXiv:2008.00022.
  * [66] N. Viaux, M. Catelan, P. B. Stetson, G. Raffelt, J. Redondo, A. A. R. Valcarce, and A. Weiss, Neutrino and axion bounds from the globular cluster M5 (NGC 5904), Phys. Rev. Lett. 111 (2013) 231301, [arXiv:1311.1669].
  * [67] M. M. Miller Bertolami, B. E. Melendez, L. G. Althaus, and J. Isern, Revisiting the axion bounds from the Galactic white dwarf luminosity function, JCAP 10 (2014) 069, [arXiv:1406.7712].
  * [68] N. Craig, A. Hook, and S. Kasko, The Photophobic ALP, JHEP 09 (2018) 028, [arXiv:1805.06538].
  * [69] C. Benso, V. Brdar, M. Lindner, and W. Rodejohann, Prospects for Finding Sterile Neutrino Dark Matter at KATRIN, Phys. Rev. D 100 (2019), no. 11 115035, [arXiv:1911.00328].
  * [70] LSST Dark Matter Group Collaboration, A. Drlica-Wagner et al., Probing the Fundamental Nature of Dark Matter with the Large Synoptic Survey Telescope, arXiv:1902.01055.
  * [71] C. Dvorkin, T. Lin, and K. Schutz, The cosmology of sub-MeV dark matter freeze-in, arXiv:2011.08186.