# ALP Dark Matter in a Primordial Black Hole Dominated Universe

Nicolás Bernal  Yuber F. Perez-Gonzalez
Yong Xu  and Óscar Zapata

###### Abstract

We investigate the phenomenological consequences of axion-like particle (ALP)
dark matter with an early matter domination triggered by primordial black
holes (PBHs). We focus on light BHs with masses smaller than $\sim 10^{9}~{}$g
which fully evaporate before Big Bang nucleosynthesis. We numerically solve
the coupled Boltzmann equations, carefully taking the greybody factors and BH
angular momentum into account. We find that the entropy injection from PBH
evaporation dilutes the ALP relic abundance originally produced via the vacuum
misalignment mechanism, opening the parameter space with larger scales $f_{a}$
or, equivalently, smaller ALP-photon couplings $g_{a\gamma}$, within the reach
of future detectors as ABRACADABRA, KLASH, ADMX, and DM-Radio. Moreover, the
ALP minicluster masses can be several orders of magnitude larger if the early
Universe features an PBH dominated epoch. For the relativistic ALPs produced
directly from Hawking radiation, we find that their contribution to the dark
radiation is within the sensitivity of next generation CMB experiments. For
the sake of completeness, we also revisit the particular case of the QCD
axion.

PI/UAN-2021-702FT, FERMILAB-PUB-21-478-T, NUHEP-TH/21-16, IPPP/21/36

##  1 Introduction

Primordial black holes (PBHs), which could have been copiously produced in the
early Universe due to large density fluctuations, are attracting intensive
investigations, for recent reviews see e.g. Refs. [1, 2, 3, 4]. Those with
masses larger than $\sim 10^{15}$ g have not fully evaporated at present and
are potential cold dark matter (DM) candidates. On the contrary, light PBHs
could have also been created, and have not a less interesting phenomenology.
For example, in their evaporation process they can source particle DM or dark
radiation (DR) [5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17], trigger
baryogenesis [18, 19, 20, 21, 22], and radiate all new degrees of freedom,
such as right-handed neutrinos [23]. Along the same research direction, in
this work we focus on phenomenology of QCD axion and axion-like particles
(ALPs) in a Universe dominated by PBHs lighter than $\sim 10^{9}$ g, which
have fully evaporated before Big Bang nucleosynthesis (BBN).

The QCD axion [24, 25], a pseudo-Nambu-Goldstone boson resulting from the
spontaneous breaking of the global $U(1)$ Peccei-Quinn symmetry [26, 27, 28],
is a well-motivated candidate for DM. In the early Universe, it can be
produced via a number of processes, being the vacuum misalignment mechanism
[29, 30, 31] and the decay of topological defects [32] the most popular in the
literature. In the standard cosmological scenario, the axion window is rather
narrow, in particular an axion with mass $m_{a}\simeq[10^{-6}-10^{-5}]$ eV (or
correspondingly to a Peccei-Quinn scale $f_{a}\simeq 10^{12}$ GeV) is
expected, if one does not want to introduce fine tuning of the initial
misalignment angle $\theta_{i}$ [33, 34, 35]. In light of the accumulated
interesting implications for PBHs [36, 37, 38, 39, 40, 41], it has been
recently shown that the axion DM window can be enlarged to a mass as low as
$m_{a}\sim\mathcal{O}(10^{-8})$ eV if the Universe features an early PBHs
dominated epoch [42].

In this paper, we go beyond the QCD axion DM and consider general ALPs, see
e.g. Refs. [43, 44, 45, 33, 46, 34]. ALPs could also arise from the
spontaneous breaking of a global $U(1)$ symmetry, similar to the QCD axion.
They are also quite ubiquitous in string theory [47, 48]. Differently from QCD
axions, ALPs do not solve the strong CP problem since they are in general not
involved in the strong interaction. However, they serve as good candidate for
DM, see e.g. Ref. [45]. In this paper, we focus on ALP DM generated via the
usual misalignment mechanism. Therefore, our results are complementary to
those presented in Ref. [41]. In the standard cosmological scenario, the ALP
relic density depends on three parameters: the decay constant $f_{a}$, its
mass $m_{a}$ and the initial misalignment angle $\theta_{i}$, leading to
strong bounds on the viable parameter space for natural values of
$\theta_{i}$. However, new regions of the parameter space become viable with a
nonstandard cosmological epoch before BBN. We note that the phenomenology of
ALPs as DM with nonstandard cosmological epochs such as kination or an early
matter phase was investigated in Ref. [49] recently, where it was shown that a
broader parameter is available and could be within the reach of several
proposed experiments. Similarly to Ref. [42], and differently from Ref. [49],
we focus on the phenomenological consequences of ALPs in an early matter epoch
triggered by PBHs.

Since PBHs effectively behave as matter, with an energy density that redshifts
slower than radiation, it could be expected that the early Universe has
undergone a PBH dominated epoch. Moreover, it is worth to mention some
differences between an early matter dominated epoch triggered by PBHs and by a
long-lived heavy particle. In the former case, the evaporation (or decay) rate
is time dependent while the later is not. Besides, an early PBH era would
inevitably give rise to gravitational wave (GW) signatures [50, 51, 52, 53],
offering an interesting avenue to constrain on the initial amount of energy
density in the early matter epoch, while such kind of constraints are absent
in case of a heavy particle.

For the case with an early PBH era, we numerically solve a system of coupled
Boltzmann equations for the background (PBHs and SM radiation) based on the
code developed in Refs. [15, 16], carefully including the greybody factor in
PBHs spectra, and considering the effect of the PBH angular momentum. We find
an increased allowed parameter space with larger $f_{a}$ or, equivalently,
smaller ALP-photon coupling $g_{a\gamma}$ due to the entropy injection from
PBH evaporation. Furthermore, in the scenario where the Peccei-Quinn symmetry
is spontaneously broken after inflation, the gravitational clump of ALP
(axion) density inhomogeneities at the time of the matter-radiation equality
gives rise to miniclusters [54, 55, 56, 57, 58]. Due to the PBH domination
epoch, miniclusters with larger masses could be formed.

Additionally, ALPs and axions are directly radiated from PBHs evaporation,
being ultra-relativistic, and thus contributing to the DR. Taking carefully
into account the effect of the PBH spin, we numerically compute the
contribution to $\Delta N_{\text{eff}}$, which is within the sensitivity of
next generation CMB experiments [59, 60, 61, 62]. We also find that for Kerr
PBHs, $\Delta N_{\text{eff}}$ turns out to be smaller than for nonrotating
Schwarzschild PBHs.

The paper is organized as follows. We first revisit PBH evaporation for both
Schwarzschild and Kerr BHs in Sec. 2. In Sec. 3, we estimate the background
evolution by setting up the coupled Boltzmann equations and the formalism for
estimating the entropy injection. Sec. 4 is devoted to the phenomenological
consequences of ALP DM. For the sake of completeness, in Sec. 5 we revisit QCD
axion DM. Additionally, in Sec. 6 we focus on the axion and ALP DR directly
produced from PBH evaporation. Finally we sum up our findings in Sec. 7. We
use natural units where $\hbar=c=k_{\rm B}=1$ throughout this manuscript.

##  2 Primordial Black Hole Evaporation

PBHs are hypothetical objects which could be generated due to inhomogeneities
of density perturbations in the early Universe [63, 64]. When these
fluctuations reenter the horizon, if above a threshold, they could collapse
and form a BH according to the Press-Schechter formalism [65]. In this paper,
we focus on the case where PBHs form in a radiation-dominated epoch. With an
initial cosmic temperature $T=T_{\text{in}}$, the initial PBH mass is given by
the whole mass within the particle horizon [66, 1],

$M_{\text{in}}\equiv
M_{\text{BH}}(T_{\text{in}})=\frac{4\pi}{3}\,\gamma\,\frac{\rho_{R}(T_{\text{in}})}{H^{3}(T_{\text{in}})}\,,$
(2.1)

where we take the efficiency factor to be $\gamma\simeq 0.2$. Additionally,
$\rho_{R}$ and $H$ correspond to the standard model (SM) energy density and
the Hubble expansion rate, respectively. For simplicity, we assume a
monochromatic mass spectrum, such that all PBHs were created with the same
mass. The extension to more realistic mass distributions will be considered
elsewhere. The PBH initial energy density $\rho_{\text{BH}}(T_{\text{in}})$ is
usually related to the SM radiation energy density at formation via the
$\beta$ parameter

$\beta\equiv\frac{\rho_{\text{BH}}(T_{\text{in}})}{\rho_{R}(T_{\text{in}})}\,.$
(2.2)

Since $\rho_{\text{BH}}$ redshifts slower than radiation, an early PBH-
dominated epoch with $\rho_{\text{BH}}>\rho_{R}$ is triggered if
$\beta>\beta_{c}$, with

$\beta_{c}\equiv\frac{T_{\text{ev}}}{T_{\text{in}}}\,,$ (2.3)

where $T_{\text{ev}}$, given by [42]

$T_{\text{ev}}\simeq\left(\frac{g_{\star}(T_{\text{in}})}{640}\right)^{1/4}\left(\frac{M_{P}^{5}}{M_{\text{in}}^{3}}\right)^{1/2},$
(2.4)

is the SM temperature at which PBHs completely evaporate. Additionally,
$M_{P}\equiv 1/\sqrt{8\pi\,G}$ denotes the reduced Planck mass, and
$g_{\star}(T)$ corresponds to the number of relativistic degrees of freedom
contributing to the SM energy density. It is interesting to note that if
$\beta>\beta_{c}$, PBHs start to dominate the total energy density of the
Universe at $T=T_{\text{eq}}$, defined as
$\rho_{R}(T_{\text{eq}})=\rho_{\text{BH}}(T_{\text{eq}})$, with

$T_{\text{eq}}=\beta\,T_{\text{in}}\left(\frac{g_{\star
s}(T_{\text{in}})}{g_{\star s}(T_{\text{eq}})}\right)^{1/3},$ (2.5)

with $g_{\star s}(T)$ being the number of relativistic degrees of freedom
contributing to the SM entropy.

Several bounds exist in the PBH parameter space spanned by the initial
fraction $\beta$ and mass $M_{\text{in}}$ [1, 2]. Here, we will focus on PBHs
that evaporated before BBN, thus having masses smaller than $\sim 10^{9}$ g.
Although there exist constraints on such light PBHs, they are typically model
dependent [1, 2, 23]. Nevertheless, recent constraints have been derived after
considering the GWs emitted from the Hawking evaporation. In particular, a
backreaction problem can be avoided if the energy contained in GWs never
overtakes the one of the background universe [50]. More importantly, a
modification of BBN predictions due to the energy density stored in GWs can be
avoided if [51]

$\beta\lesssim 1.1\times
10^{-6}\left(\frac{\gamma}{0.2}\right)^{-\frac{1}{2}}\left(\frac{M_{\text{in}}}{10^{4}~{}\text{g}}\right)^{-\frac{17}{24}}.$
(2.6)

Next, we briefly describe the time evolution properties of Schwarzschild and
Kerr PBHs relevant for our purposes. Further details can be found in Refs.
[67, 16, 15].

