ALP Dark Matter in a Primordial Black Hole Dominated Universe

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.

The QCD axion [25,26], a pseudo-Nambu-Goldstone boson resulting from the spontaneous breaking of the global U (1) Peccei-Quinn symmetry [27][28][29], 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 [30][31][32] and the decay of topological defects [33] 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 [10 −6 − 10 −5 ] eV (or correspondingly to a Peccei-Quinn scale f a 10 12 GeV) is expected, if one does not want to introduce fine tuning of the initial misalignment angle θ i [34][35][36] 2 . In light of the accumulated interesting implications for PBHs [39][40][41][42][43][44], it has been recently shown that the axion DM window can be enlarged to a mass as low as m a ∼ O(10 −8 ) eV if the Universe features an early PBHs dominated epoch [45].
In this paper, we go beyond the QCD axion DM and consider general ALPs, see e.g. Refs. [34,35,[46][47][48][49]. 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 [50,51]. 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. [48]. In this paper, we focus on ALP DM generated via the usual misalignment mechanism. Therefore, our results are complementary to those presented in Ref. [44]. 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 θ i , leading to strong bounds on the viable parameter space for natural values of θ 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. [52] recently, where it was shown that a broader parameter is available and could be within the reach of several proposed experiments. Similarly to Ref. [45], and differently from Ref. [52], 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 [53][54][55][56][57], 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γ 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 [58][59][60][61][62]. 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 ∆N eff , which is within the sensitivity of next generation CMB experiments [63][64][65][66]. We also find that for Kerr PBHs, ∆N 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 = c = k B = 1 throughout this manuscript.

Primordial Black Hole Evaporation
PBHs are hypothetical objects which could be generated due to inhomogeneities of density perturbations in the early Universe [67,68]. When these fluctuations reenter the horizon, if above a threshold, they could collapse and form a BH according to the Press-Schechter formalism [69]. In this paper, we focus on the case where PBHs form in a radiation-dominated epoch. With an initial cosmic temperature T = T in , the initial PBH mass is given by the whole mass within the particle horizon [1,70], where we take the efficiency factor to be γ 0.2. Additionally, ρ 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 ρ BH (T in ) is usually related to the SM radiation energy density at formation via the β parameter Since ρ BH redshifts slower than radiation, an early PBH-dominated epoch with where T ev , given by [45] T ev g (T in ) 640 is the SM temperature at which PBHs completely evaporate. Additionally, M P ≡ 1/ √ 8π G denotes the reduced Planck mass, and g (T ) corresponds to the number of relativistic degrees of freedom contributing to the SM energy density. It is interesting to note that if β > β c , PBHs start to dominate the total energy density of the Universe at T = T eq , defined as ρ R (T eq ) = ρ BH (T eq ), with with g 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 β and mass M in [1,2]. Here, we will focus on PBHs that evaporated before BBN, thus having masses smaller than ∼ 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 [54]. More importantly, a modification of BBN predictions due to the energy density stored in GWs can be avoided if [55]   .
(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. [15,16,71].

Schwarzschild Black Holes
For the most simple BH scenario, PBHs are only characterized by their mass, so that their horizon temperature T BH is determined according to [67] (2.7) All particles with masses smaller than T 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 µ i and internal degrees of freedom g i can be described as [15,67] , with E i (p) = p 2 + µ 2 i denoting the energy and σ 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 [72,73] with the mass evaporation function ε (M BH ) ≡ i g i ε i (z i ), where the contribution per degree of freedom ε i (z i ) is given by [15]

Kerr Black Holes
For the case of BHs with a nonzero spin, the horizon temperature has an additional dependence on the spin parameter a ∈ [0, 1] (with a = 0 corresponding to the Schwarzschild limit) where a ≡ 8π J M 2 P /M 2 BH , 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, with Ω = 4π a M 2 P /M BH 1 + 1 − a 2 −1 horizon angular velocity. The time evolution for the BH spin and mass is described by the following set of equations [71] where the evaporation functions ε (M BH , a ) and γ (M BH ) are given by . 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.

Background Evolution
The evolution of the SM entropy density s(T ) = 2π 2 45 g s (T ) T 3 can be tracked via the Boltzmann equation 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 However, PBHs can have a strong impact on the evolution of the background dynamics [74].
With an early PBH dominated epoch, cosmology can be characterized by four distinct regimes, where the Hubble expansion rate is given by [45] H(T )

3)
3 Note that the definition of x is the same as in the Schwarzschild case.
which can be simplified to for being more conveniently used in next expressions. Here, we have introduced the temperature scale T c , from which the SM radiation does not scale as free radiation (i.e., ρ R (R) ∝ R −1 with R being the scale factor) due to the entropy injected by the PBH evaporation. It is given by [45] T c g (T in ) π 5760 Additionally, during its evaporation, PBHs radiate SM particles and therefore dilute all previously produced species. The entropy injection factor is [45] for T ev ≥ T , (3.6) which can be simplified to 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 = 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 = 0 whenever we talk about Schwarzschild PBHs. Besides, in the evaporation functions ε (M BH , a ) and γ (M BH ) we have included an additional, almost massless, pseudoscalar degree of freedom, corresponding to the ALP or the axion.

ALP Dark Matter
In this section, we focus on a general light pseudoscalar, namely the ALP [46][47][48]. 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 [50,51].
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 [30][31][32]75] where θ(t) ≡ 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 osc defined as In the standard cosmological scenario, as the SM entropy is conserved, the ALP energy density ρ a at present is given by where T 0 is the CMB temperature at present, and we considered the conservation of the ALP number density in a comoving volume. Additionally, ρ a (T osc ) 1 2 m 2 a f 2 a θ 2 i in the limit in which the kinetic energy is neglected, and for a quadratic potential. Moreover, θ i denotes the initial misalignment angle. The entropy injection from PBHs evaporation would inevitably dilute the energy density so that one has 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 where ρ c /h 2 1.1×10 −5 GeV/cm 3 is the critical energy density and s(T 0 ) 2.69×10 3 cm −3 is the entropy density at present, which should match the total DM abundance Ωh 2 0.12 [76].
Knowing the Hubble expansion rate given in the previous section, the ALP oscillation temperature T osc can be estimated to be Using Eq. (4.4), one can compute the relic density in the aforementioned four regimes, which becomes or equivalently Note that the expressions for Ω a h 2 are identical for T eq ≥ T osc ≥ T c and T c ≥ T 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 = 0, red band) or Kerr (a = 0.999, blue band) PBHs. The thickness of the bands corresponds to an initial misalignment angle 0.5 ≤ θ i ≤ π/ √ 3. 4 The regions on the right (left) of the bands produce a DM underabundance (overabundance). Additionally, the gray regions are in tension with BBN (T ev < 4 MeV [77][78][79][80][81][82]), whereas the green region with GWs (i.e., Eq. (2.6)). Finally, the dashed red line corresponds to β = β 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. [15,16], including the additional ALP degree of freedom. We have checked that the numerical estimations in Eq. (4.7) fit well the analytical results.
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 β allowed by GWs. This happens for M in 5.7 × 10 8 g and β 4.6 × 10 −10 for a = 0, or M in 7.6 × 10 8 g and β 3.5 × 10 −10 for a = 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 in 5.7 × 10 8 g and β 4.6 × 10 −10 , for a = 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 eq 81 GeV and T = T ev 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 osc decreases for a PBH dominated era. Again, these results were obtained using the numerical code and fit well the analytical estimations.
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 where the coupling constant g aγ is model dependent and related to the breaking scale of the PQ symmetry as with z ≡ 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 [83,84], whereas for DFSZ models E/N = 8/3 [85,86]. The electromagnetic interaction of ALPs is by far the most exploited to look for signatures in observations and experimental searches [87].
Once one assumes that all DM is composed by ALPs, m a and f a (and therefore g aγ ) 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γ |). 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 θ i in the range [0.5, π/ √ 3], with θ i = π/ √ 3 being the top border. The regions on the top and bottom give rise to a DM underabundance and overabundance, respectively. 5 This parameter space is modified once an early PBH dominated epoch occurs. The red (blue) thick band shows the viable parameter space for a = 0 (a = 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γ |. It is interesting to note that for m a 10 −5 eV, T osc < T eq , and therefore the ALP relic density is independent on its mass. However, for m a 10 −5 eV, T osc > T eq and g aγ ∝ m 1/4 a , cf. Eq. (4.7). Figure 3 also overlays in brown current bounds on the ALP-photon coupling. Small masses (m a 4 × 10 −11 eV) are constrained astrophysical black hole spins [94], whereas high masses (m a 10 2 eV) by the X-ray background and extragalactic background light [95]. The masses in the range 10 −6 eV m a 10 −5 eV are already probed by the so-called haloscope experiments [96], ADMX [97,98], HAYSTAC [99], CAPP [100][101][102], QUAX [103,104], ORGAN [105], using a highly tuned microwave cavity that converts ALPs into photons in the presence of a static magnetic field. 6 Additionally, the light blue regions show projected sensitivities from a number of experiments that count ORGAN [105], MADMAX [106], AL-PHA [107], ADMX [108], KLASH [109], DM-Radio [110] and ABRACADABRA [111,112]. These bounds have been adapted from Ref. [93]. Even if the PBH domination tends to imply smaller values for |g aγ |, this parameter space could be potentially tested in the future by next-generation experiments, especially for light ALP masses, m a 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 [32]. The ALP field is homogeneous through various Hubble patches, with a unique value of θ 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 [113] impose the lower bound on f a which can be cast as [35] H I 0.9 × 10 7 GeV θ i π f a 10 11 GeV . (4.11) Since the BBN scale represents a lower bound for H I , it follows that f a 100 GeV (|g aγ | 10 −5 GeV −1 ) for θ i ∼ 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 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 in 6 × 10 8 g, and therefore for θ i ∼ O(1), it implies that |g aγ | 10 −12 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 [−π, π)) so that θ i = θ 2 i = π/ √ 3. 7 Additionally, the upper limit on the inflationary scale H I < 2.5 × 10 −5 M P [113] implies that f a 10 13 GeV (and in turn g aγ 10 −16 GeV −1 ), which finally can be translated into an lower bound on the ALP mass m a 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 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 [35,[58][59][60][61].
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 osc ∼ H(T osc ) −1 , at the time of the matter-radiation equality. Therefore, their mass at formation is (4.14) In a radiation dominated Universe and assuming that all relic ALPs bound up in miniclusters, it reduces to 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 [95]. However, if ALPs start to oscillate during the PBH dominated period, the initial minicluster mass can be estimated to be 8 Schwarzschild, or M in 7.6 × 10 8 g and β 3.5 × 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. [116][117][118] for related studies). Hence, searches for ALPs, such as indirect detection through gravitational microlensing [119][120][121][122], could be extended to regions of the parameter space previously thought to be not relevant.

