Primordial black hole superradiance and evaporation in the string axiverse

In the string axiverse scenario, light primordial black holes may spin up due to the Hawking emission of a large number of light (sub-MeV) axions. We show that this may trigger superradiant instabilities associated with a heavier axion during the black holes' evolution, and study the coupled dynamics of superradiance and evaporation. We find, in particular, that the present black hole mass-spin distribution should follow the superradiance threshold condition for black hole masses below the value at which the superradiant cloud forms, for a given heavy axion mass. Furthermore, we show that the decay of the heavy axions within the superradiant cloud into photon pairs may lead to a distinctive line in the black hole's emission spectrum, superimposed on its electromagnetic Hawking emission.


Introduction
There has been a growing interest in the literature on the study of primordial black holes (PBHs), originally predicted by Hawking to form in the early Universe through the direct gravitational collapse of overdense regions [1][2][3].These PBHs are natural candidates to account for at least a fraction, and potentially all, the dark matter in the Universe [4,5].Within the standard cosmological paradigm, such PBHs would be born with very little natal spin, given that the ambient radiation pressure would lead to a nearly spherical gravitational collapse.This could explain why the merging BHs recently detected with the LIGO/Virgo/Kagra gravitational wave interferometers seem to be slowly rotating [6].
PBHs, particularly light ones with sub-solar masses, also offer new avenues for research in fundamental particle physics, given their ability to produce large numbers of particles, both known and exotic, through Hawking emission [7] and rotational superradiance (see [8] and references therein).Although in both cases the total number of particles produced over the PBHs lifetime decreases with the PBH mass as (M/M P ) 2 , where M P ≃ 2.176 × 10 −8 kg denotes the Planck mass, light PBHs have larger Hawking temperatures T H ≃ M 2 P /8πM and rotate with larger angular velocities Ω H ≃ ã/2M (1 + √ 1 − ã2 ), for a given dimensionless spin parameter ã.Since Hawking emission and rotational superradiance can only efficiently produce particles of mass µ ≲ T H and µ ≲ Ω H , light PBHs may therefore generate much heavier particles than their stellar or supermassive counterparts.This is appealing since most extensions of the Standard Model predict the existence of exotic particles across many orders of magnitude in mass, well beyond the "ultra-light regime" accessible with known astrophysical BHs .
The timescales involved in Hawking evaporation and superradiance also decrease with the BH mass, so that particle production may efficiently occur at early times or even continuously throughout the cosmic history.PBHs have, in particular, been shown to be relevant for the generation of dark matter [36][37][38][39][40][41][42], including both light (but not ultralight) and heavy axions [43][44][45][46], and to probe the existence of new particles beyond the energy/luminosity reach of current particle accelerators [47,48].If they have indeed formed in the early Universe, PBHs may thus provide unique laboratories for fundamental physics.
A particularly interesting setup for PBH particle production is the string axiverse scenario, which conjectures that realistic string theory compactifications lead to hundreds or even thousands of light axion fields, whose masses are generated only through nonperturbative effects.The number of axion fields is simply dictated by the large number of non-trivial cycles in the six compact extra-dimensions supporting the higher-dimensional Neveu-Schwarz and Ramond-Ramond form-fields, each cycle yielding a pseudo-scalar axion field endowed with a perturbative shift symmetry in the effective 4-dimensional field theory.While axion masses could result from supersymmetry breaking or the mechanism(s) responsible for moduli stabilization, the authors of [9] argued that in string compactifications realizing the Peccei-Quinn solution to the strong CP problem, and which therefore include at least one light QCD axion, all axion masses result exclusively from non-perturbative effects.Such scenarios will therefore include a large number of axions spanning a broad range of mass scales.
A population of PBHs born with masses ∼ 10 12 kg will evaporate within the Universe's lifetime.Within the string axiverse scenario, as two of us have shown with March-Russell in [46], such PBHs will emit not only photons, electrons and other Standard Model degrees of freedom but also all axions with masses below a few MeV.This changes not only their lifetime (and hence the initial mass of the PBHs that evaporate away before the present day) but also their spin.While the emission of fermions, vector bosons and gravitons inevitably carries away a BH's angular momentum, scalar particles like axions are the only ones that can be emitted in the l = 0 mode, as originally shown by Taylor, Hiscock and Chambers [49,50].Scalar particle emission therefore decreases the BH mass but not its angular momentum, making it spin faster.In [46] we have simulated the evaporation of light PBHs emitting all known Standard Model particles and an arbitrary number of light axions (mass ≲ MeV) and, as we review below, we concluded that slowly-rotating PBHs may develop spin parameters ã ≳ 0.1 before evaporating away in the presence of a few hundred light axions.The present mass-spin distribution of light PBHs (at different stages of the evaporation process) depends on the total number of light axions, thus providing a unique probe of the string axiverse.This has the appeal of being a purely gravitational probe of this scenario, independent of how the individual axions interact with known particles, as well as of the details of the axion mass spectrum.
The fact that PBHs naturally develop non-negligible spin parameters through Hawking emission in this scenario motivates exploring whether this may trigger superradiant instabilities.In particular, in the string axiverse spectrum there may also exist a number of heavy axions (µ ≳ MeV), since the non-perturbative nature of the axion mass generation mechanism only implies that their masses are exponentially suppressed compared to a high mass scale such as the supersymmetry breaking scale.Hence, if the PBHs are born with low spin, the condition for superradiant particle production µ < Ω H will only be satisfied once the PBHs evaporate sufficiently and spin up due to the emission of light axions.
In this work, we thus study the dynamical generation of superradiant heavy axion clouds around PBHs born with mass ∼ 10 12 kg throughout the cosmic history, including both superradiance and Hawking emission.We will show that indeed such clouds may form with two important observational consequences.
First, the formation of superradiant clouds spins down the PBHs faster than evaporation can spin them up.This modifies the present PBH mass-spin distribution such that the lightest PBHs (which have evaporated sufficiently for superradiant clouds to form) saturate the superradiance condition, Ω H ≃ µ, while the spins of the heavier PBHs are determined solely by their evaporation stage and, hence, by the number of light axion species.
Second, the decay of the heavy axions into photon pairs leads to a characteristic gamma-ray line in the PBH-axion cloud photon emission spectrum.This is a very unique signature since the Hawking emission spectrum (including both primary and secondary photons) evolves as the PBH evaporates, while the line has a fixed energy corresponding to approximately half of the heavy axion's mass.
We begin our discussion by reviewing PBH evaporation through Hawking emission in the next section, discussing in particular the string axiverse case.In section 3 we start by reviewing the dynamics of superradiant instabilities for massive scalar fields and then bring together these two particle production mechanisms by fist considering a toy model where a BH evaporates by emitting a single light axion and a superradiant instability is induced by another heavy axion.Although unrealistic, this toy model allows one to understand the basic dynamics of the problem towards exploring a more realistic setup where all Standard Model particles are included in the PBH evaporation process alongside an arbitrary number of light axions.In section 4, we compute the PBH photon emission spectrum including primary and secondary Hawking emission as well as heavy axion decay within the superradiant clouds.We summarize our main results and conclusions in the final section.
We note in advance that axion self-interactions are assumed to play a negligible role in the dynamical evolution.This is a good approximation for the large axion decay constants typically predicted in string constructions and that we also take into account when discussing observational prospects in section 4.