###  2.1 Schwarzschild Black Holes

For the most simple BH scenario, PBHs are only characterized by their mass, so
that their horizon temperature $T_{\text{BH}}$ is determined according to [63]

$T_{\text{BH}}=\frac{M_{P}^{2}}{M_{\text{BH}}}\simeq 10^{13}\,{\rm
GeV}\left(\frac{1\,{\rm g}}{M_{\text{BH}}}\right).$ (2.7)

All particles with masses smaller than $T_{\text{BH}}$ can be emitted via
Hawking radiation. Within a time interval $dt$ and momentum $[p,\,p+dp]$, the
energy spectrum of an emitted species $i$ with spin $s_{i}$, mass $\mu_{i}$
and internal degrees of freedom $g_{i}$ can be described as [63, 15]

$\frac{\text{d}^{2}\mathcal{N}_{i}}{\text{~{}d}p~{}\text{d}t}=\frac{g_{i}}{2\pi^{2}}\,\frac{\sigma_{s_{i}}\left(M_{\text{BH}},\mu_{i},p\right)}{\exp\left[E_{i}(p)/T_{\text{BH}}\right]-(-1)^{2s_{i}}}\,\frac{p^{3}}{E_{i}(p)}\,,$
(2.8)

with $E_{i}(p)=\sqrt{p^{2}+\mu_{i}^{2}}$ denoting the energy and
$\sigma_{s_{i}}$ describing the absorption cross-section. Summing over all
possible particle species and integrating over the phase space in Eq. (2.8),
one can obtain the BH mass evolution, which is shown to be [68, 69]

$\frac{\text{d}M_{\text{BH}}}{\text{d}t}\equiv\sum_{i}\left.\frac{\text{d}M_{\text{BH}}}{\text{d}t}\right|_{i}=-\sum_{i}\int_{0}^{\infty}E_{i}\frac{\text{d}^{2}\mathcal{N}_{i}}{\text{~{}d}p\text{~{}d}t}\text{~{}d}p=-\varepsilon\left(M_{\text{BH}}\right)\frac{M_{P}^{4}}{M_{\text{BH}}^{2}}\,,$
(2.9)

with the mass evaporation function
$\varepsilon\left(M_{\text{BH}}\right)\equiv\sum_{i}g_{i}\,\varepsilon_{i}(z_{i})$,
where the contribution per degree of freedom $\varepsilon_{i}(z_{i})$ is given
by [15]

$\varepsilon_{i}(z_{i})=\frac{27}{128\pi^{3}}\int_{z_{i}}^{\infty}\frac{\psi_{s_{i}}\left(x,z_{i}\right)\left(x^{2}-z_{i}^{2}\right)}{\exp(x)-(-1)^{2s_{i}}}\,x\text{~{}d}x\,,$
(2.10)

where $x=E_{i}/T_{\text{BH}}$, $z_{i}=\mu_{i}/T_{\text{BH}}$, and
$\psi_{s_{i}}(E,\mu)\equiv\sigma_{s_{i}}(E,\mu)/(27\pi\,G^{2}\,M_{\text{BH}}^{2})$.

###  2.2 Kerr Black Holes

For the case of BHs with a nonzero spin, the horizon temperature has an
additional dependence on the spin parameter $a_{\star}\in[0,1]$ (with
$a_{\star}=0$ corresponding to the Schwarzschild limit)

$T_{\text{BH}}=\frac{2M_{P}^{2}}{M_{\text{BH}}}\,\frac{\sqrt{1-a_{\star}}}{1+\sqrt{1-a_{\star}}}\,,$
(2.11)

where $a_{\star}\equiv 8\pi\,J\,M_{P}^{2}/M_{\text{BH}}^{2}$, with $J$
denoting the total BH angular momentum. The energy spectrum is modified by the
presence of the angular momentum, introducing an explicit dependence on the
total $l$ and axial $m$ angular quantum numbers,

$\frac{\text{d}^{2}\mathcal{N}_{i}}{\text{~{}d}p\text{~{}d}t}=\frac{g_{i}}{2\pi^{2}}\sum_{l=s_{i}}\sum_{m=-l}^{l}\frac{\sigma_{s_{i}}^{lm}\left(M_{\text{BH}},p,a_{\star}\right)}{\exp\left[\left(E_{i}-m\,\Omega\right)/T_{\text{BH}}\right]-(-1)^{2s_{i}}}\,\frac{p^{3}}{E_{i}}\,,$
(2.12)

with
$\Omega=\left(4\pi\,a_{\star}\,M_{P}^{2}/M_{\text{BH}}\right)\left(1+\sqrt{1-a_{\star}^{2}}\right)^{-1}$
horizon angular velocity. The time evolution for the BH spin and mass is
described by the following set of equations [67]

$\displaystyle\frac{dM_{\text{BH}}}{dt}$
$\displaystyle=-\varepsilon\left(M_{\text{BH}},a_{\star}\right)\frac{M_{P}^{4}}{M_{\text{BH}}^{2}}\,,$
(2.13a) $\displaystyle\frac{da_{\star}}{dt}$
$\displaystyle=-a_{\star}\left[\gamma\left(M_{\text{BH}},a_{\star}\right)-2\varepsilon\left(M_{\text{BH}},a_{\star}\right)\right]\frac{M_{P}^{4}}{M_{\text{BH}}^{3}}\,,$
(2.13b)

where the evaporation functions
$\varepsilon\left(M_{\text{BH}},a_{\star}\right)$ and
$\gamma\left(M_{\text{BH}}\right)$ are given by

$\displaystyle\varepsilon_{i}\left(z_{i},a_{\star}\right)$
$\displaystyle=\frac{27}{128\pi^{3}}\int_{z_{i}}^{\infty}\sum_{l,m}\frac{\psi_{s_{i}}^{lm}\left(x,a_{\star}\right)\left(x^{2}-z_{i}^{2}\right)}{\exp\left(x^{\prime}/2f\left(a_{\star}\right)\right)-(-1)^{2s_{i}}}\,x\text{~{}d}x\,,$
(2.14a) $\displaystyle\gamma_{i}\left(z_{i},a_{\star}\right)$
$\displaystyle=\frac{27}{16\pi^{2}}\int_{z_{i}}^{\infty}\sum_{l,m}\frac{m\,\psi_{s_{i}}^{lm}\left(x,a_{\star}\right)\left(x^{2}-z_{i}^{2}\right)}{\exp\left(x^{\prime}/2f\left(a_{\star}\right)\right)-(-1)^{2s_{i}}}\,\text{d}x\,,$
(2.14b)

with $x^{\prime}=x-m\Omega^{\prime}$, being111Note that the definition of $x$
is the same as in the Schwarzschild case. $x=M_{\text{BH}}E_{i}/M_{P}^{2}$,
$\Omega^{\prime}=M_{\text{BH}}\Omega/M_{P}^{2}$, and
$f(a_{\star})=\sqrt{1-a_{\star}}/(1+\sqrt{1-a_{\star}})$. This set of time
evolution equations are the basis for describing the effect of a PBH dominated
Universe on the ALP and axion DM genesis.

##  3 Background Evolution

The evolution of the SM entropy density $s(T)=\frac{2\pi^{2}}{45}g_{\star
s}(T)\,T^{3}$ can be tracked via the Boltzmann equation

$\frac{ds}{dt}+3\,H\,s=-\frac{1}{T}\,\frac{\rho_{\text{BH}}}{M_{\text{BH}}}\left.\frac{dM_{\text{BH}}}{dt}\right|_{\text{SM}},$
(3.1)

where $H^{2}=(\rho_{R}+\rho_{\text{BH}})/(3M_{P}^{2})$, with
$\rho_{R}(T)=\frac{\pi^{2}}{30}g_{\star}(T)\,T^{4}$, and
$\frac{-1}{M_{\text{BH}}}\left.\frac{dM_{\text{BH}}}{dt}\right|_{\text{SM}}$
corresponds to the time-dependent PBH evaporation rate into SM particles. We
emphasize that Eq. (3.1) has to be numerically solved together with Eq. (2.9)
in the case of a Schwarzschild PBH, or with Eqs. (2.13) in the case of a Kerr
PBH, to extract the dynamics of the background, and in particular the
evolution of the SM temperature and the (non-conserved) SM entropy.

In a SM radiation dominated scenario, the Hubble expansion rate takes the
simple form

$H(T)=H_{R}(T)\equiv\sqrt{\frac{\rho_{R}(T)}{3\,M_{P}^{2}}}=\frac{\pi}{3}\sqrt{\frac{g_{\star}(T)}{10}}\,\frac{T^{2}}{M_{P}}\,.$
(3.2)

However, PBHs can have a strong impact on the evolution of the background
dynamics [70]. With an early PBH dominated epoch, cosmology can be
characterized by four distinct regimes, where the Hubble expansion rate is
given by [42]

$H(T)\simeq\begin{dcases}H_{R}(T)&\text{for }T\geq T_{\text{eq}}\,,\\\
H_{R}(T_{\text{eq}})\left[\frac{g_{\star s}(T)}{g_{\star
s}(T_{\text{eq}})}\left(\frac{T}{T_{\text{eq}}}\right)^{3}\right]^{1/2}&\text{for
}T_{\text{eq}}\geq T\geq T_{\text{c}}\,,\\\
H_{R}(T_{\text{ev}})\left[1-\frac{720}{\pi\,g_{\star}(T_{\text{in}})}\frac{M_{\text{in}}^{3}}{M_{P}^{4}}\frac{H_{R}^{2}(T_{\text{ev}})-H_{R}^{2}(T)}{H_{R}(T_{\text{ev}})}\right]&\text{for
}T_{\text{c}}\geq T\geq T_{\text{ev}}\,,\\\ H_{R}(T)&\text{for
}T_{\text{ev}}\geq T\,,\end{dcases}$ (3.3)

which can be simplified to

$H(T)\simeq\begin{dcases}H_{R}(T)&\text{ for }T\geq T_{\text{eq}}\,,\\\
H_{R}(T_{\text{eq}})\left(\frac{T}{T_{\text{eq}}}\right)^{3/2}&\text{ for
}T_{\text{eq}}\gg T\gg T_{\text{c}}\,,\\\
H_{R}(T_{\text{ev}})\left(\frac{T}{T_{\text{ev}}}\right)^{4}&\text{ for
}T_{\text{c}}\gg T\gg T_{\text{ev}}\,,\\\ H_{R}(T)&\text{ for
}T_{\text{ev}}\geq T\,,\end{dcases}$ (3.4)

for being more conveniently used in next expressions. Here, we have introduced
the temperature scale $T_{\text{c}}$, from which the SM radiation does not
scale as free radiation (i.e., $\rho_{R}(R)\propto R^{-1}$ with $R$ being the
scale factor) due to the entropy injected by the PBH evaporation. It is given
by [42]

$T_{\text{c}}\simeq\left[\frac{g_{\star}(T_{\text{in}})\,\pi}{5760}\,\frac{M_{P}^{10}\,T_{\text{eq}}}{M_{\text{in}}^{6}}\right]^{1/5}\simeq\left(T_{\text{eq}}\,T_{\text{ev}}^{4}\right)^{1/5}.$
(3.5)