Axion Dark Matter
The production of QCD axion in a cosmology dominated by PBHs was analytically studied in Ref. [45]. 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 [123] with T QCD 150 MeV. As the axion mass is temperature dependent, its abundance in the standard cosmology is 4.8 × 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 [45,124,125]. 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. [45]. 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 = 0) and Kerr (a = 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 ev > 4 MeV) and GWs (cf. Eq. (2.6)). Once we restrict to initial misalignment angles in the range θ i ∈ [0.5, π/ √ 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 ∼ [10 −6 , 10 −5 ] eV to ∼ [10 −8 , 10 −5 ] eV. 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 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 Contrary to the ALP case, there is no significant lower bound on f a when M 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 10 13 GeV, and hence the lower bound on the axion mass m a 6 × 10 −7 eV. In the standard cosmological scenario, axions start oscillating at the temperature T osc 1.1 GeV. Therefore, if the misalignment mechanism is the unique method to produce axions, they should have a mass m a 1.8 × 10 −5 eV (for θ i = π/ √ 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 osc )/S(T 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. 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 radiationmatter equality leads to miniclusters having a mass M 0 2.1 × 10 −11 M 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 T ev for T eq ≥ T osc ≥ T c , π 1/2 10 g (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 .

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,35,44]. We note that ALP DR directly emitted from nonrotating PBHs was analyzed in Ref. [44] 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 where T equality 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 with ε a (T ev ) being the axion evaporation function, and ε R (T ev ) the total SM radiation contribution to the evaporation function. For the case of Schwarzschild PBHs, ∆N eff is simply which is independent from β, and consistent with the result obtained in Ref. [44]. For Kerr PBHs, the time depletion of the angular momentum makes it more difficult to obtain a simple analytical form. Left (right) panel of Fig. 6 shows the contribution to ∆N eff as function of BH mass and β for Schwarzschild (Kerr) PBHs. Let us note that in a PBH-dominated era, ∆N eff tends to be independent of β, as noted before, c.f. Eq. (6.3). On the contrary, when β < β c , larger M in is needed with smaller β in order to yield a fixed ∆N eff . We also note that when M in 2 × 10 7 g (correspondingly T ev 0.1 GeV), there are some features in the lines shown in Fig. 6. This is because the sharp change of g and g s due to QCD phase transition [5,23]. Comparing the resulting values of ∆N 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 [71]. 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 ∆N eff is in fact reduced [12,126]. The current upper bound on ∆N eff comes from the Planck collaboration, and reads N eff = 2.99 ± 0.17 [76]. However, ∆N eff due to emitted hot axions/ALPs could reach the sensitivity of next generation CMB experiments: CMB-S4 [63] (∆N eff ∼ 0.06), PICO [64] (∆N eff ∼ 0.06) and SPT-3G/SO [65,66] (∆N eff ∼ 0.1), and therefore axion/ALP DR could be tested in the near future. 9

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 [10 −6 −10 −5 ] eV gets broaden to ∼ [10 −8 − 10 −5 ] eV (for misalignment angles of order 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 ∆N 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.