Hawking emission and black hole evaporation
In curved space-time, different observers do not necessarily agree in their definition of what is the quantum vacuum state, i.e. the state with the lowest possible energy.This is due to the use of different time coordinates to perform the separation between the positive and negative frequency modes that underlies the field quantization procedure.While in flat Minkowski space all inertial observers perform this mode separation in an equivalent way due to Lorentz invariance, in space-time manifolds that include regions with non-negligible curvature, particularly event horizons, this is typically not the case.
In particular, Hawking showed in 1974 [51,52] that a stationary (and therefore noninertial) observer standing far away from a BH horizon will measure an outgoing flux of particles with a nearly thermal spectrum, if the associated quantum field is in the vacuum state as defined by an observer freely falling into the BH, or equivalently the vacuum state as defined in the asymptotic past well before the collapsing matter formed the BH.This Hawking radiation is therefore a purely gravitational effect, such that a BH essentially emits all particle species with masses below its Hawking temperature, T H ≃ M 2 P /8πM for a slowly rotating BH.These remove the mass and angular momentum of the BH (as measured by an asymptotic observer), thus leading to its evaporation.
In this section we will compute in detail how light PBHs evaporate by emitting not only the known Standard Model particles but also a large number of scalar axions, as mentioned above, and we are particularly interested in describing the evolution of its mass and spin.This necessarily involves a precise numerical calculation of the spectrum of emitted particles with different spin, namely the associated "gray-body" factors that quantify the deviations from a purely Bose-Einstein or Fermi-Dirac distribution.These are a consequence of the effective potential probed by the field modes propagating in the BH space-time, which we will assume to be described by the Kerr solution (since any primordial electric charge is radiated away well before it loses a significant amount of mass or angular momentum [53]).These gray-body factors are, in fact, related to the transmission coefficients Γ s l,m for the associated wave scattering problem in the same BH effective potential, and which can be obtained by solving the radial Teukolsky equation [54][55][56][57] describing massless waves of arbitrary spin s.We will use a shooting method to numerically solve this equation (see e.g.[58]).We will then use these results to study the dynamical evolution of a rotating black-hole using the formalism described in [49,50,[59][60][61][62].
For simplicity, in this section we consider geometrized units such that ℏ = c = G = 1 (M P = 1).

Quantum fields in the Kerr space-time
We consider a rotating black hole described by the Kerr solution, which in Boyer-Lindquist coordinates (t, r, θ, φ) reads where M is the hole mass, a = J/M is the BH angular momentum, ∆ = r 2 + a 2 − 2M r, and Σ = r 2 a 2 cos 2 θ.This solution has an inner Cauchy horizon and an outer event horizon at r = r ± (∆(r ± ) = 0).Despite the different equations governing the dynamics of massless (test) fields in curved space-time, in the case of the Kerr metric a unified description can be obtained with the aid of the Newman-Penrose (NP) formalism [63][64][65], and one can condense all these equations into the so-called Teukolsky master equation, which in Boyer-Lindquist coordinates is given by: where the functions Υ s encode the NP scalars, obtained by contraction of the original tensor fields with the Kinnersley tetrad null vectors [64].This equation thus encodes the dynamics of Klein-Gordon scalar fields (s = 0), Weyl fermions (s = ±1/2), vector fields (s = ±1), Rarita-Schwinger fields (s = 3/2) and gravitational perturbations (s = ±2) in the massless limit.Moreover, the Teukolsky equation can be solved by separation of variables, with the NP scalars admitting a mode decomposition of the form: where ω is the perturbation frequency, m is the azimuthal angular momentum quantum number.The angular functions S s (θ) are the so-called spin-weighted spheroidal harmonics, which satisfy the equation [66-69] These functions reduce to scalar spherical harmonics for s = 0 and a = 0, and generalize conventional spin-weighted spherical harmonics to the more general case of axial-symmetry.
s A m l = s A m l (aω) are the eigenvalues of (2.4) and cannot be expressed analytically in terms of the spherical angular momentum quantum numbers l, m.Nevertheless, for aω ≪ 1 they can be computed using a perturbative expansion, yielding: The functions R s are the radial part of the NP scalars and satisfy the radial equation where s Q m l = s A m l + a 2 ω 2 − 2aωm and K = (r 2 + a 2 )ω − ma.These functions take the following form far away and near the BH horizon: where r * is the tortoise coordinate, defined via dr * /dr = (r 2 + a 2 )/∆, and we have imposed ingoing boundary conditions at the horizon.We note that the solutions R s and R −s for the same spin are, in general, distinct, but are nevertheless related through the Teukolsky-Starobinsky identities [70][71][72][73][74] that can be derived from the original field equations: where D = ∂ r − i K ∆ .The Starobinsky constants C s (omitting the l, m quantum numbers for simplicity) for the fields of interest are given by: (2.12) Note that C 0 , C 1/2 , C 1 , and C 3/2 are real, while C 2 is imaginary.
Through the redefinition Y s = ∆ s/2 (r 2 +a a ) 1/2 R s we may also write (2.6) in a Schrödingerlike form: ( where the effective potential V = s V m l vanishes both at the horizon, where r * → −∞, and at infinity, for r * → +∞.The form of Eq. (2.16) guarantees that the Wronskian does not change if calculated at different radial positions.In particular, there is a conserved current: (2.17) Substituting the asymptotic and near-horizon solutions in the Teukolsky-Starobinsky identities yields the relations between R +|s| and R −|s| , while (2.17) yields an energy conservation law, with equal energy flux at the horizon and at infinity.In the wave scattering problem, the transmission coefficient is then given by the ratio between the energy flux into the BH horizon and the incoming energy flux at infinity: .18)Note that this coefficient depends on the frequency, spin and angular momentum quantum numbers of each field mode, as well as on the BH spin parameter, Γ = s Γ l m (a, ω).One can show that the same coefficient yields the corresponding gray-body factor for Hawking emission, since in the latter case it quantities the filtering of field modes by the BH effective potential as they propagate away from the event horizon.

Numerical computation of gray-body factors
An analytical computation of the transmission coefficients is only possible under very stringent approximations [75], so numerical methods are in general required to compute them for different wave modes.Here, we will use a shooting method similar to the one employed in e.g.[58] and first proposed by Starobinsky [75].The first step is to write Eq. (2.6) in terms of the re-scaled radial coordinate x = (r − r + )/r + : where the effective radial potential can be written as: with k = (2 − τ )(ω − mΩ H )r + + x(x + 2)ωr + , and τ = (r + − r − )/r + .Imposing ingoing boundary conditions at the horizon, the near-horizon solutions of Eq. (2.19) can then be expressed in a Taylor expansion [58,76] of the form where ϖ = (2 − τ )(ω − mΩ H )r + and the coefficients a n can be determined by substituting the power series (2.21) in (2.19) and solving iteratively the resulting algebraic equations.The near-horizon solution is then used as a boundary condition for numerically integrating the radial Teukolsky equation up to large distances, where the general form of the solution is known and reads: where ω = ωr + .It is then possible to extract the coefficient s R lm in (ω) in order to evaluate the transmission coefficient.The normalization of the scattering problem is set by setting e.g. a 0 = 1 which is equivalent to This then yields the transmission coefficients for the different spin fields: or, in a more compact way, with . (2.30)

Hawking evaporation in the string axiverse
We determine the evolution of PBHs following the formalism described in [59][60][61] and later in [49,50,62].The PBH mass and spin evolution is determined by the functions F ≡ −M 2 dM/dt and G ≡ −(M/ã)dJ/dt, which remove the dependence on the BH mass.
Here ã = a/M is the BH dimensionless spin parameter.These are given by: where the sum is taken over all particle species i and angular momentum quantum numbers (l, m), x = ωM , k = ω − mΩ H and κ = √ 1 − ã2 /2r + is the surface gravity of the Kerr BH, with Ω H denoting the angular velocity at the event horizon, located at r + .The upper/lower sign corresponds to fermion/boson fields.The function determines whether a black hole spins up or down during its evolution, taking into account the relative magnitude of mass and angular momentum loss rates.If there is a value ã * for which H(ã * ) = 0, the PBH spin parameter will tend to this stable value provided that ∂ ãH| ã=ã * > 0. We note that for ∂ ãH| ã=ã * ≤ 0 the equilibrium point is unstable but that we will not find such cases in our analysis.The differential equations governing the PBH spin and mass evolution can be written in terms of dimensionless variables useful for numerical integration: such that where M i is the initial BH mass, with initial conditions z(t = 0) = 0 and τ (t = 0) = 0.The numerical method described in section 2.2 allow us to compute the gray-body factors for massless fields.In principle one may compute these for massive fields, but given that the emission of particles with masses, µ, above the Hawking temperature T H = κ/2π ≃ 1 GeV(10 10 kg/M ) is exponentially suppressed we work in the approximation where particles are considered massless for T H > µ and are otherwise absent from the emission spectrum.Massless particles as photons, gravitons are emitted since the BH forms alongside all particles with mass µ < T H ∼ few MeV like neutrinos and electrons/positrons, given that this is the natal temperature of PBHs with a lifetime comparable to the age of the Universe, for which M i ∼ 10 12 kg.
As a PBH evaporates its Hawking temperature increases, allowing for the emission of more and more massive degrees of freedom, like muons, tau particles, etc, above the corresponding mass thresholds.The only hadrons with mass below the QCD scale are pions (π 0 and π ± ), these being the only hadronic states included directly in the BH emission spectrum.Temperatures above the QCD scale allow for the direct emission of elementary quarks and gluons that subsequently hadronize.Following [77][78][79][80][81][82][83][84][85][86], we have considered the effective quark and gluon QCD masses given in [87], taking these as threshold values above which each particle is included in the PBH emission spectrum.We note that our results do not change significantly if we consider other values for the effective quark and gluon masses given in the literature, such as in [88].
In order to reproduce the results first obtained in [46], in addition to the Standard Model particles we have considered an arbitrary number of light axions, N a , corresponding to the fraction of the string axiverse with mass below a few MeV.We show, in Fig. 1, our results for the present spin of PBHs, ã0 , as a function of their present mass, M 0 , for different numbers of axions.We consider two limiting cases for the natal PBH spin: ãi = 0.01 (solid curves) and ãi = 0.99 (dashed curves), corresponding to PBH formation in the radiation-dominated era [89][90][91][92] or in an early matter-dominated era [93], respectively.The reader should note that the PBHs of Fig. 1 correspond to the remnants of an initial population of PBHs with nearly the same mass and which are presently at different stages of their evolution, justifying the assumption of a common initial spin.This means that Fig. 1 also depicts the time evolution of the PBH spin, with time flowing from right to left.Note also that the initial mass of PBHs with a lifetime matching the age of the Universe of 13.8 Gyrs depends on the number of emitted species, in particular the number of light axions.In particular, this critical initial mass ranges from 5 × 10 11 kg in the absence of axions to ∼ 2.7 × 10 12 kg for N a = 1000, scaling as N 1/3 a for N a ≳ 10.As one can see in this figure, in the absence of axions (N a = 0, black curves) PBHs lose their spin quite quickly, such that any PBHs with present mass ≲ 10 11 kg should have negligible spin, i.e. spin parameters well below the percent level.A drastic change in this picture occurs in the string axiverse for N a ≳ 100, with PBHs initially spinning up due to the emission of a large number of light scalars (or spin loss being initially halted for initially near-extremal PBHs).When the PBH mass approaches ∼ 10 10 kg and the corresponding Hawking temperature exceeds the QCD scale, the large number of spin-1/2 and spin-1 degrees of freedom emitted starts counteracting the light scalar emission, effectively spinning down the PBH as it evaporates for N a ≲ 400.Above this number of light axions, the PBH spin asymptotes to a non-vanishing value, which tends to the critical value ã ≃ 0.555 originally found in [49] for pure scalar emission as N a → ∞.
Fig. 1 shows that, independently of their natal spin, PBHs with present mass M 0 ≲ 10 11 kg should have a non-negligible spin in the string axiverse scenario (for N a ≳ 100).For the case of initially slowly spinning PBHs, this is particularly relevant, since this spin up due to Hawking emission may render them unstable with respect to superradiant particle creation, namely if the string axiverse includes (as one may expect) heavier axions.As we will analyze in the next sections, this may have a dramatic effect on the PBH spin evolution, making the present PBH mass-spin distribution an even powerful probe of the string axiverse spectrum, with potential directly observable signatures.

Basics of black hole superradiance
Before analyzing how superradiant instabilities may be triggered by PBH evaporation, as suggested by the analysis of [46] reviewed in the previous section, we begin by discussing the basic dynamical features of black hole superradiance neglecting the effects of Hawking emission.Consider then a massive scalar field minimally coupled to gravity, of mass µ, described the action: from which we may derive the corresponding equation of motion in the Kerr metric Eq.(2.1): which reduces to the Teukolsky equation (2.2) for s = 0 in the massless limit.Similarly to the latter, the massive Klein-Gordon equation admits a mode decomposition of the form: where S l,m (θ)e imφ denote scalar spheroidal harmonic functions and now the radial function, R n,l (r) obeys a "massive" Teukolsky equation It is well known that this equation admits quasi-bound state solutions with complex frequencies ω = ω R + iω I , where ω R < µ such that the field is trapped in the gravitational potential well created by the BH [8,76,[94][95][96][97][98][99][100][101][102][103][104].In the non-relativistic limit, where the dimensionless mass coupling is small, the real part of the quasi-bound state spectrum approaches a Hydrogen-like form where α plays the role of the fine-structure constant, a simple consequence of the fact that, when written in the Schrödinger-like form Eq. (2.16) the potential is essentially Coulomblike at large distances from the event horizon (where the scalar field finds support in the α ≲ 1 regime), V (r) ≃ −α/r.As for an electron in a Hydrogen atom, the typical velocity is ∼ α, which justifies denoting α ≪ 1 as the non-relativistic regime.
The imaginary part, ω I , reflects the instability of the bound-states (hence the use of the prefix "quasi-"), with ω I < 0 corresponding to a decay or absorption of the scalar field by the BH, and ω I > 0 to an exponential amplification of the field and of the associated particle number.For α ≪ 1, one finds an approximate analytical expression for the imaginary part of the frequency given by: The superradiant instability thus occurs whenever ω R < mΩ H (α < ã/4 for slowly spinning BHs), leading to an extremely efficient production of particles forming a bound superradiant cloud around the spinning BH.We may regard this as a kind of stimulated emission (even though the process is classical) since all produced particles have the same quantum numbers.In particular, for the fastest growing "2p-state" (n = 2, l = m = 1): Note that the number of particles grows twice as fast, since N ∝ Φ2 .Also, the "2pstate" grows exponentially faster than all others, so we may neglect any other modes in the dynamics of superradiance (neglecting self-interactions, as we discuss below).Thus, each particle produced by superradiance carries one unit of spin from the BH, along side a mass µ, so energy and angular momentum conservation yield: such that the dimensionless spin parameter evolves according to: where in the last step we considered the limit α ≪ 1.The number of particles within the superradiant cloud then follows: where Γ s = 2ω I .It will be useful to note that, for slowly rotating BHs: Note that, strictly speaking, the instability growth rates are computed assuming a fixed BH mass and spin parameter, but since µ ≪ M and M ≫ M P in the regime of interest to our discussion, we may take this a good approximation.The same is true for the semiclassical calculation of the Hawking emission rate, and we may for similar reasons take the two particle production processes as independent, specially since they typically involve different particle species as we discuss below.

A toy model
Given the discussion in the previous subsection, we may now consider the full evolution of a PBH mass and spin taking into account the effects of both superradiance and Hawking evaporation, given by: ) As previously discussed, we are interested in PBHs with a lifetime close to the age of the Universe, i.e. with an initial mass in the range 5×10 11 −10 12 kg, and particularly those born in the radiation era, with initial spins at or below the percent level.Despite their low spins, such PBHs may be superradiantly unstable already at formation, provided there are axions within the string axiverse in the right mass range.In particular, for PBHs with such mass and spin, superradiant instabilities may be triggered for axions with mass µ ≲ 1 MeV, but the axion mass cannot be too low, since the instability growth rate is proportional to µ 9 as given approximately in Eq. (3.12).Note, furthermore, that a significant amount of spin is only extracted from the PBH once the number of particles within the superradiant cloud N ∼ ãM 2 /M 2 P ∼ 10 37 (ã/0.01)(M/10 12kg) 2 , requiring O(100) e-folds of superradiant amplification.This means that superradiance is only efficient for axions roughly in the 0.1-1 MeV mass range for PBHs born in the radiation-era, as illustrated in Fig. 2.
Although there may be string theory compactifications including one or possibly more axions in this mass range, this is certainly not a generic expectation, since axion masses are exponentially sensitive to the magnitude of the non-perturbative effects that generate them.The hundreds or even thousands of light axions expected in realistic string compactifications should have masses distributed throughout a wide range of mass scales.Hence, scenarios with an axion in the mass range shown in Fig. 2 are certainly possibly but not necessarily the most likely, so we will focus our discussion henceforth in scenarios where nearly all axions have masses well below the MeV scale (contributing to the Hawking emission spectrum already at PBH formation), with possibly one extra axion above the MeV scale.The latter will not contribute to the initial Hawking spectrum (although it will once the PBH becomes hot enough), nor will it be produced via the superradiant instability until the PBH spin increases sufficiently as a result of evaporation.
To better understand the dynamical interplay between evaporation and superradiance, we start by considering a toy model where a PBH evaporates through the emission of a single light axion (well below the MeV mass scale), while superradiant instabilities may be triggered for a heavy axion of mass µ ≫ 1 MeV.Although unrealistic, this will help us identifying the main qualitative features of the problem without the intricacies of adding the Standard Model particles across different mass thresholds.
In Fig. 3  The dynamics is thus expected to develop as follows.Initially, while the PBH spin is below the superradiance threshold for a given heavy axion mass, the latter is nonsuperradiant and any quantum fluctuations in the corresponding field are damped by the PBH.However, light scalar emission through the Hawking effect increases the PBH spin until at some point it crosses the threshold for superradiant heavy axion production.Any subsequent quantum fluctuation in the heavy axion field is then expected to be exponentially amplified via the superradiant instability, leading to the growth of a heavy axion cloud around the PBH.We note that the timescales for superradiance and Hawking emission above the threshold differ by several orders of magnitude.For instance, as one can see in Fig. 3, for µ ∼ 100 MeV the superradiance threshold is attained when ã ∼ 0.1 and M ∼ 10 11 kg.Such a PBH will evaporate in ∼ 10 8 years (F(ã = 0.1) ∼ 10 −4 ), while the superradiance e-folding time when e.g. the spin exceeds the critical value by 1% is ∼ 10 −14 s.We illustrate this in Fig. 4, where we plot the Hawking evaporation and superradiance timescales for a heavy axion with µ = 100 MeV and a given PBH spin, as a function of the PBH mass.As one can see in this figure, superradiance is a much faster process for the larger values of the PBH mass, implying that this will be the dominant process determining the PBH mass and spin after the critical spin value yielding ω < Ω H is reached.This difference in the timescales of the two processes may pose a numerical challenge for solving Eqs.(3.13) and (3.14) alongside dN/dt = Γ s N for the number of particles within the superradiant cloud.Nevertheless, we have found that the numerical tools available in e.g.Mathematica are sufficiently accurate for this purpose.An alternative possibility is to artificially reduce the superradiant growth rate via a tunable multiplicative factor and then extrapolate the obtained results to the realistic case.We find that these two methodologies yield results consistent with each other.
A further numerical difficulty is crossing the superradiance threshold, since N decreases exponentially fast in the non-superradiant regime, thus quickly reaching values below numerical precision before the PBH attains the critical spin value through light scalar Hawking emission.This, however, does not correspond to a realistic approach, since it discards the quantum nature of the heavy axion field.Although the development of superradiant instabilities from quantum field fluctuations has not, to our knowledge, been studied in detail so far, it is widely believed that superradiance will amplify any quantum field fluctuations, quickly increasing the corresponding occupation number in the quasi-bound state, so that a classical description is then sufficient to describe the dynamics.
In fact, Kofman showed [105] that Hawking emission populates not only free states, with ω > µ, but also quasi-bound states ω < µ, in a semi-classical calculation similar to the original computation by Hawking.While Kofman's analysis considered only a static BH, so that bound particles produced by Hawking emission are quickly reabsorbed by the BH, in principle it should extend also to the rotating case.The difference for a Kerr BH should reside in the exponential amplification of the bound state occupation number for spin parameters above the superradiance threshold.
In our numerical analysis, we assume this to be the case, and we simulate the effect of bound state quantum emission by first setting N = 0 in the differential equations for the PBH mass and spin evolution, Eqs.(3.13) and (3.14), until just after the superradiance threshold is crossed within our numerical precision.We then take the obtained mass and spin values as initial conditions for the subsequent evolution, where we include the heavy axion cloud starting with N = 1 (changing this initial value somewhat does not significantly affect our results).In the example shown in Fig. 5, we begin with M i = 10 12 kg and ã = 0.01, while the second part of the simulation including a heavy axion with µ = 100 MeV starts with M = 6.30× 10 10 kg and ã = 0.097.
As one can see in this figure, once the superradiant instability is triggered after the critical spin value is attained, the PBH follows closely the superradiance threshold.This is simply due to the fact that the latter occurs on much shorter timescales, quickly depleting the PBH spin until ω = Ω H and superradiant heavy axion production is halted.However, this condition is never fully attained since Hawking emission continuously spins up the PBH due to light axion emission.To better illustrate this, we show in Fig. 6 the time evolution of the PBH spin parameter and of the number of heavy axions in the superradiant cloud for the same example.This shows that the number of heavy axions produced by superradiance grows exponentially fast after the instability is triggered, quickly decreasing the PBH spin back to close to the critical value.As one can observe in Fig. 6, this does not constitute a very significant decrease in the PBH spin, since the superradiant instability is triggered just above the critical value at ã ≃ 4α.In fact, the number of heavy axions increases only until the superradiant term in Eq. (3.14) becomes comparable to the Hawking emission term.At this stage the system reaches a quasi-equilibrium, in which the spin-down effect of superradiance is nearly compensated by the spin-up due to Hawking evaporation.Setting dã/dt ≃ 0 and ã ≃ 4α in Eq. (3.14) then yields the quasi-equilibrium condition: This is analogous to the condition found in [41], although in the latter case the opposite effect was observed since, in the absence of scalar emission, Hawking evaporation tends to spin down the PBH, leading to a reabsorption of the (initially superradiant) cloud in the Γ s < 0 regime.In the present case the cloud remains in the superradiant regime, i.e. with Γ s > 0, so that superradiance produces more and more heavy axions within the cloud as evaporation continues to spin up the PBH.Since the product Γ s N is approximately constant, the number of particles grows linearly in this phase at a rate 4µ(2F −G)) ≃ 4×10 19  axions per second in this example.This quasi-equilibrium configuration is maintained only while the number of heavy axions within the superradiant cloud does change significantly, in this example up until ∼ 10 13 s.After this, superradiance efficiently spins down the PBH, keeping the spin parameter very close to the critical value.
Despite the large number of heavy axions produced until this stage, superradiance has little effect on the PBH mass, which only begins to decrease after ∼ 10 15 s (∼ 30 Myrs), corresponding to the remaining lifetime of the PBH when the superradiant cloud forms.
The subsequent decrease in the PBH mass has two important effects, since it decreases the dimensionless mass coupling α = µM/M 2 P .First, it lowers the critical spin value for which ω = Ω H ; second, it damps the superradiance growth rate Γ s ∝ α 8 .The first effect makes the PBH follow a trajectory in the Regge plane corresponding to the superradiance threshold, as observed in Fig. 5.This holds while superradiance remains faster than evaporation despite the decreasing PBH mass, i.e. down to masses ∼ 10 7 kg.This means that in its final hour (literally in this example) the PBH spins up once more as light scalar emission takes over in the last stages of evaporation.Asymptotically the PBH reaches the stable value ã * = 0.555 yielding H(ã * ) = 0 for pure scalar Hawking emission, as discussed in Section 2. This is not visible in Fig. 5, since it is only attained in the very last stages of the PBH evaporation, beyond the reach of the numerical precision of our simulation.
To summarize our findings in this toy model, a PBH formed with a mass ∼ 10 12 kg evaporates through light axion emission and consequently spins up.After nearly ∼ 14 billion years, its spin surpasses the critical value for triggering a superradiant instability, producing a cloud of heavy axions around it.For most of its remaining lifetime, the PBH is in a quasi-equilibrium configuration with the heavy axion cloud, with evaporation spinning up the PBH nearly at the same rate superradiance spins it down.In our working example the PBH remains in this stage for about 30 million years.At the end of its life, its mass starts decreasing and the PBH follows a Regge trajectory along the superradiance threshold up until its very last stages where evaporation once more increases its spin.
We note that once superradiance becomes inefficient the number of heavy axions within the superradiant cloud stabilizes near the maximum possible value: where the subscript 'c' indicates the PBH mass and spin parameter when superradiance is triggered.This corresponds to converting most of the PBH's angular momentum into heavy axions via the superradiant instability (but fueled by the spin up produced by light scalar Hawking emission).In our example this yields nearly 10 36 axions.
Although the critical PBH-mass spin values for superradiance change for different values of the heavy axion mass, we observe the same qualitative behaviour for all µ > few MeV (recalling that in the 0.1 − 1 MeV mass range superradiance is triggered at PBH formation for ã = 0.01 as discussed earlier).
Our toy model should be an accurate description when the PBH can emit N a ≫ 1 light axions, up to an overall rescaling of the PBH lifetime by a factor ∼ N 1/3 a .