Additionally, during its evaporation, PBHs radiate SM particles and therefore
dilute all previously produced species. The entropy injection factor is [42]

$\frac{S(T)}{S(T_{\text{ev}})}\simeq\begin{dcases}\frac{g_{\star
s}(T_{\text{eq}})}{g_{\star
s}(T_{\text{ev}})}\,\frac{g_{\star}(T_{\text{ev}})}{g_{\star}(T_{\text{eq}})}\,\frac{T_{\text{ev}}}{T_{\text{eq}}}&\text{for
}T\geq T_{\text{c}}\,,\\\ \frac{g_{\star s}(T)}{g_{\star
s}(T_{\text{ev}})}\left(\frac{T}{T_{\text{ev}}}\right)^{3}\left[1-\frac{720}{\pi\,g_{\star}(T_{\text{in}})}\frac{M_{\text{in}}^{3}}{M_{P}^{4}}\frac{H_{R}^{2}(T_{\text{ev}})-H_{R}^{2}(T)}{H_{R}(T_{\text{ev}})}\right]^{-2}&\text{for
}T_{\text{c}}\geq T\geq T_{\text{ev}}\,,\\\ 1&\text{for }T_{\text{ev}}\geq
T\,,\end{dcases}$ (3.6)

which can be simplified to

$\frac{S(T)}{S(T_{\text{ev}})}\simeq\begin{dcases}\frac{T_{\text{ev}}}{T_{\text{eq}}}&\text{
for }T\geq T_{\text{c}}\,,\\\ \left(\frac{T_{\text{ev}}}{T}\right)^{5}&\text{
for }T_{\text{c}}\geq T\geq T_{\text{ev}}\,,\\\ 1&\text{ for
}T_{\text{ev}}\geq T\,,\end{dcases}$ (3.7)

again to ease further analytical estimations. We note that in Eqs. (3.6) and
(3.7) there are three distinctive regimes, and not four as in Eq. (3.3).

Finally, it is important to emphasize that all these analytical estimations
help to understand the dynamics of the cosmological evolution. However,
hereafter full numerical solutions of the background will be used. In our
code, we solve such Boltzmann equations together with the time evolution
equations for the PBHs, given in Eqs. (2.13). Notice that taking as initial
value $a_{\star}=0$, the evolution equations in (2.13) reduce to the ones for
the Schwarzschild scenario, i.e., Eq. (2.9). In other words, such scenario is
included in the more general Kerr case. Thus, we take in general the equations
for the Kerr case, putting $a_{\star}=0$ whenever we talk about Schwarzschild
PBHs. Besides, in the evaporation functions
$\varepsilon\left(M_{\text{BH}},a_{\star}\right)$ and
$\gamma\left(M_{\text{BH}}\right)$ we have included an additional, almost
massless, pseudoscalar degree of freedom, corresponding to the ALP or the
axion.

##  4 ALP Dark Matter

In this section, we focus on a general light pseudoscalar, namely the ALP [43,
44, 45]. Similarly to the QCD axion, ALPs could arise as consequence of the
spontaneous breaking of a global $U(1)$ symmetry or, alternatively, they could
emerge from string theory [47, 48].

In the so-called misalignment mechanism, the present energy density stored in
the zero mode of an ALP field $a$ of mass $m_{a}$ can be obtained by solving
the equation of motion [29, 71, 30, 31]

$\ddot{\theta}+3\,H(t)\,\dot{\theta}+m_{a}^{2}\,\sin\theta=0\,,$ (4.1)

where $\theta(t)\equiv a(t)/f_{a}$, and $f_{a}$ is the Peccei-Quinn (PQ)
symmetry breaking energy scale. In this scenario, the ALP field starts to roll
about the minimum of the potential once the Hubble friction is overcome by the
potential term. Coherent oscillations of the ALPs are set around the
temperature $T_{\text{osc}}$ defined as

$3\,H(T_{\text{osc}})\equiv m_{a}\,.$ (4.2)

In the standard cosmological scenario, as the SM entropy is conserved, the ALP
energy density $\rho_{a}$ at present is given by

$\rho_{a}(T_{0})=\rho_{a}(T_{\text{osc}})\,\frac{s(T_{0})}{s(T_{\text{osc}})}\,,$
(4.3)

where $T_{0}$ is the CMB temperature at present, and we considered the
conservation of the ALP number density in a comoving volume. Additionally,
$\rho_{a}(T_{\text{osc}})\simeq\frac{1}{2}m_{a}^{2}\,f_{a}^{2}\,\theta_{i}^{2}$
in the limit in which the kinetic energy is neglected, and for a quadratic
potential. Moreover, $\theta_{i}$ denotes the initial misalignment angle. The
entropy injection from PBHs evaporation would inevitably dilute the energy
density so that one has

$\rho_{a}(T_{0})=\rho_{a}(T_{\text{osc}})\frac{s(T_{0})}{s(T_{\text{osc}})}\times\frac{S(T_{\text{osc}})}{S(T_{\text{ev}})}\simeq\frac{1}{2}m_{a}^{2}\,f_{a}^{2}\,\theta_{i}^{2}\,\frac{s(T_{0})}{s(T_{\text{osc}})}\times\frac{S(T_{\text{osc}})}{S(T_{\text{ev}})}\,,$
(4.4)

where $S(T)=s(T)\,R^{3}$ corresponds to the total SM entropy. With such an
energy density, one can compute the ALP DM relic density via

$\Omega_{a}h^{2}\equiv\frac{\rho_{a}(T_{0})}{\rho_{c}/h^{2}}\,,$ (4.5)

where $\rho_{c}/h^{2}\simeq 1.1\times 10^{-5}$ GeV/cm3 is the critical energy
density and $s(T_{0})\simeq 2.69\times 10^{3}$ cm-3 is the entropy density at
present, which should match the total DM abundance $\Omega h^{2}\simeq 0.12$
[72].

Knowing the Hubble expansion rate given in the previous section, the ALP
oscillation temperature $T_{\text{osc}}$ can be estimated to be

$T_{\text{osc}}\simeq\begin{dcases}\left(\frac{1}{\pi}\sqrt{\frac{10}{g_{\star}}}\,M_{P}\,m_{a}\right)^{1/2}&\text{for
}T\geq T_{\text{eq}}\,,\\\
\left(\frac{1}{\pi}\sqrt{\frac{10}{g_{\star}}}\,\frac{M_{P}\,m_{a}}{\sqrt{T_{\text{eq}}}}\right)^{2/3}&\text{for
}T_{\text{eq}}\geq T\geq T_{\text{c}}\,,\\\
\left(\frac{1}{\pi}\sqrt{\frac{10}{g_{\star}}}\,M_{P}\,m_{a}\,T_{\text{ev}}^{2}\right)^{1/4}&\text{for
}T_{\text{c}}\geq T\geq T_{\text{ev}}\,,\\\
\left(\frac{1}{\pi}\sqrt{\frac{10}{g_{\star}}}\,M_{P}\,m_{a}\right)^{1/2}&\text{for
}T_{\text{ev}}\geq T\,.\end{dcases}$ (4.6)

Using Eq. (4.4), one can compute the relic density in the aforementioned four
regimes, which becomes

$\Omega_{a}h^{2}\simeq\frac{45\,\theta_{i}^{2}}{4\pi^{2}\,g_{\star
s}}\,\frac{s(T_{0})}{\rho_{c}/h^{2}}\times\begin{dcases}\left(\sqrt{\frac{g_{\star}}{10}}\,\frac{\pi}{M_{P}}\right)^{3/2}\frac{T_{\text{ev}}}{T_{\text{eq}}}\,m_{a}^{1/2}\,f_{a}^{2}&\text{
for }T_{\text{osc}}\geq T_{\text{eq}}\,,\\\
\left(\sqrt{\frac{g_{\star}}{10}}\,\frac{\pi}{M_{P}}\right)^{2}T_{\text{ev}}\,f_{a}^{2}&\text{
for }T_{\text{eq}}\geq T_{\text{osc}}\geq T_{\text{ev}}\,,\\\
\left(\sqrt{\frac{g_{\star}}{10}}\,\frac{\pi}{M_{P}}\right)^{3/2}m_{a}^{1/2}\,f_{a}^{2}&\text{
for }T_{\text{ev}}\geq T_{\text{osc}}\,,\end{dcases}$ (4.7)

or equivalently

$\frac{\Omega_{a}h^{2}}{0.12}\simeq\begin{dcases}\left(\frac{\theta_{i}}{0.5}\right)^{2}\left(\frac{T_{\text{ev}}}{4~{}\text{MeV}}\right)\left(\frac{81~{}\text{GeV}}{T_{\text{eq}}}\right)\left(\frac{m_{a}}{10^{-3}~{}\text{eV}}\right)^{1/2}\left(\frac{f_{a}}{5\times
10^{14}~{}\text{GeV}}\right)^{2}&\text{for }T_{\text{osc}}\geq
T_{\text{eq}},\\\
\left(\frac{\theta_{i}}{0.5}\right)^{2}\left(\frac{T_{\text{ev}}}{4~{}\text{MeV}}\right)\,\left(\frac{f_{a}}{10^{15}~{}\text{GeV}}\right)^{2}&\text{for
}T_{\text{eq}}\geq T_{\text{osc}}\geq T_{\text{ev}},\\\
\left(\frac{\theta_{i}}{0.5}\right)^{2}\left(\frac{m_{a}}{10^{-15}~{}\text{eV}}\right)^{1/2}\,\left(\frac{f_{a}}{3\times
10^{15}~{}\text{GeV}}\right)^{2}&\text{for }T_{\text{ev}}\geq
T_{\text{osc}}.\end{dcases}$ (4.8)

Note that the expressions for $\Omega_{a}h^{2}$ are identical for
$T_{\text{eq}}\geq T_{\text{osc}}\geq T_{\text{c}}$ and $T_{\text{c}}\geq
T_{\text{ev}}$. Considering the full numerical solution of the Boltzmann
equations, we present in Fig. 1 the parameter space generating the whole
observed ALP DM abundance for benchmark values of $m_{a}=10^{-7}$ eV and
$f_{a}=10^{14}$ GeV, while assuming Schwarzschild ($a_{\star}=0$, red band) or
Kerr ($a_{\star}=0.999$, blue band) PBHs. The thickness of the bands
corresponds to an initial misalignment angle
$0.5\leq\theta_{i}\leq\pi/\sqrt{3}$.222The lower bound comes from assuming an
$\mathcal{O}(1)$ misalignment angle that avoids introducing fine tuning of the
initial conditions, whereas the upper bound is an effective average
misalignment angle, typical from post-inflationary scenarios. The regions on
the right (left) of the bands produce a DM underabundance (overabundance).
Additionally, the gray regions are in tension with BBN ($T_{\text{ev}}<4$ MeV
[73, 74, 75]), whereas the green region with GWs (i.e., Eq. (2.6)). Finally,
the dashed red line corresponds to $\beta=\beta_{c}$, limiting the region
where PBH energy density is subdominant with respect to SM radiation, and
therefore one has a standard cosmological scenario (below), with the region
where a nonstandard cosmological expansion triggered by PBHs occurs (above).
It is important to emphasize that these results where obtained using the
numerical code developed in Ref. [16, 15], including the additional ALP degree
of freedom. We have checked that the numerical estimations in Eq. (4.7) fit
well the analytical results.

Figure 1: Parameter space compatible with the whole observed ALP DM abundance,
for $m_{a}=10^{-7}$ eV and $f_{a}=10^{14}$ GeV. The thickness of the bands
correspond to $0.5\leq\theta_{i}\leq\pi/\sqrt{3}$. The red bands show
$a_{\star}=0$, whereas the blue $a_{\star}=0.999$.

As seen in Fig. 1, the presence of an era dominated by PBHs can have a strong
impact on the ALP production in the early Universe. The maximal deviation from
the standard cosmological scenario corresponds to long-lived PBHs evaporating
just before the onset of BBN, and the maximal value for $\beta$ allowed by
GWs. This happens for $M_{\text{in}}\simeq 5.7\times 10^{8}$ g and
$\beta\simeq 4.6\times 10^{-10}$ for $a_{\star}=0$, or $M_{\text{in}}\simeq
7.6\times 10^{8}$ g and $\beta\simeq 3.5\times 10^{-10}$ for
$a_{\star}=0.999$, and corresponds to the upper right white corners in Fig. 1.
This maximal deviation from the standard cosmological scenario is shown in
Fig. 2, for the benchmark $M_{\text{in}}\simeq 5.7\times 10^{8}$ g and
$\beta\simeq 4.6\times 10^{-10}$, for $a_{\star}=0$. The left panel shows the
evolution of the radiation and PBH energy densities as a function of the SM
temperature. PBHs dominate the total energy density between
$T=T_{\text{eq}}\simeq 81$ GeV and $T=T_{\text{ev}}\simeq 4$ MeV.
Additionally, the right panel shows the oscillation temperature for radiation
domination (solid), and PBH domination (dotted). As expected from Eq. (4.6),
for a given ALP mass, $T_{\text{osc}}$ decreases for a PBH dominated era.
Again, these results were obtained using the numerical code and fit well the
analytical estimations.

Figure 2: Left: Evolution of the energy densities for SM radiation and PBHs.
Right: Oscillation temperature of ALPs, for radiation domination (solid) and
PBH domination (dotted). For both panels, the maximal allowed deviation from
the standard cosmological scenario for Schwarzschild PBHs was assumed, i.e.
$M_{\text{in}}\simeq 5.7\times 10^{8}$ g and $\beta\simeq 4.6\times 10^{-10}$.

It is customary to assume that the global PQ symmetry is also anomalous
respect to the electromagnetic gauge group thank to the existence of exotic
vector-like charged fermions. Consequently, ALPs will couple to two photons,
via the effective dimension-5 operator

$\mathcal{L}_{a\gamma}=-\frac{1}{4}\,g_{a\gamma}\,a\,F_{\mu\nu}\,\tilde{F}^{\mu\nu}=g_{a\gamma}\,a\,\vec{E}\cdot\vec{B}\,,$
(4.9)

where the coupling constant $g_{a\gamma}$ is model dependent and related to
the breaking scale of the PQ symmetry as

$g_{a\gamma}=\frac{\alpha}{2\pi\,f_{a}}\left(\frac{E}{N}-\frac{2}{3}\,\frac{4+z}{1+z}\right)\simeq
10^{-13}~{}\text{GeV}^{-1}\left(\frac{10^{10}~{}\text{GeV}}{f_{a}}\right),$
(4.10)

with $z\equiv m_{u}/m_{d}$, and $E$ and $N$ are the electromagnetic and color
anomalies associated with the ALP anomaly. For KSVZ models $E/N=0$ [76, 77],
whereas for DFSZ models $E/N=8/3$ [78, 79]. The electromagnetic interaction of
ALPs is by far the most exploited to look for signatures in observations and
experimental searches [80].

Once one assumes that all DM is composed by ALPs, $m_{a}$ and $f_{a}$ (and
therefore $g_{a\gamma}$) are not longer independent, for a given initial
misalignment angle. In Fig. 3 we display the parameter space compatible with
the DM relic density constraint, in the plane $(m_{a},\,|g_{a\gamma}|)$. The
diagonal band dubbed “ALPs in radiation domination” corresponds to the viable
region within the standard cosmology, with the thickness of the band
representing the possible values for $\theta_{i}$ in the range
$[0.5,\,\pi/\sqrt{3}]$, with $\theta_{i}=\pi/\sqrt{3}$ being the top border.
The regions on the top and bottom give rise to a DM underabundance and
overabundance, respectively.333Notice that we are not considering further
production mechanisms of ALPs such as the decay of cosmic strings or domain
walls [81, 82, 83, 84]. This parameter space is modified once an early PBH
dominated epoch occurs. The red (blue) thick band shows the viable parameter
space for $a_{\star}=0$ ($a_{\star}=0.999$), enhanced by the presence of PBHs.
The impact of the PBH spin is very mild, and therefore the two bands are
almost superimposed. The PBH entropy injection dilutes the ALP abundance, and
therefore higher values for the PQ breaking scale are allowed, which
corresponds to lower values for $|g_{a\gamma}|$. It is interesting to note
that for $m_{a}\lesssim 10^{-5}$ eV, $T_{\text{osc}}<T_{\text{eq}}$, and
therefore the ALP relic density is independent on its mass. However, for
$m_{a}\gtrsim 10^{-5}$ eV, $T_{\text{osc}}>T_{\text{eq}}$ and
$g_{a\gamma}\propto m_{a}^{1/4}$, cf. Eq. (4.7).

Figure 3: Viable parameter space for ALPs and the QCD axion, for radiation
dominated or PBH dominated scenarios. The brown and light blue shaded regions
are the current exclusions and the projected sensitivities of the different
experiments described in the text. The figure has been adapted from Ref. [85].

Figure 3 also overlays in brown current bounds on the ALP-photon coupling.
Small masses ($m_{a}\lesssim 4\times 10^{-11}$ eV) are constrained
astrophysical black hole spins [86], whereas high masses ($m_{a}\gtrsim
10^{2}$ eV) by the $X$-ray background and extragalactic background light [87].
The masses in the range $10^{-6}$ eV $\lesssim m_{a}\lesssim 10^{-5}$ eV are
already probed by the so-called haloscope experiments [88], ADMX [89, 90],
HAYSTAC [91], CAPP [92, 93, 94], QUAX [95, 96], ORGAN [97], using a highly
tuned microwave cavity that converts ALPs into photons in the presence of a
static magnetic field.444We note that the parameter space shown in Fig. 3 is
consistent with an ALP lifetime larger than the age of the Universe.
Additionally, the light blue regions show projected sensitivities from a
number of experiments that count ORGAN [97], MADMAX [98], ALPHA [99], ADMX
[100], KLASH [101], DM-Radio [102] and ABRACADABRA [103, 104]. These bounds
have been adapted from Ref. [85]. Even if the PBH domination tends to imply
smaller values for $|g_{a\gamma}|$, this parameter space could be potentially
tested in the future by next-generation experiments, especially for light ALP
masses, $m_{a}\lesssim 10^{-7}$ eV.

Before closing this section, it is interesting to note that further
constraints on $f_{a}$ appear depending on the scale at which the spontaneous
breakdown of the PQ symmetry happened. If the PQ symmetry was broken during
inflation ($f_{a}>H_{I}$, with $H_{I}$ being the inflationary scale) and never
restored, the ALP field would be homogenized through the Hubble patch. On the
other hand, if the inflationary epoch ends before the PQ transition
($f_{a}<H_{I}$) the initial misalignment angle would vary along the different
patches of the Universe. Let us consider in more detail these two scenarios
and their possible relation with PBHs.

#### Pre-inflationary scenario

In this case, the PQ symmetry is spontaneously broken during inflation (i.e.
$H_{I}<f_{a}$), and it is not restored afterwards [31]. The ALP field is
homogeneous through various Hubble patches, with a unique value of
$\theta_{i}$ characterizing the whole observable Universe. In this scenario,
the ALP is present during inflation and therefore ALP isocurvature
fluctuations (converted into curvature perturbations) are expected to leave a
imprint on the CMB. However, since the CMB measurements do not allow for
sizable isocurvature modes, the scale of inflation is pushed to relatively low
values. Accordingly, the isocurvature bounds obtained from Planck data [105]
impose the lower bound on $f_{a}$ which can be cast as [34]

$H_{I}\lesssim 0.9\times
10^{7}~{}\mathrm{GeV}\left(\frac{\theta_{i}}{\pi}\right)\left(\frac{f_{a}}{10^{11}~{}\mathrm{GeV}}\right).$
(4.11)

Since the BBN scale represents a lower bound for $H_{I}$, it follows that
$f_{a}\gtrsim 100$ GeV ($|g_{a\gamma}|\lesssim 10^{-5}$ GeV-1) for
$\theta_{i}\sim\mathcal{O}(1)$, which in turn does not have an impact on the
parameter space displayed in Fig. 3. A stronger limit appears when examining
the connection with PBHs. Combining the above expression with Eq. (2.1) we
obtain a lower limit for the initial mass of the PBH

$M_{\text{in}}=4\pi\,\gamma\,\frac{M_{P}^{2}}{H(T_{\text{in}})}\gtrsim
4\pi\,\gamma\,\frac{M_{P}^{2}}{H_{I}}\simeq 8.6\times
10^{6}~{}\mathrm{g}\left(\frac{1}{\theta_{i}}\right)\left(\frac{10^{11}~{}\mathrm{GeV}}{f_{a}}\right),$
(4.12)

which reflects the fact that PBHs are created after inflation, in a radiation-
dominated epoch. For PBHs fully evaporating before the onset of BBN,
$M_{\text{in}}\lesssim 6\times 10^{8}$ g, and therefore for
$\theta_{i}\sim\mathcal{O}(1)$, it implies that

$|g_{a\gamma}|\lesssim 10^{-12}~{}\textrm{GeV}^{-1},$ (4.13)

which is again automatically satisfied for our relevant parameter space.

#### Post-inflationary scenario

For this case, the PQ symmetry is spontaneously broken after the inflationary
period, i.e. $f_{a}<H_{I}$. The misalignment angle takes random values along
different Hubble patches, in which case it is averaged out over many patches
(in the range $[-\pi,\,\pi)$) so that
$\theta_{i}=\sqrt{\langle\theta_{i}^{2}\rangle}=\pi/\sqrt{3}$.555When taking
into account the anharmonicity of the potential, the average becomes
$\theta_{i}=\sqrt{\langle\theta_{i}^{2}\rangle}\simeq 2.15$ [106].
Additionally, the upper limit on the inflationary scale $H_{I}<2.5\times
10^{-5}~{}M_{P}$ [105] implies that $f_{a}\lesssim 10^{13}$ GeV (and in turn
$g_{a\gamma}\gtrsim 10^{-16}$ GeV-1), which finally can be translated into an
lower bound on the ALP mass $m_{a}\gtrsim 10^{-7}$ eV, for ALPs being the
whole DM.