Realistic string axiverse scenarios
With the basic understanding of the main dynamical features of the interplay between superradiance and evaporation in the simplified toy model, we now perform more realistic simulations, with a finite number N a of light axions in the Hawking emission spectrum alongside all the Standard Model degrees of freedom.As described in Section 2, each particle species is included in the emission spectrum once the Hawking temperature exceeds its mass (or effective mass as in the case of quarks and gluons above the QCD scale).
Given our understanding of the evaporation dynamics in the absence of superradiance, the main difference expected between the toy model and more realistic scenarios is the fact that most Standard Model particles have a non-zero spin, therefore carrying away part of the angular momentum of the BH.This means that Hawking emission is overall less efficient in spinning up the BH, and unless the number of light axions is sufficiently large the BH may actually spin down, as discussed in Section 2. The expectation is therefore that superradiant instabilities can only be triggered above a minimum number of light axion species N a .This is illustrated in Fig. 7, where we show the results of our numerical simulations for different numbers of light axions and a heavy axion with µ = 100 MeV.
As one can see in this figure, for different numbers of light axions the superradiant instability is triggered for different values of the PBH mass and spin, although converging to those found in the toy model in the limit N a → ∞.For N a ≲ 100 the superradiant threshold is not crossed for heavy axions with µ ≳ 10 MeV, but since the string axiverse generically predicts hundreds or even thousands of light axions we typically expect instabilities to occur during the PBH evolution if axions in this mass range exist.
As for the toy model, superradiance is initially much faster than Hawking emission in changing the PBH spin, so that after the instability is triggered the PBH follows a trajectory in the Regge plane corresponding to the superradiance threshold ω = Ω H (ã ≃ 4α = 4µM/M 2 P for slowly rotating PBHs).The main difference in realistic scenarios is the fact that we do not observe a spin up of the PBH for low masses, i.e. at the end of its lifetime, as also clear in the time evolution plots shown in Fig. 8, given that Hawking emission is in this case much less efficient in increasing ã than for single scalar emission.Although this may occur when the PBH reaches masses below those that our numerical precision can probe, we may safely conclude that for present PBH masses ≳ 10 6 kg (lifetime exceeding ∼1 s), the PBH distribution in the mass-spin Regge plane should exhibit a single peak at the values (M c , ãc ) at which the instability is triggered and which depend on the string axiverse parameters µ and N a .In particular, the mass of the heavy axion can be inferred from the superradiance threshold condition: (3.17) The dependence on the number of light axions, N a , emitted through the Hawking process is less trivial since it depends on the PBH evaporation dynamics, which has to be computed numerically.In Fig. 9 we show the critical spin contours in the (µ, N a ) plane, from which one can determine N a upon computing µ from Eq. (3.17).We thus find a very unique signature of the string axiverse with hundreds of light (≲0.1 MeV) axions and a single heavy axion (≳ few MeV), corresponding to a sharply peaked spin distribution as a function of mass, with a nearly linear relation between PBH mass and spin for masses below the peak.Moreover, as shown above, the number of light axions and the mass of the heavy axion can be determined from the position of this peak in the Regge plane, so that the full PBH distribution need not be probed across many orders of magnitude in mass.
This shows that measuring the present mass-spin distribution of PBHs below 10 12 kg may have a very significant impact on finding (or excluding) new physics.Methodologies to determine both the mass and spin of a PBH from its photon Hawking emission spectrum have been developed by two of us in [106].Although these may be challenging from the experimental perspective, since they require measuring the PBH photon spectrum close to the primary emission peak energy (where the photon flux is lower than for the secondary component at lower energies), they may be within the reach of future gamma-ray telescopes, as we discuss in section 5.
We note that, in the presence of multiple heavy axions (> few MeV), the first instability to be triggered during the evolution of a PBH corresponds to the lightest of these.The growth of this first heavy axion superradiant cloud will quickly spin down the black hole close to the corresponding superradiant threshold, as we have observed.This will therefore inhibit superradiant instabilities for heavier axions (except in the last fractions of a second of a PBH's lifetime where evaporation may still spin up the PBH).Hence, the shape of the present PBH mass-spin distribution is determined only by the lightest of the heavy axions, being largely insensitive to the existence of other axions.
We also note that our distinction between light axions and the heavy axion refers to the Hawking temperature of ∼ 10 12 kg PBHs at formation.As they evaporate towards their present day mass, the Hawking temperature of these PBHs increases, such that at some stage the heavy axion can also be efficiently emitted.Since we are considering scenarios with N a ≫ 1, the inclusion of one (or even a few) more axion(s) does not significantly change the dynamics, and for simplicity we have kept N a fixed throughout the numerical evolution of the PBH mass and spin.
In our numerical simulations we have considered only free axions, i.e. we have neglected the effects of axion self-interactions, which have been analyzed in detail in [107,108] and also [109] (see also [110][111][112][113][114][115]).The latter considered superradiant axion production around rotating PBHs, although heavier than the ones considered in the present work so that the effects of Hawking emission could be neglected.Axion self-interactions lead, in particular, to 2-2 scattering processes that populate other superradiant and non-superradiant levels in the "gravitational atom" corresponding to the spectrum of BH-axion quasi-bound states.Some axions are "ionized" in these processes, escaping the BH's gravitational potential, which slows down the growth of the dominant 2p-superradiant cloud and may, in fact, prevent its occupation number from growing beyond a maximum number.
The results obtained in [108,109] cannot be easily extrapolated to the case of PBHs with a lifetime comparable to the age of the Universe, given how significant a role we have found PBH evaporation to play in the development of superradiant clouds.We may, nevertheless, try to estimate the parametric regimes in which it is a good approximation to neglect the effects of axion self-interactions, based on the analyses of [108,109].Since we are mostly interested in the non-relativistic regime, we may consider the effects of the leading non-linear term in the axion potential in the resulting Schrödinger-like equation, which has the Gross-Pitaevskii form: where the axion field Φ = (ψe −iµt + c.c.)/ √ 2µ and λ = µ 2 /f 2 a , with f a denoting the axion decay constant.In the limit λ → 0 this corresponds to a Schrödinger equation for a Coulomb-like potential, yielding a Hydrogen-like spectrum of (quasi-)bound states as pre-viously discussed.Self-interactions may thus play an important role when the non-linear term becomes comparable to the energy eigenvalue of the linear Hamiltonian, i.e. when λ|ψ| 2 ∼ α 2 µ 3 for the 2p-state (Φ ∼ αf a ).Since |ψ| 2 represents the axion number density, and the 2p-cloud has approximately a toroidal shape with volume V cloud = 50π 2 /(αµ) 3  (with (αµ) −1 yielding the gravitational Bohr radius) [43,109], we conclude that selfinteractions can be neglected for N ≲ 50π 2 /(αλ).We may then derive a lower bound on f a by taking the maximum number of heavy axions produced in the 2p-cloud when the superradiant instability is triggered by PBH evaporation, N max ≃ ãc (M c /M P ) 2 : where we recall that the subscript 'c' refers to the PBH parameters when the superradiant instability is triggered by evaporation, with 4α c ≃ ãc in the slowly rotating limit.Since ãc ≲ 0.5, given that Hawking emission cannot spin up a PBH beyond this value (which is only achieved for pure light scalar emission), we conclude that for heavy axions with decay constants above the grand unification scale, f a ≳ 10 16 GeV, we may safely neglect the effects of self-interactions in the development of superradiant instabilities.Such large decay constants are, in fact, generic for string axions (see e.g.[9]), thus justifying the free-axion approximation in this context.We note that our dynamical simulations are applicable to any heavy scalar field (µ ≳ 1 MeV) and not only axion-like fields, but the above arguments show that only for very feeble self-interactions may the dynamical effects of the latter be neglected.For instance, neutral pions are similar to heavy axions but interact quite strongly, with λ ≃ 1, as already analyzed in detail in [113].

Direct detection of superradiant axion clouds
In the previous section we have shown that PBH evaporation in the string axiverse may trigger superradiant instabilities for heavy axions due to the emission of hundreds (or even thousands) of light scalar axions and the consequent spin up of (initially slowly-rotating) PBHs.In addition to the unique imprint this leaves on the present mass-spin distribution of PBHs with masses ≲ 10 12 kg, the formation of superradiant clouds may leave a much more direct observational signature, since the produced axions decay into photon pairs.In particular, as we will now describe in detail, an evaporating PBH surrounded by a heavy axion cloud will emit photons as a result of both Hawking emission and heavy axion decay, yielding a unique spectrum.
Hawking emission leads to two types of photons in a PBH emission spectrum.Primary photons are directly emitted by the PBH with a nearly-thermal spectrum (up to the graybody factors discussed in section 2) given by [59,60]: where ω = E γ is the mode frequency (see section 2).In addition, charged particles produced via the Hawking effect also emit photons as they travel away from the PBH, and additional photons also result from the decay of unstable particles like the neutral pion 1 .Such secondary photons are less energetic than their primary counterparts but may nevertheless dominate the emission spectrum at energies below the primary emission peak.
Although the primary spectrum can be computed using semi-analytical methods (computing the gray-body factors numerically as described in Section 2), determining the secondary spectrum typically requires numerical methods of convoluting the primary emission rate for each particle species (analogous to Eq. (4.1)) with their corresponding photon emission rate.We have used the publicly available BlackHawk code [116][117][118][119] to compute both the primary and secondary emission spectra of PBHs with mass and spin satisfying the superradiance threshold condition ω = Ω H , corresponding to the trajectory followed by a PBH after the formation of a heavy axion superradiant cloud of a given mass µ ≳ 1 MeV.We have nevertheless checked that our semi-analytical calculation of the primary emission spectrum agrees with the results obtained using this code.
The latest version of BlackHawk uses two well-known particle physics codes to compute the number of photons radiated by primary particles, namely Hazma [120,121] for primary particle energies below a few GeV and PYTHIA [122] for energies > 5 GeV.PYTHIA code may operate in an extended range via extrapolation tables, but as reported in [123] this may lead to unreliable spectra due to its failure in describing physical processes as the neutral pion decay, π 0 → γγ which should cause a symmetric emission peak centered at half of the pion's mass.We note that the primary emission peak corresponds to photon energies ∼ 5 times the Hawking temperature.For this reason, and taking into account the limits of validity of Hazma and PYTHIA, we employ PYTHIA for PBH masses M < 2.5 × 10 10 kg, while for M > 2.5 × 10 10 kg we use Hazma.
The heavy axions within the superradiant cloud decay into photon pairs with a rate (see e.g.[124,125]): where α EM ≃ 1/137 is the electromagnetic fine structure constant and C aγγ = O(1 − 10) is a model-dependent numerical factor (possibly reaching larger values in some axion models (see e.g.[126]).We may write this as: Note that for the heavier axions this may exceed the present Hubble rate H 0 ∼ 10 −33 eV, i.e. yield axions with a lifetime shorter than the age of the Universe.However, the heavy axion cloud is only formed after ≃ 14 Gyrs, once evaporation spins up the PBH sufficiently to trigger the superradiant instability.It is easy to check that the axion decay rate is always smaller than the PBH evaporation rate when the cloud forms: since the function characterizing the PBH mass loss rate (see Section 2) F ≳ 10 −2 for N a ≳ 100 and, as discussed in the previous section, f a ≳ 10 16 GeV for string axions, taking also into account that superradiance is triggered for α c < 0.15 in all axiverse scenarios.This means that axion decay does not play a significant role in the formation and evolution of the superradiant clouds.It may, however, yield an observable signal as we now show.The corresponding photon emission spectrum is given by: since each of the two photons has approximately half of the axion rest energy (up to subleading gravitational binding energy corrections) and, in the last step, we have replace the monochromatic spectrum by a Gaussian function of width ∆E in order to take into account the effects of a detector's resolution.We then obtain for the maximum photon emission rate from the superradiant axion cloud (at E γ ≃ µ/2), considering the maximum number of axions produced in the evolution, as computed in the previous section: The energy of the axion line is always smaller than the peak of the primary photon emission spectrum of the PBH, since the latter occurs for E γ ≃ 5T H ≃ 0.2α −1 µ ≳ µ.Hence, whether the axion line is detectable depends on the magnitude of the secondary photon emission spectrum from Hawking evaporation.
In Fig. 10 we give two examples illustrating the effect of a superradiant cloud with heavy axions with µ = 100 MeV and 1 GeV on the emission spectrum of PBHs with three different masses and spin.The heaviest PBHs in each case correspond to a PBH where the heavy axion cloud has just formed (and reached its maximum mass), after evaporating with N a = 400 light axions for nearly the age of the Universe.The other two mass and spin values correspond to subsequent stages of the same PBH as it evaporates further and follows the Regge trajectory given by the superradiance threshold condition ω = Ω H (ã ≃ 4α) as discussed in the previous section.We note that in practice one would aim to observe three distinct PBHs presently at different stages of the evaporation process (already dressed with a heavy axion cloud), and not the same PBH at different times, since the evaporation timescale in this mass range is still very large (∼ 1 million years).
In these examples we have chosen f a = 10 16 GeV and |C aγγ | = 10 (or equivalently any combination with g aγγ ≃ 10 −18 GeV −1 ), which maximizes the intensity of the heavy axion line given the constraints obtained from neglecting axion self-interactions discussed above and the typical values of axion decay constants of string compactifications.We see that with a ≃ 2% peak energy resolution the axion line is clearly visible above the secondary In each case the heaviest PBH (black curve) corresponds to the mass and spin values at which the superradiant cloud forms, while the red and blue curves correspond to subsequent stages in the latter's evolution (along the superradiance threshold curve).The energy resolution for the axion line is taken to be 1% of the axion mass in each case.
photon emission from the PBH evaporation for the PBH mass values considered.Although these examples may be somewhat optimistic, it is quite remarkable that such an axion line is observable for such low values of the axion-photon coupling.Note that in the case of heavier axions, for which the instability is triggered at higher spin values, the axion line is more pronounced as given by Eq. (4.6).
It is worth remarking that detecting a slowly rotating black hole with a mass ≲ 10 11 kg, which must in principle be a remnant of the evaporation of a heavier PBH2 , exhibiting a monochromatic line in its electromagnetic emission spectrum would be evidence for the existence not only of a heavy axion but also of hundreds of light axions, as otherwise it could not have developed a large enough spin to trigger the superradiant instability (recall that for N a = 0 any natal spin is quickly lost, as can be seen in Fig. 1).

Conclusions
In this work we have considered the evaporation of PBHs in the context of the string axiverse, following on the seminal work in [46].The generic prediction of hundreds or even thousands of light scalar axions in realistic string scenarios has a tremendous impact on the dynamics of small PBHs, since light scalar emission tends to spin up a BH, as opposed to the emission of particles with non-zero spin.This is due to spin zero particles being the only particles that can be emitted in the spherically symmetric l = 0 mode, i.e. without carrying away the BH's angular momentum.As shown in [46] and revised in detail in Section 2 of the present work, an initially slowly rotating PBH (ã ≲ 0.01) can spin up up to values ã ∼ 0.1 − 0.5 for N a ≳ 100 light axions.
This increase in a PBH's angular velocity, which for PBHs born with ∼ 10 12 kg occurs on timescales comparable to the age of the Universe, has an important consequence that we have explored in detail in this work -it may trigger superradiant instabilities.The string axiverse typically includes axions with masses spread out over several orders of magnitude [9], most of which are likely below the MeV scale and hence included in the PBH Hawking emission spectrum for the above-mentioned natal mass range.However, there may be one or more axions with a larger mass, and which can be produced via the superradiant instability once a PBH reaches a critical spin value as a result of evaporation.
The dynamical interplay between Hawking evaporation (with light sub-MeV axions) and the superradiant instability (producing heavy super-MeV axions in clouds gravitationally bound to the PBH) is quite interesting, given in particular the very different timescales of the two particle production processes.As we have shown in this work, once evaporation spins up a PBH above a certain critical spin, the superradiant instability quickly amplifies any quantum fluctuation in the heavy axion field, and the expense of reducing the PBH's spin back to the critical value.On a longer timescale, the PBH continues to spin up due to light axion emission, therefore feeding the superradiant instability and the heavy axion cloud.These two opposing effects keep the PBH-axion cloud system in a quasi-equilibrium state with nearly constant spin for a long time, and as the PBH mass decreases it follows a very simple Regge trajectory (mass-spin plane) corresponding to the superradiance threshold for the heavy axion ã ≃ 4µM/M 2 P .Towards the end of the PBH's lifetime superradiance becomes less and less efficient in extracting the PBH spin, as a consequence of the decreasing dimensionless mass coupling α = µM/M 2 P .The number of heavy axions in the cloud stabilizes near the maximum value N max ≃ ãc (M c /M P ) 2 , where the subscript 'c' denotes the PBH parameters when the superradiant instability is triggered, as supported by our numerical simulations.Evaporation then takes over as the main mechanism driving the PBH evolution and therefore increasing its spin for a sufficiently large number of light axions.Numerically, we can only observe this final spin up in the toy model with pure scalar Hawking emission, given that numerical precision limits the considered PBH mass range to M 0 ≳ 10 6 kg3 in a toy model with pure scalar emission.This toy model mimics what happens in the limit N a → ∞.
This thus leads to a striking prediction for the present mass-spin distribution of PBHs in the range 10 6 − 10 12 kg.On the one hand, for the heavier ones that are still spinning up due to light axion emission, the spin parameter should decrease with the mass (the exact function depending on the number of light axion species, N a ).On the other hand, for the lighter PBHs that have already formed a heavy axion cloud, the spin parameter should increase linearly with the PBH mass, along the Regge trajectory corresponding to the superradiance threshold ã ≃ 4µM/M 2 P .This gives a peaked mass-spin distribution (see Fig. 7), the mass and spin of the most rapidly rotating PBH depending on the number of light axions, N a , and the mass of the heavy axion µ (see Fig. 9 and associated discussion).
In addition to this indirect signature of the string axiverse, the presence of a superradiant axion cloud can in principle be directly detected as a single emission line on top of the PBH's Hawking emission spectrum, located at approximately half of the heavy axion's mass (since axions decay into photon pairs).Although we have not performed a detailed analysis of the detectability of this axion line, we have shown that its intensity can be comparable to that of the PBH's (secondary) Hawking emission for axion-photon couplings as low as g aγγ ∼ 10 −18 GeV, corresponding to axion decay constants of the order of the grand unification scale, f a ∼ 10 16 GeV, typical of string axions, up to an O(10) model-dependent coefficient C aγγ .This feature is quite unique, since PBHs with different masses and spins, presently at distinct evaporation stages, should exhibit the same axion line despite their different Hawking emission spectra if they have grown a superradiant cloud around them.
Both indirect and indirect signatures of the string axiverse depend intrinsically on detecting and accurately measuring a PBH's photon emission spectrum, since in addition to the axion line this allows for a determination of both its mass and spin, following e.g. the methodologies devised in [106].In particular, the latter require determining specific features in the spectrum close to the primary emission peak, where, as illustrated in Fig. 10 the emission rate is lower.As discussed in [106], the sensitivity of planned gamma-ray telescopes such as the All-Sky-ASTROGAM [128,129] or AMEGO [130,131] missions may not be sufficient for these purposes unless we can find light PBHs (≲ 10 12 kg) at a distance below 100 AU of the Earth, which although not impossible is unlikely given current bounds on their abundance [132].However, the proposed MAST mission [133], with an unprecedentedly large detector area, may potentially reach enough sensitivity.
An important question also comes out of our analysis in this work − what happens to the heavy axion clouds once the PBHs evaporate away?The analysis of superradiant dark matter production by light PBHs (< 10 6 kg) performed in [41] has suggested (although not rigorously proven), that superradiant clouds may survive black hole evaporation as self-gravitating, microscopic boson stars.The main idea is that, as a PBH evaporates, its gravitational potential (which bounds the scalar cloud) decreases in time, first adiabatically (compared to the timescale of the Hydrogen-like wave function), but speeding up towards the end of the PBH's lifetime so that the PBH suddenly vanishes -much like a quantum quench.Using the results obtained in [41], we find that PBH evaporation should only become non-adiabatic when the PBH reaches a value: where M cloud = µN is the total mass of the axion cloud.This then suggests that the heavy axion field profile should slowly evolve from a superradiant cloud around a PBH to an essentially self-gravitating configuration well before the PBH fully evaporates away.PBH evaporation could thus leave behind microscopic axion stars!Note that the cloud expands from an initial size of a few times the gravitational Bohr radius ∼ M 2 P /M c µ 2 to the much larger size of the self-gravitating configuration, ∼few×M 2 P /M cloud µ 2 , given that the axion cloud only contains in general a small fraction of the PBH mass when it forms, i.e.M cloud ≪ M c .Note also that this should result in a rotating boson star by angular momentum conservation, but that these configurations are unstable and end up decaying into non-rotating spherical stars [134][135][136].
Showing that superradiant clouds may indeed become self-gravitating states requires dedicated numerical simulations, given the intrinsically non-linear nature of the problem, and which are beyond the scope of this work.Nevertheless, it is interesting to speculate about the possibility of directly observing such a transition, since after its final Hawking explosion, a PBH could leave behind a compact object (the "axion star") with a monochromatic gamma-ray spectrum, as computed in the previous section.
Whether or not sufficiently sensitive telescopes will become available within the foreseeable future to detect all the effects proposed in this work, this demonstrates the enormous potential that evaporating PBHs can have as probes of beyond the Standard Model physics, in particular the string axiverse.We can only hope that the Universe has been kind enough to provide us with a sufficiently large number of these fascinating compact objects.

Figure 1 .
Figure1.Present PBH spin, ã0 , as a function of the corresponding present mass, M 0 , for an initial population with spin ãi = 0.01 (solid curves) or ãi = 0.99 (dashed curves).Curves are labelled (and coloured) according to the number of light axions (≲ few MeV) included in the PBH emission spectrum.

Figure 2 .
Figure 2. Value of the axion mass for which superradiant instabilities are triggered at PBH formation, for PBHs with a lifetime comparable to the age of the Universe born with ãi = 0.01.
we show the PBH spin as a function of its mass considering only the effects of single scalar Hawking emission, obtained by solving numerically Eqs.(3.13) and (3.14) for an initial PBH mass M i = 10 12 kg and spin ãi = 0.01.In this figure we also give the curves in the PBH mass-spin plane corresponding to the superradiance threshold ω = Ω H for different heavy axion masses.

Figure 3 .
Figure 3. PBH Regge trajectory for single-scalar Hawking emission (HE) for M i = 10 12 kg and ãi = 0.01, and superradiance threshold curves for different heavy axion masses, as labelled.

Figure 5 .
Figure 5. Black hole trajectory in the mass-spin "Regge" plane through Hawking emission and superradiance (solid orange curve) for M i = 10 12 kg, ãi = 0.01 and a heavy axion with mass µ = 100 MeV.Also shown are the trajectory in the absence of superradiance (dotted blue curve) and the superradiance threshold (dashed blue curve).

10 -9 10 -Figure 6 .
Figure 6.Numerical evolution of the PBH spin (left) and number of heavy axions in the superradiant cloud (right) for the same parameters of the example shown in Fig.5.In these plots time is measured from the onset of the superradiant instability.

Figure 7 .
Figure 7. Black hole trajectory in the mass-spin "Regge" plane through Hawking emission and superradiance (solid curves) for M i = 10 12 kg, ãi = 0.01 and a heavy axion with mass µ = 100 MeV, for different values of the number of light axions N a in the Hawking emission (HE) spectrum, as labelled.Also shown are the trajectory in the absence of superradiance (dotted curves) and the superradiance threshold (dashed curve).

Figure 8 .
Figure 8. Numerical evolution of the PBH spin (left) and number of heavy axions in the superradiant cloud (right) for the same parameters of the example shown in Fig.7.In these plots time is measured from the onset of the superradiant instability.

Figure 9 .
Figure 9. Contours of the critical spin parameter ãc at which the superradiant instability is triggered as a function of the number of light axions N a and the mass of the heavy axion µ.This corresponds to the present maximum spin of PBHs in the range 10 6 − 10 12 kg (slowly rotating at birth).In the black region no superradiant instabilities are triggered.

= 7 ×= 1 . 4 ×ΔE = 10 MeVFigure 10 .
Figure 10.Photon emission spectrum of PBHs in the presence of a heavy axion with 100 MeV (left) or 1 GeV (right) mass and f a /C aγγ = 10 15 GeV (g aγγ ≃ 10 −18 GeV).In each case the heaviest PBH (black curve) corresponds to the mass and spin values at which the superradiant cloud forms, while the red and blue curves correspond to subsequent stages in the latter's evolution (along the superradiance threshold curve).The energy resolution for the axion line is taken to be 1% of the axion mass in each case.