An important feature of the post-inflationary scenario is that the value of
the initial misalignment angle changes by a factor of $\mathcal{O}(1)$ from
one causal patch to the next. Accordingly, the density of cold ALPs is
characterized by sizable inhomogeneities. Their free streaming length is too
short to erase these inhomogeneities before the matter-radiation equality, so
that the density perturbations decouple from the Hubble expansion and start
growing by gravitational instability, rapidly forming gravitationally bound
objects, called miniclusters [54, 56, 55, 57, 34].

The formation time of ALP miniclusters is sensitive to the cosmology prior to
BBN and, in particular, to the temperature at which ALP oscillations take
place. The scale of ALP minicluster mass is set by the total mass of ALPs
within one Hubble volume of radius $R_{\text{osc}}\sim
H(T_{\text{osc}})^{-1}$, at the time of the matter-radiation equality.
Therefore, their mass at formation is

$M_{0}=\frac{4\pi}{3}\,\rho_{\rm
DM}(T_{0})\,\frac{s(T_{\text{osc}})}{s(T_{0})}\left(\frac{1}{H(T_{\text{osc}})}\right)^{3}\times\frac{S(T_{\text{ev}})}{S(T_{\text{osc}})}\,.$
(4.14)

In a radiation dominated Universe and assuming that all relic ALPs bound up in
miniclusters, it reduces to

$M_{0}=\frac{8\,g_{\star
s}(T_{\text{osc}})}{5}\left(\frac{10}{g_{\star}(T_{\text{osc}})}\right)^{3/2}\frac{\Omega_{\rm
DM}h^{2}\,\rho_{c}/h^{2}}{s(T_{0})}\,\frac{M_{P}^{3}}{T_{\text{osc}}^{3}}\simeq
2\times
10^{-16}\left(\frac{10^{-5}~{}\text{eV}}{m_{a}}\right)^{3/2}M_{\odot}\,,$
(4.15)

as shown in Fig. 4 with a solid black line. The vertical bands are excluded
because $f_{a}>H_{I}$, or in tension with the $X$-ray background and the
extragalactic background light [87].

Figure 4: Minicluster mass in the post-inflationary scenario, at the time of
formation. ALPs: The thick black line corresponds to the the standard
cosmological scenario, whereas the purple band to the PBH domination. Axions:
The star corresponds to the the standard cosmological scenario, whereas the
purple line to the PBH domination. The vertical bands are excluded because
$f_{a}>H_{I}$, or in tension with $X$-rays measurements.

However, if ALPs start to oscillate during the PBH dominated period, the
initial minicluster mass can be estimated to be666It is worth mentioning that
in either a radiation or early matter domination era may be possible the
formation of ALP miniclusters even in pre-inflationary scenarios [107, 108].

$M_{0}\simeq\frac{16\pi\,g_{\star s}}{g_{\star}}\,\frac{\Omega_{\rm
DM}h^{2}\,\rho_{c}/h^{2}}{s(T_{0})}\times\begin{dcases}\pi^{1/2}\left(\frac{g_{\star}}{10}\right)^{1/4}\left(\frac{M_{P}}{m_{a}}\right)^{3/2}\frac{T_{\text{eq}}}{T_{\text{ev}}}&\text{for
}T_{\text{osc}}\geq T_{\text{eq}}\,,\\\
\frac{M_{P}^{2}}{m_{a}\,T_{\text{ev}}}&\text{for }T_{\text{eq}}\geq
T_{\text{osc}}\geq T_{\text{ev}}\,,\\\
\pi^{1/2}\left(\frac{g_{\star}}{10}\right)^{1/4}\left(\frac{M_{P}}{m_{a}}\right)^{3/2}&\text{for
}T_{\text{ev}}\geq T_{\text{osc}}\,.\end{dcases}$ (4.16)

Figure 4 also shows the maximal impact of the PBH domination era with a broad
purple band. Again, we have taken the PBH benchmark values of
$M_{\text{in}}\simeq 5.7\times 10^{8}$ g and $\beta\simeq 4.6\times 10^{-10}$
for Schwarzschild, or $M_{\text{in}}\simeq 7.6\times 10^{8}$ g and
$\beta\simeq 3.5\times 10^{-10}$ for Kerr. The spin of the PBHs has a very
limited impact, and therefore the differences in the plot are barely visible.
However, it is clear that heavier ALP miniclusters can be formed due to the
entropy injection. In other words, the entropy injection works as an
enhancement factor, cf. Eq. (4.14) (see e.g. Refs. [108, 109, 110] for related
studies). Hence, searches for ALPs, such as indirect detection through
gravitational microlensing [111, 112, 113, 114], could be extended to regions
of the parameter space previously thought to be not relevant.

##  5 Axion Dark Matter

The production of QCD axion in a cosmology dominated by PBHs was analytically
studied in Ref. [42]. For the sake of completeness, here we revisit that
scenario from a numerically perspective. In particular, we: $i)$ include the
greybody factors in the PBHs emission properties, $ii)$ take the effect of the
PBH angular momentum into account, and $iii)$ numerically solve the evolution
of the PBH and radiation energy density. This is particularly pertinent
because it allows to take into account the time-dependent BH evaporation rate.

For the case of QCD axions, the scale $f_{a}$ at which the PQ $U(1)$ symmetry
is spontaneously broken determines the axion mass through the topological
susceptibility of QCD [115]

$\chi(T)\simeq 0.0245~{}\text{fm}^{-4}\times\begin{dcases}1&\text{ for }T\leq
T_{\text{QCD}}\,,\\\ \left(\frac{T}{T_{\text{QCD}}}\right)^{-8.16}&\text{ for
}T\geq T_{\text{QCD}}\,,\end{dcases}$ (5.1)

as

$m_{a}(T)=\frac{\sqrt{\chi(T)}}{f_{a}}\,,$ (5.2)

with $T_{\text{QCD}}\simeq 150$ MeV. As the axion mass is temperature
dependent, its abundance in the standard cosmology is

$\Omega_{a}h^{2}\simeq
0.12\left(\frac{\theta_{i}}{10^{-3}}\right)^{2}\times\begin{dcases}\left(\frac{m_{a}}{m_{a}^{\text{QCD}}}\right)^{-\frac{3}{2}}&\text{for
}m_{a}\leq m_{a}^{\text{QCD}},\\\
\left(\frac{m_{a}}{m_{a}^{\text{QCD}}}\right)^{-\frac{7}{6}}&\text{for
}m_{a}\geq m_{a}^{\text{QCD}},\end{dcases}$ (5.3)

with $m_{a}^{\text{QCD}}\equiv m_{a}(T_{\text{osc}}=T_{\text{QCD}})\simeq
4.8\times 10^{-11}$ eV. However, if a PBH-dominated epoch occurs, the axion
relic abundance gets modified due to the reduction of its oscillation
temperature and its dilution induced by the entropy injection [116, 42, 117].
Again, we numerically solve the coupled Boltzmann equations, determine the
dilution factor and further calculate the modified QCD axion relic density. In
Fig. 5, we present our numerical results, which agrees well with the
analytical estimations presented in Ref. [42]. The purple and blue bands of
Fig. 5 show the misalignment angle required to reproduce the whole observed
axion DM abundance for the case with Schwarzschild ($a_{\star}=0$) and Kerr
($a_{\star}=0.999$) PBH domination, respectively. The black thick line
corresponds to the standard cosmological scenario. The thickness of the band
brackets all possible PBH scenarios compatible with BBN ($T_{\text{ev}}>4$
MeV) and GWs (cf. Eq. (2.6)). Once we restrict to initial misalignment angles
in the range $\theta_{i}\in[0.5,\,\pi/\sqrt{3}]$, we have that the axion mass
can take lower values than those allowed in a purely radiation dominated
epoch. Concretely, the axion mass range gets broaden towards low values from
$\sim[10^{-6},\,10^{-5}]$ eV to $\sim[10^{-8},\,10^{-5}]$ eV.

Figure 5: Initial misalignment angle of QCD axions required to explain the
observed DM abundance for $a_{\star}=0$ (red band) and $a_{\star}=0.999$ (blue
band). The corresponding values in the standard cosmological scenario are
represented by the black thick line.

This enlargement of the viable mass range is also shown in Fig. 3, now taking
into account the axion coupling to photons. The thickness of the band
correspond to the KSVZ and DFSZ models for the axion-photon coupling. It is
interesting to see that this parameter space is within the reach of future
detectors as ABRACADABRA, KLASH, ADMX, and DM-Radio.

#### Pre-inflationary scenario

When axion is present during inflation (and the PQ symmetry is not restored
after it), the requirement of axions to constitute all DM along with the
standard cosmological scenario along with the isocurvature bound (4.11) leads
to

$H_{I}\lesssim\begin{dcases}3.2\times
10^{9}~{}\mathrm{GeV}\left(\frac{f_{a}}{10^{17}~{}\mathrm{GeV}}\right)^{1/4}&\text{
for }f_{a}\gtrsim 1.2\times 10^{17}~{}\text{GeV},\\\ 2.6\times
10^{7}~{}\mathrm{GeV}\left(\frac{f_{a}}{10^{12}~{}\mathrm{GeV}}\right)^{5/12}&\text{
for }f_{a}\lesssim 1.2\times 10^{17}~{}\text{GeV}.\end{dcases}$ (5.4)

Consequently, under the PBH domination slightly higher inflation scales can be
reached since $f_{a}$ is allowed to take higher values. In other words, the
relevant parameter space displayed in Figs. 3 and 5 does not get affected.
Translating the above bounds into the initial mass for the PBHs, we have

$M_{\text{in}}\gtrsim\begin{dcases}7.5\times
10^{3}~{}\mathrm{g}\left(\frac{10^{17}~{}\mathrm{GeV}}{f_{a}}\right)^{1/4}&\text{
for }f_{a}\gtrsim 1.2\times 10^{17}~{}\text{GeV},\\\ 9.4\times
10^{5}~{}\mathrm{g}\left(\frac{10^{12}~{}\mathrm{GeV}}{f_{a}}\right)^{5/12}&\text{
for }f_{a}\lesssim 1.2\times 10^{17}~{}\text{GeV}.\end{dcases}$ (5.5)

Contrary to the ALP case, there is no significant lower bound on $f_{a}$ when
$M_{\text{in}}$ approaches to the maximum value allowed by BBN.

#### Post-inflationary scenario

Since the PQ symmetry remains conserved during the inflation period
($f_{a}<H_{I}$), the axion field does not get homogenized. The upper limit on
the inflationary scale implies that $f_{a}\lesssim 10^{13}$ GeV, and hence the
lower bound on the axion mass $m_{a}\gtrsim 6\times 10^{-7}$ eV.

In the standard cosmological scenario, axions start oscillating at the
temperature $T_{\text{osc}}\simeq 1.1$ GeV. Therefore, if the misalignment
mechanism is the unique method to produce axions, they should have a mass
$m_{a}\simeq 1.8\times 10^{-5}$ eV (for $\theta_{i}=\pi/\sqrt{3}$) in order to
account for the whole of the DM abundance. This case is identified by a star
in the upper right region in the Fig. 5. This result is remarkably modified by
the dilution factor $S(T_{\text{osc}})/S(T_{\text{ev}})$ that emerges within
the PBH domination, which now allows for a mass range instead of a single
point, reaching values up to two orders of magnitude lower, i.e.
$\mathcal{O}(10^{-7})$ eV. This new range is represented by the thick solid
horizon line in Fig. 5.

The gravitational clump of axion density inhomogeneities at the time of
radiation-matter equality leads to miniclusters having a mass $M_{0}\simeq
2.1\times 10^{-11}~{}M_{\odot}$ in the standard cosmological scenario (see the
star mark in Fig. 4). However, in a PBH dominated scenario, such a mass gets
enhanced due to the entropy injection by BH evaporation. The mass $M_{0}$ of
the axion miniclusters at formation becomes

$M_{0}\simeq\frac{16\,g_{\star s}}{g_{\star}}\,\frac{\Omega
h^{2}\,\rho_{c}/h^{2}}{s(T_{0})}\times\begin{dcases}\pi^{1/2}\left(\frac{10}{g_{\star}}\right)^{1/4}\frac{M_{P}^{5/2}}{T_{\text{QCD}}^{2}\,m_{a}^{1/2}}\frac{T_{\text{eq}}}{T_{\text{ev}}}&\text{for
}T_{\text{osc}}\geq T_{\text{eq}}\,,\\\
\pi^{3/11}\left(\frac{10}{g_{\star}}\right)^{4/11}\frac{M_{P}^{30/11}}{T_{\text{QCD}}^{12/11}\,m_{a}^{3/11}\,T_{\text{eq}}^{4/11}\,T_{\text{ev}}}&\text{for
}T_{\text{eq}}\geq T_{\text{osc}}\geq T_{\text{c}}\,,\\\
\pi^{1/2}\left(\frac{10}{g_{\star}}\right)^{1/4}\frac{M_{P}^{5/2}}{T_{\text{QCD}}^{2}\,m_{a}^{1/2}}&\text{for
}T_{\text{c}}\geq T_{\text{osc}}\,.\end{dcases}$ (5.6)

In Fig. 4 we show the enhancement on $M_{0}$ (purple line) generated by the
PBH domination. It is remarkable that now $M_{0}$ can reach values as high as
$10^{-8}~{}M_{\odot}$.

##  6 Dark Radiation

Axions and ALPs emitted from PBH evaporation are ultra-relativistic, thus
behaving as DR and contributing to the effective number of neutrinos [5, 23,
34, 41]. We note that ALP DR directly emitted from nonrotating PBHs was
analyzed in Ref. [41] recently. Here we revisit ALP DR from PBHs with a more
general perspective, namely we will study both nonrotating and rotating PBHs.
The contribution to effective number of neutrinos can be simply estimated,
obtaining

$\Delta
N_{\text{eff}}=\left[\frac{8}{7}\left(\frac{11}{4}\right)^{\frac{4}{3}}+N_{\text{eff}}^{\text{SM}}\right]\frac{\rho_{a}(T_{\text{ev}})}{\rho_{R}(T_{\text{ev}})}\left(\frac{g_{\star}(T_{\text{ev}})}{g_{\star}(T_{\text{equality}})}\right)\left(\frac{g_{\star
s}(T_{\text{equality}})}{g_{\star s}(T_{\text{ev}})}\right)^{\frac{4}{3}},$
(6.1)

where $T_{\text{equality}}\simeq 0.75$ eV denotes the temperature at (late)
matter-radiation equality. In the PBH-dominated scenario, the ratio between
the axion (or ALP) and the radiation energy densities is related to their
contributions to the evaporation function

$\frac{\rho_{a}(T_{\text{ev}})}{\rho_{R}(T_{\text{ev}})}=\frac{\varepsilon_{a}(z_{i},a_{\star})}{\varepsilon_{R}(z_{i},a_{\star})}\,,$
(6.2)

with $\varepsilon_{a}(T_{\text{ev}})$ being the axion evaporation function,
and $\varepsilon_{R}(T_{\text{ev}})$ the total SM radiation contribution to
the evaporation function. For the case of Schwarzschild PBHs, $\Delta N_{\rm
eff}$ is simply

$\Delta N_{\text{eff}}\simeq
0.237\,\frac{g_{\star}\left(T_{\text{ev}}\right)}{g_{\star
s}\left(T_{\text{ev}}\right)^{\frac{4}{3}}}\,,$ (6.3)

which is independent from $\beta$, and consistent with the result obtained in
Ref. [41]. For Kerr PBHs, the time depletion of the angular momentum makes it
more difficult to obtain a simple analytical form.

Figure 6: Contribution to the effective number of neutrinos $\Delta
N_{\text{eff}}$ coming from axions/ALPs directly emitted via Hawking radiation
for $a_{\star}=0$ (left) and $a_{\star}=0.999$ (right). The blue line
indicates the region within the reach of the future CMB-S4 experiment [59].

Left (right) panel of Fig. 6 shows the contribution to $\Delta N_{\text{eff}}$
as function of BH mass and $\beta$ for Schwarzschild (Kerr) PBHs. Let us note
that in a PBH-dominated era, $\Delta N_{\rm eff}$ tends to be independent of
$\beta$, as noted before, c.f. Eq. (6.3). On the contrary, when
$\beta<\beta_{c}$, larger $M_{\text{in}}$ is needed with smaller $\beta$ in
order to yield a fixed $\Delta N_{\text{eff}}$. We also note that when
$M_{\text{in}}\simeq 2\times 10^{7}$ g (correspondingly $T_{\text{ev}}\simeq
0.1~{}\text{GeV}$), there are some features in the lines shown in Fig. 6. This
is because the sharp change of $g_{\star}$ and $g_{\star s}$ due to QCD phase
transition [5, 23]. Comparing the resulting values of $\Delta N_{\rm eff}$ for
Schwarzschild and Kerr PBH, we observe that for the latter case the final
values are smaller than those for nonrotating PBHs. This is simply related to
the properties of Kerr PBHs: if the BH spin is close-to-maximal, the emission
of higher spin particles is enhanced in the early stages of the evaporation
[67]. Thus, the particle injection to SM radiation density is larger than in
the Schwarzschild case. Meanwhile, the axion/ALP emission is not affected by
the BH spin, so their final contribution to $\Delta N_{\rm eff}$ is in fact
reduced [118, 12]. The current upper bound on $\Delta N_{\text{eff}}$ comes
from the Planck collaboration, and reads $N_{\mathrm{eff}}=2.99\pm 0.17$ [72].
However, $\Delta N_{\text{eff}}$ due to emitted hot axions/ALPs could reach
the sensitivity of next generation CMB experiments: CMB-S4 [59] ($\Delta
N_{\mathrm{eff}}\sim 0.06$), PICO [60] ($\Delta N_{\mathrm{eff}}\sim 0.06$)
and SPT-3G/SO [61, 62] ($\Delta N_{\mathrm{eff}}\sim 0.1$), and therefore
axion/ALP DR could be tested in the near future.777It is interesting to note
that extra radiation, at the level of $0.2\lesssim\Delta
N_{\text{eff}}\lesssim 0.5$, can alleviate the tension between measurements of
the Hubble parameter at early and late times [119, 72, 5, 120, 121, 122].
However, such large levels can not be reached by PBHs.

##  7 Conclusions

PBHs are without doubt one of the most interesting objects that could exist in
nature. As the PBH energy density scale like non-relativistic matter, they can
naturally dominate the expansion rate of the Universe, triggering a
nonstandard cosmological epoch. Moreover, they also source all types of
particles, thus injecting large quantities of entropy to the primordial
plasma. Hence, these two specific properties, namely, the possibility of
dominating expansion of the Universe and the large injection of entropy, make
the phenomenology of a PBH dominated Universe unique among the possible
scenarios for a nonstandard cosmology.

In this paper, we studied the phenomenological consequences of axion-like
particles (ALPs) as dark matter candidates in an early PBH-dominated epoch. To
that end, we numerically solved the set of Boltzmann equations for the
background dynamics, carefully taking into account both the greybody factors
and the PBH angular momentum. PBHs have a strong impact on the misalignment
ALP production. On the one hand, as the Hubble expansion rate is enhanced, the
oscillation temperature decreases, which corresponds to a delay in the
oscillations. On the other hand, the entropy injection due to the PBH
evaporation has to be compensated by a larger spontaneous breaking scale value
of the Peccei-Quinn symmetry or, equivalently, by a smaller ALP-photon
coupling. An equivalent trend happens for QCD axions, with the particularity
that its standard mass window $m_{a}\simeq[10^{-6}-10^{-5}]$ eV gets broaden
to $\sim[10^{-8}-10^{-5}]$ eV (for misalignment angles of order
$\mathcal{O}(1)$), mainly due to the entropy injection. It is interesting to
note that both for ALPs and axions, the new viable parameter space is within
the projected sensitivities of detectors as ABRACADABRA, KLASH, ADMX, and DM-
Radio.

Additionally, we showed that at production the ALP/axion minicluster mass
could increase by several orders of magnitude, due to the enhancement required
by the PBH entropy injection. This could give new prospects for the indirect
detection via gravitational microlensing of ALP and axion dark matter, with
respect to the standard cosmological scenario.

Finally, we analyzed in detail the dark radiation arising from relativistic
axions and ALPs directly emitted from PBHs evaporation. For Kerr PBH, the
contribution $\Delta N_{\text{eff}}$ are smaller than those for Schwarzschild
PBHs. Interestingly, this contribution is within the projected reach of future
CMB Stage 4 experiments and could help relaxing the tension between late and
early-time Hubble determinations.

## Acknowledgments

NB received funding from the Spanish FEDER/MCIU-AEI under grant
FPA2017-84543-P, and the Patrimonio Autónomo - Fondo Nacional de
Financiamiento para la Ciencia, la Tecnología y la Innovación Francisco José
de Caldas (MinCiencias - Colombia) grant 80740-465-2020. The work of OZ is
supported by Sostenibilidad-UdeA and the UdeA/CODI Grants 2017-16286 and
2020-3317. This project has received funding/support from the European Union’s
Horizon 2020 research and innovation programme under the Marie Skłodowska-
Curie grant agreement No 860881-HIDDeN. Fermilab is operated by the Fermi
Research Alliance, LLC under contract No. DE-AC02-07CH11359 with the United
States Department of Energy.

## References

  * [1] B. Carr, K. Kohri, Y. Sendouda and J. Yokoyama, _Constraints on Primordial Black Holes_ , 2002.12778.
  * [2] B. Carr and F. Kuhnel, _Primordial Black Holes as Dark Matter: Recent Developments_ , _Ann. Rev. Nucl. Part. Sci._ 70 (2020) 355 [2006.02838].
  * [3] P. Villanueva-Domingo, O. Mena and S. Palomares-Ruiz, _A brief review on primordial black holes as dark matter_ , _Front. Astron. Space Sci._ 8 (2021) 87 [2103.12087].
  * [4] B. Carr and F. Kuhnel, _Primordial Black Holes as Dark Matter Candidates_ , 2110.02821.
  * [5] D. Hooper, G. Krnjaic and S.D. McDermott, _Dark Radiation and Superheavy Dark Matter from Black Hole Domination_ , _JHEP_ 08 (2019) 001 [1905.01301].
  * [6] I. Masina, _Dark matter and dark radiation from evaporating primordial black holes_ , _Eur. Phys. J. Plus_ 135 (2020) 552 [2004.04740].
  * [7] I. Baldes, Q. Decant, D.C. Hooper and L. Lopez-Honorez, _Non-Cold Dark Matter from Primordial Black Hole Evaporation_ , _JCAP_ 08 (2020) 045 [2004.14773].
  * [8] P. Gondolo, P. Sandick and B. Shams Es Haghi, _Effects of primordial black holes on dark matter models_ , _Phys. Rev. D_ 102 (2020) 095018 [2009.02424].
  * [9] N. Bernal and Ó. Zapata, _Self-interacting Dark Matter from Primordial Black Holes_ , _JCAP_ 03 (2021) 007 [2010.09725].
  * [10] N. Bernal and Ó. Zapata, _Gravitational dark matter production: primordial black holes and UV freeze-in_ , _Phys. Lett. B_ 815 (2021) 136129 [2011.02510].
  * [11] N. Bernal and Ó. Zapata, _Dark Matter in the Time of Primordial Black Holes_ , _JCAP_ 03 (2021) 015 [2011.12306].
  * [12] I. Masina, _Dark matter and dark radiation from evaporating Kerr primordial black holes_ , 2103.13825.
  * [13] A. Arbey, J. Auffinger, P. Sandick, B. Shams Es Haghi and K. Sinha, _Precision calculation of dark radiation from spinning primordial black holes and early matter-dominated eras_ , _Phys. Rev. D_ 103 (2021) 123549 [2104.04051].
  * [14] S. Jyoti Das, D. Mahanta and D. Borah, _Low scale leptogenesis and dark matter in the presence of primordial black holes_ , 2104.14496.
  * [15] A. Cheek, L. Heurtier, Y.F. Perez-Gonzalez and J. Turner, _Primordial Black Hole Evaporation and Dark Matter Production: I. Solely Hawking radiation_ , 2107.00013.
  * [16] A. Cheek, L. Heurtier, Y.F. Perez-Gonzalez and J. Turner, _Primordial Black Hole Evaporation and Dark Matter Production: II. Interplay with the Freeze-In/Out Mechanism_ , 2107.00016.
  * [17] P. Sandick, B.S. Es Haghi and K. Sinha, _Asymmetric Reheating by Primordial Black Holes_ , 2108.08329.
  * [18] J.D. Barrow, E.J. Copeland, E.W. Kolb and A.R. Liddle, _Baryogenesis in extended inflation. 2. Baryogenesis via primordial black holes_ , _Phys. Rev. D_ 43 (1991) 984.
  * [19] Y. Hamada and S. Iso, _Baryon asymmetry from primordial black holes_ , _PTEP_ 2017 (2017) 033B02 [1610.02586].
  * [20] D. Hooper and G. Krnjaic, _GUT Baryogenesis With Primordial Black Holes_ , _Phys. Rev. D_ 103 (2021) 043504 [2010.01134].
  * [21] Y.F. Perez-Gonzalez and J. Turner, _Assessing the tension between a black hole dominated early universe and leptogenesis_ , 2010.03565.
  * [22] S. Datta, A. Ghosal and R. Samanta, _Baryogenesis from ultralight primordial black holes and strong gravitational waves from cosmic strings_ , _JCAP_ 08 (2021) 021 [2012.14981].
  * [23] C. Lunardini and Y.F. Perez-Gonzalez, _Dirac and Majorana neutrino signatures of primordial black holes_ , _JCAP_ 08 (2020) 014 [1910.07864].
  * [24] S. Weinberg, _A New Light Boson?_ , _Phys. Rev. Lett._ 40 (1978) 223.
  * [25] F. Wilczek, _Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons_, _Phys. Rev. Lett._ 40 (1978) 279.
  * [26] R.D. Peccei and H.R. Quinn, _CP Conservation in the Presence of Instantons_ , _Phys. Rev. Lett._ 38 (1977) 1440.
  * [27] R.D. Peccei and H.R. Quinn, _Some Aspects of Instantons_ , _Nuovo Cim. A_ 41 (1977) 309.
  * [28] R.D. Peccei and H.R. Quinn, _Constraints Imposed by CP Conservation in the Presence of Instantons_ , _Phys. Rev. D_ 16 (1977) 1791.
  * [29] J. Preskill, M.B. Wise and F. Wilczek, _Cosmology of the Invisible Axion_ , _Phys. Lett. B_ 120 (1983) 127.
  * [30] L.F. Abbott and P. Sikivie, _A Cosmological Bound on the Invisible Axion_ , _Phys. Lett. B_ 120 (1983) 133.
  * [31] M. Dine and W. Fischler, _The Not So Harmless Axion_ , _Phys. Lett. B_ 120 (1983) 137.
  * [32] P. Sikivie, _Axion Cosmology_ , _Lect. Notes Phys._ 741 (2008) 19 [astro-ph/0610440].
  * [33] D.J.E. Marsh, _Axion Cosmology_ , _Phys. Rept._ 643 (2016) 1 [1510.07633].
  * [34] L. Di Luzio, M. Giannotti, E. Nardi and L. Visinelli, _The landscape of QCD axion models_ , _Phys. Rept._ 870 (2020) 1 [2003.01100].
  * [35] P. Sikivie, _Invisible Axion Search Methods_ , _Rev. Mod. Phys._ 93 (2021) 015004 [2003.02206].
  * [36] S. Bird, I. Cholis, J.B. Muñoz, Y. Ali-Haïmoud, M. Kamionkowski, E.D. Kovetz et al., _Did LIGO detect dark matter?_ , _Phys. Rev. Lett._ 116 (2016) 201301 [1603.00464].
  * [37] M. Sasaki, T. Suyama, T. Tanaka and S. Yokoyama, _Primordial Black Hole Scenario for the Gravitational-Wave Event GW150914_ , _Phys. Rev. Lett._ 117 (2016) 061101 [1603.08338].
  * [38] S. Clesse and J. García-Bellido, _The clustering of massive Primordial Black Holes as Dark Matter: measuring their mass distribution with Advanced LIGO_ , _Phys. Dark Univ._ 15 (2017) 142 [1603.05234].
  * [39] S. Clesse and J. García-Bellido, _Seven Hints for Primordial Black Hole Dark Matter_ , _Phys. Dark Univ._ 22 (2018) 137 [1711.10458].
  * [40] A.D. Gow, C.T. Byrnes, A. Hall and J.A. Peacock, _Primordial black hole merger rates: distributions for multiple LIGO observables_ , _JCAP_ 01 (2020) 031 [1911.12685].
  * [41] F. Schiavone, D. Montanino, A. Mirizzi and F. Capozzi, _Axion-like particles from primordial black holes shining through the Universe_ , _JCAP_ 08 (2021) 063 [2107.03420].
  * [42] N. Bernal, F. Hajkarim and Y. Xu, _Axion Dark Matter in the Time of Primordial Black Holes_ , _Phys. Rev. D_ 104 (2021) 075007 [2107.13575].
  * [43] J. Jaeckel and A. Ringwald, _The Low-Energy Frontier of Particle Physics_ , _Ann. Rev. Nucl. Part. Sci._ 60 (2010) 405 [1002.0329].
  * [44] P. Arias, D. Cadamuro, M. Goodsell, J. Jaeckel, J. Redondo and A. Ringwald, _WISPy Cold Dark Matter_ , _JCAP_ 06 (2012) 013 [1201.5902].
  * [45] A. Ringwald, _Exploring the Role of Axions and Other WISPs in the Dark Universe_ , _Phys. Dark Univ._ 1 (2012) 116 [1210.5081].
  * [46] I.G. Irastorza and J. Redondo, _New experimental approaches in the search for axion-like particles_ , _Prog. Part. Nucl. Phys._ 102 (2018) 89 [1801.08127].
  * [47] P. Svrcek and E. Witten, _Axions In String Theory_ , _JHEP_ 06 (2006) 051 [hep-th/0605206].
  * [48] A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper and J. March-Russell, _String Axiverse_ , _Phys. Rev. D_ 81 (2010) 123530 [0905.4720].
  * [49] N. Blinov, M.J. Dolan, P. Draper and J. Kozaczuk, _Dark matter targets for axionlike particle searches_ , _Phys. Rev. D_ 100 (2019) 015049 [1905.06952].
  * [50] T. Papanikolaou, V. Vennin and D. Langlois, _Gravitational waves from a universe filled with primordial black holes_ , _JCAP_ 03 (2021) 053 [2010.11573].
  * [51] G. Domènech, C. Lin and M. Sasaki, _Gravitational wave constraints on the primordial black hole dominated early universe_ , _JCAP_ 04 (2021) 062 [2012.08151].
  * [52] G. Domènech, V. Takhistov and M. Sasaki, _Exploring Evaporating Primordial Black Holes with Gravitational Waves_ , 2105.06816.
  * [53] G. Domènech, _Scalar induced gravitational waves review_ , 2109.01398.
  * [54] C.J. Hogan and M.J. Rees, _Axion Miniclusters_ , _Phys. Lett. B_ 205 (1988) 228.
  * [55] E.W. Kolb and I.I. Tkachev, _Nonlinear axion dynamics and formation of cosmological pseudosolitons_ , _Phys. Rev. D_ 49 (1994) 5040 [astro-ph/9311037].
  * [56] E.W. Kolb and I.I. Tkachev, _Axion miniclusters and Bose stars_ , _Phys. Rev. Lett._ 71 (1993) 3051 [hep-ph/9303313].
  * [57] E.W. Kolb and I.I. Tkachev, _Large amplitude isothermal fluctuations and high density dark matter clumps_ , _Phys. Rev. D_ 50 (1994) 769 [astro-ph/9403011].
  * [58] B. Eggemeier, J. Redondo, K. Dolag, J.C. Niemeyer and A. Vaquero, _First Simulations of Axion Minicluster Halos_ , _Phys. Rev. Lett._ 125 (2020) 041301 [1911.09417].
  * [59] CMB-S4 collaboration, _CMB-S4 Science Book, First Edition_ , 1610.02743.
  * [60] NASA PICO collaboration, _PICO: Probe of Inflation and Cosmic Origins_ , 1902.10541.
  * [61] SPT-3G collaboration, _SPT-3G: A Next-Generation Cosmic Microwave Background Polarization Experiment on the South Pole Telescope_ , _Proc. SPIE Int. Soc. Opt. Eng._ 9153 (2014) 91531P [1407.2973].
  * [62] Simons Observatory collaboration, _The Simons Observatory: Science goals and forecasts_ , _JCAP_ 1902 (2019) 056 [1808.07445].
  * [63] S.W. Hawking, _Particle Creation by Black Holes_ , _Commun. Math. Phys._ 43 (1975) 199.
  * [64] B.J. Carr and S.W. Hawking, _Black holes in the early Universe_ , _Mon. Not. Roy. Astron. Soc._ 168 (1974) 399.
  * [65] W.H. Press and P. Schechter, _Formation of galaxies and clusters of galaxies by selfsimilar gravitational condensation_ , _Astrophys. J._ 187 (1974) 425.
  * [66] B.J. Carr, K. Kohri, Y. Sendouda and J. Yokoyama, _New cosmological constraints on primordial black holes_ , _Phys. Rev. D_ 81 (2010) 104019 [0912.5297].
  * [67] D.N. Page, _Particle Emission Rates from a Black Hole: Massless Particles from an Uncharged, Nonrotating Hole_ , _Phys. Rev. D_ 13 (1976) 198.
  * [68] J.H. MacGibbon and B.R. Webber, _Quark and gluon jet emission from primordial black holes: The instantaneous spectra_ , _Phys. Rev. D_ 41 (1990) 3052.
  * [69] J.H. MacGibbon, _Quark and gluon jet emission from primordial black holes. 2. The Lifetime emission_ , _Phys. Rev. D_ 44 (1991) 376.
  * [70] R. Allahverdi et al., _The First Three Seconds: a Review of Possible Expansion Histories of the Early Universe_ , _Open J. Astrophys._ 4 (2020)  [2006.16182].
  * [71] F.W. Stecker and Q. Shafi, _The Evolution of Structure in the Universe From Axions_ , _Phys. Rev. Lett._ 50 (1983) 928.
  * [72] Planck collaboration, _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6 [1807.06209].
  * [73] S. Sarkar, _Big bang nucleosynthesis and physics beyond the standard model_ , _Rept. Prog. Phys._ 59 (1996) 1493 [hep-ph/9602260].
  * [74] S. Hannestad, _What is the lowest possible reheating temperature?_ , _Phys. Rev. D_ 70 (2004) 043506 [astro-ph/0403291].
  * [75] P.F. de Salas, M. Lattanzi, G. Mangano, G. Miele, S. Pastor and O. Pisanti, _Bounds on very low reheating scenarios after Planck_ , _Phys. Rev. D_ 92 (2015) 123534 [1511.00672].
  * [76] J.E. Kim, _Weak Interaction Singlet and Strong CP Invariance_ , _Phys. Rev. Lett._ 43 (1979) 103.
  * [77] M.A. Shifman, A.I. Vainshtein and V.I. Zakharov, _Can Confinement Ensure Natural CP Invariance of Strong Interactions?_ , _Nucl. Phys. B_ 166 (1980) 493.
  * [78] A.R. Zhitnitsky, _On Possible Suppression of the Axion Hadron Interactions. (In Russian)_ , _Sov. J. Nucl. Phys._ 31 (1980) 260.
  * [79] M. Dine, W. Fischler and M. Srednicki, _A Simple Solution to the Strong CP Problem with a Harmless Axion_ , _Phys. Lett. B_ 104 (1981) 199.
  * [80] P.W. Graham, I.G. Irastorza, S.K. Lamoreaux, A. Lindner and K.A. van Bibber, _Experimental Searches for the Axion and Axion-Like Particles_ , _Ann. Rev. Nucl. Part. Sci._ 65 (2015) 485 [1602.00039].
  * [81] R.L. Davis, _Goldstone Bosons in String Models of Galaxy Formation_ , _Phys. Rev. D_ 32 (1985) 3172.
  * [82] D. Harari and P. Sikivie, _On the Evolution of Global Strings in the Early Universe_ , _Phys. Lett. B_ 195 (1987) 361.
  * [83] R.A. Battye and E.P.S. Shellard, _Axion string constraints_ , _Phys. Rev. Lett._ 73 (1994) 2954 [astro-ph/9403018].
  * [84] T. Hiramatsu, M. Kawasaki, K. Saikawa and T. Sekiguchi, _Production of dark matter axions from collapse of string-wall systems_ , _Phys. Rev. D_ 85 (2012) 105020 [1202.5851].
  * [85] C. O’Hare, _cajohare/axionlimits: Axionlimits_ , July, 2020.  10.5281/zenodo.3932430.
  * [86] V.M. Mehta, M. Demirtas, C. Long, D.J.E. Marsh, L. Mcallister and M.J. Stott, _Superradiance Exclusions in the Landscape of Type IIB String Theory_ , 2011.08693.
  * [87] D. Cadamuro and J. Redondo, _Cosmological bounds on pseudo Nambu-Goldstone bosons_ , _JCAP_ 02 (2012) 032 [1110.2895].
  * [88] P. Sikivie, _Experimental Tests of the Invisible Axion_ , _Phys. Rev. Lett._ 51 (1983) 1415.
  * [89] ADMX collaboration, _A SQUID-based microwave cavity search for dark-matter axions_ , _Phys. Rev. Lett._ 104 (2010) 041301 [0910.5914].
  * [90] ADMX collaboration, _Extended Search for the Invisible Axion with the Axion Dark Matter Experiment_ , _Phys. Rev. Lett._ 124 (2020) 101303 [1910.08638].
  * [91] HAYSTAC collaboration, _Results from phase 1 of the HAYSTAC microwave cavity axion experiment_ , _Phys. Rev. D_ 97 (2018) 092001 [1803.03690].
  * [92] S. Lee, S. Ahn, J. Choi, B.R. Ko and Y.K. Semertzidis, _Axion Dark Matter Search around 6.7 $\mu$eV_, _Phys. Rev. Lett._ 124 (2020) 101802 [2001.05102].
  * [93] J. Jeong, S. Youn, S. Bae, J. Kim, T. Seong, J.E. Kim et al., _Search for Invisible Axion Dark Matter with a Multiple-Cell Haloscope_ , _Phys. Rev. Lett._ 125 (2020) 221302 [2008.10141].
  * [94] CAPP collaboration, _First Results from an Axion Haloscope at CAPP around 10.7 $\mu$eV_, _Phys. Rev. Lett._ 126 (2021) 191802 [2012.10764].
  * [95] D. Alesini et al., _Galactic axions search with a superconducting resonant cavity_ , _Phys. Rev. D_ 99 (2019) 101101 [1903.06547].
  * [96] D. Alesini et al., _Search for invisible axion dark matter of mass m ${}_{a}=43~{}\mu$eV with the QUAX–$a\gamma$ experiment_, _Phys. Rev. D_ 103 (2021) 102004 [2012.09498].
  * [97] B.T. McAllister, G. Flower, E.N. Ivanov, M. Goryachev, J. Bourhill and M.E. Tobar, _The ORGAN Experiment: An axion haloscope above 15 GHz_ , _Phys. Dark Univ._ 18 (2017) 67 [1706.00209].
  * [98] S. Beurthey et al., _MADMAX Status Report_ , 2003.10894.
  * [99] M. Lawson, A.J. Millar, M. Pancaldi, E. Vitagliano and F. Wilczek, _Tunable axion plasma haloscopes_ , _Phys. Rev. Lett._ 123 (2019) 141802 [1904.11872].
  * [100] I. Stern, _ADMX Status_ , _PoS_ ICHEP2016 (2016) 198 [1612.08296].
  * [101] D. Alesini, D. Babusci, D. Di Gioacchino, C. Gatti, G. Lamanna and C. Ligi, _The KLASH Proposal_ , 1707.06010.
  * [102] S. Chaudhuri, K. Irwin, P.W. Graham and J. Mardon, _Optimal Impedance Matching and Quantum Limits of Electromagnetic Axion and Hidden-Photon Dark Matter Searches_ , 1803.01627.
  * [103] Y. Kahn, B.R. Safdi and J. Thaler, _Broadband and Resonant Approaches to Axion Dark Matter Detection_ , _Phys. Rev. Lett._ 117 (2016) 141801 [1602.01086].
  * [104] 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) 121802 [1810.12257].
  * [105] Planck collaboration, _Planck 2018 results. X. Constraints on inflation_ , _Astron. Astrophys._ 641 (2020) A10 [1807.06211].
  * [106] G. Grilli di Cortona, E. Hardy, J. Pardo Vega and G. Villadoro, _The QCD axion, precisely_ , _JHEP_ 01 (2016) 034 [1511.02867].
  * [107] E. Hardy, _Miniclusters in the Axiverse_ , _JHEP_ 02 (2017) 046 [1609.00208].
  * [108] A.E. Nelson and H. Xiao, _Axion Cosmology with Early Matter Domination_ , _Phys. Rev. D_ 98 (2018) 063516 [1807.07176].
  * [109] L. Visinelli and J. Redondo, _Axion Miniclusters in Modified Cosmological Histories_ , _Phys. Rev. D_ 101 (2020) 023008 [1808.01879].
  * [110] N. Blinov, M.J. Dolan and P. Draper, _Imprints of the Early Universe on Axion Dark Matter Substructure_ , _Phys. Rev. D_ 101 (2020) 035002 [1911.07853].
  * [111] E.W. Kolb and I.I. Tkachev, _Femtolensing and picolensing by axion miniclusters_ , _Astrophys. J. Lett._ 460 (1996) L25 [astro-ph/9510043].
  * [112] K.M. Zurek, C.J. Hogan and T.R. Quinn, _Astrophysical Effects of Scalar Dark Matter Miniclusters_ , _Phys. Rev. D_ 75 (2007) 043511 [astro-ph/0607341].
  * [113] M. Fairbairn, D.J.E. Marsh and J. Quevillon, _Searching for the QCD Axion with Gravitational Microlensing_ , _Phys. Rev. Lett._ 119 (2017) 021101 [1701.04787].
  * [114] M. Fairbairn, D.J.E. Marsh, J. Quevillon and S. Rozier, _Structure formation and microlensing with axion miniclusters_ , _Phys. Rev. D_ 97 (2018) 083502 [1707.03310].
  * [115] S. Borsanyi et al., _Calculation of the axion mass based on high-temperature lattice quantum chromodynamics_ , _Nature_ 539 (2016) 69 [1606.07494].
  * [116] M. Venegas, _Relic Density of Axion Dark Matter in Standard and Non-Standard Cosmological Scenarios_ , 2106.07796.
  * [117] P. Arias, N. Bernal, D. Karamitros, C. Maldonado, L. Roszkowski and M. Venegas, _New opportunities for axion dark matter searches in nonstandard cosmological models_ , 2107.13588.
  * [118] D. Hooper, G. Krnjaic, J. March-Russell, S.D. McDermott and R. Petrossian-Byrne, _Hot Gravitons and Gravitational Waves From Kerr Black Holes in the Early Universe_ , 2004.00618.
  * [119] A.G. Riess et al., _A 2.4% Determination of the Local Value of the Hubble Constant_ , _Astrophys. J._ 826 (2016) 56 [1604.01424].
  * [120] M. Escudero, D. Hooper, G. Krnjaic and M. Pierre, _Cosmology with A Very Light L μ $-$ Lτ Gauge Boson_, _JHEP_ 03 (2019) 071 [1901.02010].
  * [121] S. Vagnozzi, _New physics in light of the $H_{0}$ tension: An alternative view_, _Phys. Rev. D_ 102 (2020) 023518 [1907.07569].
  * [122] J. Alcaniz, N. Bernal, A. Masiero and F.S. Queiroz, _Light dark matter: A common solution to the lithium and $H_{0}$ problems_, _Phys. Lett. B_ 812 (2021) 136008 [1912.05563].