Astrophysical Searches and Constraints on Ultralight Bosonic Dark Matter

Starting from the evidence that dark matter indeed exists and permeates the entire cosmos, various bounds on its properties can be estimated. Beginning with the cosmic microwave background and large scale structure, we summarize bounds on the ultralight bosonic dark matter (UBDM) mass and cosmic density. These bounds are extended to larger masses by considering galaxy formation and evolution, and the phenomenon of black hole superradiance. We then discuss the formation of different classes of UBDM compact objects including solitons/axion stars and miniclusters. Next, we consider astrophysical constraints on the couplings of UBDM to Standard Model particles, from stellar cooling (production of UBDM) and indirect searches (decays or conversion of UBDM). Throughout, there are short discussions of"hints and opportunities"in searching for UBDM in each area.

Without going into the specifics (see Ref. [1]), it suffices to say that only channel 1 produces UBDM with the required properties as outlined in the chapter "Introduction to Dark Matter." Production channels 2 and 3 produce hot DM, or indeed dark radiation, each of which are strongly constrained by the CMB anisotropies [2,3].
In channel 1a (vacuum realignment) the UBDM relic density is a function of two parameters, ( , ), where is the initial field displacement, i.e. the location of the field in its potential relative to the minimum at "the initial time" (in practice, at the end of inflation). In this scenario, the initial field displacement is taken to be completely uniform throughout space, this state of affairs having been arranged by the same mechanism that causes the large scale observed homogeneity of the cosmic microwave background (CMB), inflation, or otherwise. The correct relic abundance can be achieved across many orders of magnitude, covering all the masses of interest (10 −33 eV, 10 −1 eV) for ≤ . For an axion-like particle (ALP) the fundamental parameter from theory is : the scale of spontaneous symmetry breaking, also called the axion decay constant. The parameter , defined via ≡ , is the initial angle that the axion field takes (recall that the axion is the phase of a complex field). At early times the axion possesses a shift symmetry, → + constant, and thus has no preferred value and can be considered a free random variable (although very small values, or values very close to are considered fine tuned). Because is undetermined there is a wide range of allowed values for the fundamental parameters ( , ) consistent with the required relic density. In particular, in this channel large values of the decay constant at the grand unified scale (∼ 10 16 GeV) or the reduced Planck scale (∼ 10 18 GeV) are allowed.
Production via channel 1b (topological defect decay) is possible only for UBDM that is a Goldstone boson of a spontaneously broken global symmetry (the "Kibble-Zurek mechanism" [4,5] described in Section 3.2). In particular it applies to the QCD axion and other ALPs, where topological strings and domain walls are formed when the global U(1) symmetry is spontaneously broken. If symmetry breaking occurs after inflation, then the defects cannot be smoothed out and inflated away, and the axion field takes on a very inhomogenous distribution (in contrast to the case of vacuum realignment). The defects later decay when non-perturbative effects give the ALP a mass. This process must be simulated using classical lattice field theory, and has only been studied in detail for the QCD axion [6][7][8]. Large numerical uncertainties related to extrapolation to physical couplings prevent an agreed estimation of the relic density. The correct relic abundance can be achieved within numerical and model uncertainty (extrapolation, domain wall number, explicit symmetry breaking) for all values of 10 12 GeV [9]. The production mechanism channel 1b works for < max , where max is the maximum thermalization temperature of the Universe, and the bound arises since defects only form if symmetry breaking occurs during the ordinary thermal history of the Universe. max is bounded from above due to observational constraints on the theory of inflation. In particular , the inflationary Hubble rate, is bounded from above by the fact that tensortype CMB anisotropies have relative amplitude 0.1 compared to scalar-type perturbations leading to the constraint 10 14 GeV. sets the temperature of the Universe during inflation to be the Gibbons-Hawking temperature, GH = /2 . The maximum thermalization temperature could actually be larger than this, which can easily be seen from the Friedmann equation during radiation domination, 3 2 2 = 2 ★ 4 /30, where the quantity ★ counts the effective number of relativistic degrees of freedom [10]: where is the degrees of freedom of species (e.g. two polarizations for the photon) and is the temperature of species , and is the photon bath temperature. The value of ★ at very high temperatures is bounded from below by the Standard Model contribution, ★,SM = 106.75. monotonically decreases, and so max = . If reheating after inflation is instantaneous and 100% efficient, we find an upper bound for max 8 × 10 15 GeV. ALPs with values of larger than this upper bound on max cannot be produced by mechanism 1b, and must be produced by mechanism 1a. The observational lower bound on max arises from demanding successful Big Bang nucleosythesis, max 1 MeV. For values of in this very large range of allowed max values, it is not = 1/ √ 8 is the reduced Planck mass, related to the mass scale given in the "Units and Conversions" section by the factor of √ 8 coming from Einstein's equation in general relativity.

The CMB and linear structure formation
Considering how the gravitational effects of DM dominate the formation of structure in the Universe, one can derive bounds on the UBDM properties from the theory of cosmological structure formation in general relativity [11]. Consider a flat, homogeneous, and isotropic spacetime described by the Friedmann-Robertson-Walker metric: The scale factor is ( ), which obeys Friedmann's equation for the Hubble rate ( ) = / : where¯is the total, spatially averaged, energy density. is composed of photons, "baryons" (by convention in cosmology we do not separately consider the small mass density of electrons), neutrinos, DM, and the cosmological constant or dark energy. Objects "on the Hubble flow", i.e., feeling negligible local gravitational potentials, appear to recede from an observer at the origin with a velocity =ˆ, where andˆare the  Fig. 1 The evolution of cosmic quantities as a function of scale factor or temperature. We show the evolution of the Hubble parameter (red line, left axis) and the comoving Hubble radius (blue line, right axis) together with various relevant cosmological events. The blue shaded area approximately encompasses the large scale structure (LSS) of the Universe, while grey shaded areas indicate where QCD axion (with ∈ [10 6 GeV, 10 18 GeV]) and fuzzy dark matter (FDM) start to become dynamical. Note that the temperature scale on the top is not exactly regular due to the scaling with the number of relativistic degrees of freedom for entropy, ★, . The quantity ★, gives the number of effective relativistic degrees of freedom contributing to the entropy density; ★, takes the same form as Eq. (1) with the fourth powers replaced by cubes (see, e.g., Ref. [10], Chapter 3).
distance and direction from the observer to the object, respectively. We begin with a Newtonian approximation to cosmology (see e.g. Ref. [12]). Consider an observer at the origin, and a single particle of UBDM on the Hubble flow. The UBDM de Broglie wavelength is = 1/( ) = 1/( ), which gives the radial position uncertainty, Δ . A net gravitational force in the positive direction along the line of centres between the observer and the UBDM requires Δ ⇒ ( ) −1/2 , which defines a critical separation crit = ( ) −1/2 . On average, UBDM separations larger than crit undergo gravitational clustering, and those smaller than it do not. The cosmological horizon size is approximated by the Hubble length scale = −1 . In order for UBDM to have any inhomogeneous gravitational effect within this radius requires crit < . We show the cosmological evolution of = −1 , and the related comoving Hubble radius ( ) −1 , as functions of temperature and scale factor in Fig. 1. The bounds and other cosmological events mentioned in what follows can often be read off directly from that figure, and we will occasionally highlight this fact going forward.
UBDM violating this bound does not cluster within our cosmological horizon, is thus indistinguishable from the cosmological constant, and will not concern us in this book.
Very light scalar fields which are homogeneous on the scale of the cosmological horizon provide models for dark energy. The simplest such models are described by a canonical kinetic term in the Lagrangian, and a scalar potential ( ), and are known as "quintessence" [13]. An ultralight bosonic field with a mass less than the bound Eq. (4) is one such very simple model, with ( ) = 2 2 /2. More complex models invoke different potentials, or more fields, or even generalize the kinetic term, at which point they cross over into theories of "modified gravity" and "beyond Horndeski" scalar-tensor theory [14].
Assuming that UBDM constitutes the entirety of the DM, we can extend the bound to any redshift of interest where we know that DM exerted a discernible gravitational effect by simply substituting the Hubble parameter at that redshift. For temperatures below about 1 MeV we can use the expression for the Hubble parameter [11]: where the second equality defines the energy function ( ). The quantities Ω and Ω Λ are the density parameters of matter and the cosmological constant, defined as the density divided by the critical density, i.e. Ω =¯/ crit and crit = 3 2 2 0 . The last term in the brackets arises from the radiation energy density, which is defined relative to the matter density via the redshift of matter-radiation equality, eq . The epoch of matter-radiation equality can be found via the relative redshifting of matter and radiation components: (1 + eq ) 3 = (1 + eq ) 4 , with the density parameters defined today. CMB observations fix eq ≈ 3390, and it is thus slightly earlier in cosmic history than decoupling, dec ≈ 1100.
The baryon acoustic oscillations (BAO, see Sec. 1.1) observed in the CMB and galaxy surveys like the Sloan Digital Sky Survey [15] require that DM was gravitationally relevant at and before matter-radiation equality: if it was not, because baryons are coupled to the photons at early times and perturbations in them cannot grow in the radiation era, the amplitude of galactic fluctuations on scales of order 1 Mpc would not not be consistent with the amplitude and scale dependence of the CMB anisotropies. Again assuming that UBDM is all the DM and substituting ( eq ) we arrive at the tighter bound (cf. Fig. 1): where we have neglected the small contribution of Ω Λ at equality, and taken reference parameters from the CMB+BAO combination in Ref. [3]. The matter-radiation equality bound, Eq. (6), is the UBDM equivalent of saying that DM is not "hot" [16]: gravitational clustering is required before matter-radiation equality in order for bottom-up hierarchical structure formation (rather than top-down fragmentation) of galaxies, consistent with observations of extremely high redshift galaxies. We could progress further with such estimates (and we will in due course), but now we must make our model more precise.

Tutorial: The Growth of Cosmic Structure
The challenge in cosmological perturbation theory [17] is to compute the transfer function, ( , ) for the mode evolution of each cosmological species (baryons, photons, neutrinos, dark matter) with Fourier wavenumber , which fully specifies linear evolution of cosmological fields from Gaussian initial conditions. That is: where , ( ) is the initial condition of the field and is a Gaussian random field defining the initial correlation functions of the field .
The codes [18] and [19] are the standards for numerical computation for CDM (and many other things), while CAMB [20] can be used for UBDM that is a real scalar field with the self-interaction potential approximated by ( ) = 2 2 /2. This tutorial gives a brief overview of the most relevant aspects of cosmological perturbation theory for UBDM constraints.
Using these reference parameters further assumes that UBDM is sufficiently CDM-like that we can use the standard CMB parameters (which are derived under the assumption of ΛCDM).
Cosmological perturbation theory deals with the evolution of fluctuations relative to a homogeneous and isotropic background. Background quantities are labeled with an overbar, since they represent the spatial average, and thus depend only on cosmic time . The perturbation modes have spatial dependence captured by their wavenumber, and perturbations at the initial time all have relative amplitude much less than one with the respect to the background quantities. The fields of interest are the components of the energy momentum tensor, written as 0 0 = −(¯+ ), = (¯+ ) + Σ , 0 = (¯+¯) , which defines the energy density, , pressure, , and heat flux, = ∇ · , and we assume anisotropic stresses Σ vanish. This gives the fields = /¯, and , while pressure is typically described it terms of a sound speed, 2 = / .
Next, perturb the metric Eq. (2), and switch to conformal time, , via d = d . The Newtonian gauge considers only scalar metric perturbations: The potential Φ is the usual Newtonian potential, and Ψ is the curvature perturbation: they are equal in the non-relativistic limit. The energy momentum tensor is coupled to the metric degrees of freedom by the Einstein equation: where is the Einstein tensor, and depends on the metric potentials and their derivatives. This is the dynamical equation determining the evolution of the metric.
The equation of motion for the UBDM field with self-interaction potential ( ) is: where the d'Alembertian ( ) is: and where and are the metric determinant and the inverse of the metric, respectively. Setting = 1   2   2 2 for simplicity, this leads to the equations of motion for the UBDM background field,¯and fluctuation mode :¯ where primes denote derivatives with respect to conformal time, and H = / = . For the UBDM field, we find = /( ) by variation of the action with respect to the metric tensor, giving: Working to first order in the metric perturbations and , and with potential = 2 2 /2 the components are: CDM is defined as a collisionless and uncoupled fluid, = 2 = 0. Baryons have a sound speed, 2 ≠ 0 (computed from the evolution of the baryon temperature), equation of state = 0 (on average the baryons have negligible pressure), and are coupled to photons via Thomson scattering. The photon equation of motion is derived from the Boltzmann equation, which is expanded in Legendre polynomials to capture the dependence on the angle between the momentum coordinate on phase space and the wavevector. The hierarchy of moment equations are labeled by the order ( ) of the Legendre polynomial: the zeroth moment gives the equation of motion for the density, the first, for the velocity, the second, for the anisotropic stress, and so on (a recursion relation can be used to approximately close the hierarchy above some max ). Truncating this Boltzmann hierarchy at the velocity moment, the photons resemble a fluid with = 2 = 1/3, collisionally coupled to the baryons. We consider perturbations to the energy density = /¯and heat flux , defined via¯(1 + ) = ( 0 ) , where ( ) is the energy momentum tensor.
Let us now consider a number of limits of the full equations of motion, which can be found in Ref. [17]. At early times, photons have enough energy to keep hydrogen and other atoms ionized, giving rise to a large free electron density. Thus, the photons and baryons are tightly coupled by Thomson scattering and can be treated as a single fluid with = . Considering only sub-horizon modes ( ), and using the Poisson equation and the pressure component of the Einstein equation Eq. (9), the photon fluid at early times obeys the equation of motion: where the photon sound speed is , = 1/ √ 3 (speed of pressure perturbations in a gas of photons in thermodynamic equilibrium). At very early times all in the driving term on the right hand side can be neglected. Then this equation has sound wave solutions for > (16 2 ) 1/2 = √ 6 . This defines the Jeans scale of the photon-baryon fluid, which is of order the comoving horizon size. Perturbations with wavelength shorter than the Jeans scale undergo coherent, pressure supported oscillations. Perturbations with wavelength longer than the Jeans scale grow due to gravitational instability. The sound waves prevent the formation of gravitationally bound structures in the photon-baryon fluid and lead to BAO. At recombination temperatures of around 0.2 eV (redshift ≈ 1100) [10], the energy of the ambient photon fluid is no longer sufficient to keep neutral hydrogen from forming. At this time, the free electron density drops to zero, the photon-baryon fluid decouples, and the sound wave stalls. This sound horizon for the BAO is given by: where , is the baryon sound speed in the plasma, 0 is the time today, and rec is the time at recombination when , drops rapidly from , to zero. The BAO scale leads to oscillations in the CMB angular power spectrum, which we have seen already in Chapter 1. The gauge invariant temperature anisotropy of the CMB is given by: where is the Thomson scattering opacity, is the baryon velocity,ˆis a unit vector giving the sky position, and the integral is along the line of sight. The four terms in Eq. (22) correspond, respectively, to: the gravitational redshift, the photon anisotropy, the Doppler effect, and the final term gives rise to the integrated Sachs-Wolfe effect, which is an additional form of gravitational redshift. Decoupling occurs at a redshift dec ≈ 1100, which gives the angular scale of the first CMB acoustic peak. The driving term on the right hand side of Eq. (20) is dominant for < eq ≈ 3400, corresponding to angular scales slightly smaller than the first peak, including the second and third peak. Thus, the relative heights of these peaks can be used to measure the matter content and its behaviour near matter-radiation equality.
How do UBDM perturbations evolve? The first transition in behaviour is in the equation of state, which becomes zero (i.e. pressureless) shortly after ( osc ) = (this defines the value of the scale factor osc when the background UBDM field,¯, begins to undergo coherent oscillations, see Problem 1). Prior to this time the UBDM is relativistic and perturbations cannot grow. For the UBDM perturbations Eq. (13) can be approximated as a fluid with sound speed [24] The non-relativistic limit of this expression is derived later on in this Chapter from the Schrödinger-Poisson equation, see Section 2.2 and Problem 2. Now compare the behaviour of UBDM and CDM+baryons for subhorizon modes in the matter-dominated era. The baryon sound speed can be neglected after decoupling, so CDM and baryons can be combined into a single pressureless fluid. In the sub-horizon , super-Compton limit, the CDM+baryon and UBDM fluids obey the coupled equations of motion: This equation ignores the effect of gravitational lensing along the line of sight. This second-order effect is important at high multipoles and is sensitive to the UBDM sound speed and structure growth. See Refs. [21,22].
For an axion-like potential, the equation of state is = −1 prior to osc . For a scalar field with potential = 2 2 /2 + 4 , the equation of state is = 1/3 at early times for large initial conditions. For a complex scalar, the early time equation of state is = 1 due to the conserved charge and Goldstone mode [23]. In each case, perturbations are suppressed relative to pressureless CDM.
This expression is exact in the UBDM comoving gauge. Additional terms due to the gauge transformation to a standard gauge, e.g., Newtonian or synchronous, decay on sub-horizon scales as all gauge artifacts do in cosmological perturbation theory [20]. Setting¯U BDM = 0 and substituting the Poisson equation Eq. (26) into Eq. (24) gives the solution + = + ( ) + − ( ) −3/2 . The growing mode initial conditions sets − ( ) = 0, and the inflationary initial conditions and matter transfer function fix + ( ). Due to the zero pressure and sound speed of CDM, all the -dependence in the solution is fixed by the initial conditions, and the dynamics are scale invariant. Now consider a UBDM-dominated universe by taking¯+ =¯ ¯U BDM (i.e., no CDM and treating the baryons as sub-dominant) in Eq. (25), and again substituting the Poisson equation. The substitution of the Poisson equation gives rise to a negative contribution on the left hand side proportional to UBDM , which drives growth of UBDM , while the positive contribution from the sound speed term leads to acoustic oscillations. The sign of the term proportional to UBDM depends on and as such different modes evolve differently. That is, we find Eq. (25) exhibits a Jeans scale, , separating growing/decaying and oscillating modes. The exact solution for pure UBDM is UBDM = + ( ) + ( , ) + − ( ) − ( , ), where the growth functions are: Consider the evolution of three wavenumbers in the pure UBDM case: the horizon size, ★ = ; the Jeans scale, = √ ; and the Compton scale, = . The Compton scale defines relativistic modes where 2 UBDM = 1; increases with time, and more modes become non-relativistic. If ★ < , then a mode is non-relativistic when it enters the horizon and behaves as CDM ("long modes"). If a mode is relativistic when it enters the horizon ("short modes") then the sound speed cannot be neglected, and modes will not grow until the later time when the Jeans wavenumber enters the horizon. The evolution of these three modes is illustrated in Fig. 2. All modes intersect at the time osc , which defines the special mode , the horizon size when the UBDM background becomes non-relativistic. All < evolve similarly to CDM. All > have suppressed growth.
The scale that determines suppression of growth compared to CDM is the Jeans scale at matter-radiation equality. Using Eq. (25) in the pure UBDM limit with UBDM ≈ 2 /4 2 2 , substituting the Poisson equation, and solving for where the effective mass term in the oscillator equation for the overdensity vanishes, we find: Recall that by definition CDM has zero sound speed. Thus CDM possesses no Jeans scale (the growing mode solution above is scale invariant), and we see that UBDM is only equivalent to CDM exactly in the limit → ∞. In practice, they are equivalent as long as does not play a role in any observation. An observable related to the matter clustering is the matter power spectrum defined by is the total matter (baryon+CDM+UBDM+neutrino) overdensity, and is the Dirac delta distribution. The presence of the sound speed and consequent Jeans scale for UBDM leads to a suppression of ( ) relative to CDM at large wavenumbers. A fit for the relative suppression in ( ) for UBDM with ( ) = 2 2 /2 versus CDM is [26]: ( ) = 1.61 For the mixed CDM-UBDM system, the behaviour of ( ) can also be derived [27,28]: perturbations with > experience a finite amplitude suppression which increases with the ratio Ω UBDM /Ω . intersect at osc when the field begins to oscillate. This defines the scale of power suppression as the comoving horizon size at this time, = osc osc = ( osc. ) −1 . Due to the slow evolution of with , and the logarithmic growth of density perturbations during the radiation epoch, the suppression scale is also approximated by the Jeans scale at matter-radiation equality. Adapted from Ref. [25].

End of Tutorial
As we have just seen in the above tutorial, two effects distinguish UBDM from other ingredients in the ΛCDM model: (1) the background expansion rate, ( ), driven by the transition in the equation of state UBDM at the epoch osc , and (2) the growth of perturbations, driven by the gradient energy in the Klein-Gordon equation and manifested as an effective sound speed, 2 UBDM . Depending on the value of osc , the change in ( ) affects different CMB multipoles. This can be understood by considering Eqs. (20) and (21) in the tutorial. First consider UBDM violating the bound Eq. (6). We know such UBDM must be a sub-dominant component of the DM. How does the CMB tell us this? Such UBDM changes the expansion rate after matter-radiation equality. This changes the distance to the surface of last scattering, and the angular size of the BAO in the CMB: it moves the first acoustic peak from its observed position ℓ ≈ 200. This can be compensated by a change in the value of the Hubble constant, 0 . After such a compensation there is a residual integrated Sachs-Wolfe effect which differs from ΛCDM. If ≠ 0 in the post-recombination Universe, then the gravitational potential Φ ≠ 0 into Eq. (22). Due to the fact that the equation of state UBDM ≠ 0, −1 (the two available equations of state in ΛCDM), the evolution of Φ is different in the presence of a small contribution of UBDM, and the shape of the ℓ < 200 CMB multipoles is very sensitive to the value of Ω UBDM . Now consider UBDM satisfying the bound Eq. (6). The change in the expansion rate compared to ΛCDM now occurs during the radiation dominated epoch. The horizon size at the time osc was smaller than one degree on the sky, corresponding to multipoles ℓ > 200, i.e. the higher acoustic peaks. UBDM changes the distance This is one of the ways the CMB is used to constrain the equation of state of dark energy. UBDM effects on the CMB temperature power spectrum. UBDM changes the expansion rate compared to CDM in the early radiation dominated epoch, 3000, which affects the damping of the BAO, visible through the heights of the power spectrum peaks at large multipoles. By eye, it is clear that the Planck data strongly exclude UBDM with ≤ 10 −26 eV. Combining the temperature data with polarization, lensing, and cross-correlations [22], tightens the bound to be roughly consistent with our estimate, Eq. (34). On the other hand, UBDM with ≥ 10 −24 eV is indistinguishable from the black best-fit CDM curve. Note that this plot rescales the -axis in the upper panel by one power of ℓ compared to the usual convention, to enhance the visibility of high-ℓ features, and that the -axis begins at ℓ = 50, since the large scales are not sensitive to this particular physics.
scales for sound waves in the photon-baryon plasma, and alters the radiation driving term by changing the relative densities of matter (including UBDM) and radiation. These effects change the relative heights of the CMB acoustic peaks. An additional effect occurs in the diffusion damping (Silk damping) at larger multipoles, since the diffusion scale depends on the expansion rate during the radiation era. Due to the above mentioned effects, the CMB is sensitive to the relative contribution of Ω UBDM ( osc ). However, any fluid component with < 1/3 becomes increasingly sub-dominant to the radiation at early times (as is the case for axion-like UBDM) and so Ω UBDM decreases moving deeper into the radiation era. Because of this decrease in UBDM / , the CMB is unable to distinguish between axion-like UBDM and CDM for osc 10 −5 [29]. Plugging = 10 5 in Eq. (5) and requiring > ( osc ) gives the bound (see Fig. 1 > 2.6 × 10 −25 eV (primary CMB anisotropies) , using the same reference parameters as Eq. (6). UBDM effects on the CMB are illustrated in Fig. 3. A detailed study of these effects on the Planck CMB anisotropies constrains axion-like UBDM violating Eq. (34) (but satisfying Eq. 4) to be at most a few percent of the total DM density [20,22,29]. We have spent a considerable time deriving what will turn out to be a rather weak lower bound on . However, this bound is extremely rigorous in practice, in a way that our later bounds are not. The bound Eq. (34) relies only on linear physics, and on the extremely well understood statistics of the CMB that give us our most rigorous evidence for the existence of DM in the first place.
A complex scalar with = 1, 1/3 prior to osc increases its energy density relative to radiation at early times. The effect in the expansion rate is similar to adding additional neutrino species, which are also strongly constrained by the CMB [23].

UBDM Hints: Precision Cosmology and ALPs from the GUT Scale
The realignment production mechanism of ALPs gives the relic density Ω as a function of mass and initial field value, . Taking to be near the GUT scale, ∈ [10 15 , 10 17 ] GeV gives a DM relic density compatible with the observed value Ω ℎ 2 ≤ 0.12 for all masses 10 −18 eV. At lower masses, a sub-dominant population is predicted, with the fraction of ALP DM saturating at around 0.1%. Upcoming cosmological surveys, including lensing tomography and intensity mapping, will greatly increase the sensitivity to sub-dominant components of the DM. The CMB is a 2D probe, and the number of modes measured with a cosmic variance precision is ℓ 2 max . An intensity mapping survey is 3D, measuring in the line-of-sight redshift direction, and thus has many more modes. The combination of a next generation CMB survey like the Simons Observatory or CMB-S4 with an intensity mapping survey by the Square Kilometer Array [30] or HIRAX [31] could make significant inroads into the GUT scale predictions [25], as will next generation Lyman-forest surveys (see below) and "pulsar timing arrays" [32,33]

Schrödinger-Poisson equations
The UBDM condensate coupled to general relativity obeys the Einstein-Klein-Gordon equations, derived from variation of the relevant fundamental action. In the non-relativistic limit (Newtonian approximation), for all forms of UBDM (be they ALPs, real, or complex scalars) these equations reduce to the Schrödinger-Poisson equations (SPEs): where we are using the convention that the Newtonian potential is dimensionless, and the field has canonical mass dimension one such that the average number density is: The subtraction of the background density in the Poisson equation follows from the background-perturbation split of the Einstein equations on the Friedmann background. Equations (35)-(36) are a nonlinear Schrödinger equation for the UBDM condensate, with Gross-Pitaevski self-coupling GP , which can be computed from the relativistic self interaction potential, . The SPEs fully describe the nonlinear, non-relativistic, structure formation in most astrophysical environments at low redshifts ( osc , 1/ , 1, Φ 1), i.e. the gravitational structure of UBDM at the coherence scale. One should avoid letting the name "Schrödinger" cause confusion; these equations have nothing quantum about them: is not a probability density, and there is no measurement problem or wavefunction collapse. The SPEs are simply the non-relativistic limit of the classical field equations, valid whenever the particle number is large: they are the UBDM equivalent of Maxwell's equations.

• ? Problem 2: Derivation of the Schrödinger-Poisson equations for UBDM
Take the metric Eq. (8) in the non-relativistic limit (Φ = Ψ) on a non-expanding background ( = 1). Evaluate the d'Alembertian, Eq. (11), to first order in Ψ. Substitute the ansatz: In the sense that all classical fields can be thought of as condensates.

Solution on page 40.
An instructive change of variables on the SPEs makes use of the Madelung transformation, = √ / to write the wave function as a fluid with density and velocity = ∇ . Substitution into the SPEs yields the continuity and Euler equations: where ≡ − 1 2 2 2 The continuity and Euler equations differ from those of CDM by the presence of the so-called "quantum pressure" -a misleading term, as it is neither quantum, nor a pressure. Expanding these equations to first order in UBDM and going to Fourier space, one can to verify that they are equivalent to the fluid equation Eq. (25) for pure UBDM: in the non-relativistic and linearised limit, the quantum pressure and sound speed are equivalent. For UBDM, the SPEs replace the normal Newtonian dynamics of particle DM. Solving gravitational collapse and dynamics in generality requires methods of solution of nonlinear partial differential equations. The challenge in this system is the non-local interaction from the Newtonian potential, the wide range of scales in gravitational collapse, and the need to accurately resolve the phase of the field in low density and large cosmic voids. Common numerical methods include lattice field theory (discretizing derivatives in real space), spectral methods (numerical Fourier analysis), or finite elements (alternative real space discretizations). A public code is [34]. Particle-based hydrodynamics using Eqs. (39)- (41) is also useful on some scales, but it fails to resolve interference fringes (as can be seen from the co-ordinate singularity in when → 0) and vortex lines, which appear generically in complex fields (the fluid has ∇ × = 0). On scales larger than the UBDM de Broglie wavelength, standard Newtonian particle mechanics is accurate e.g. the public code [35]. The convergence of the SPEs to the ordinary collisionless limit of CDM on super-de Broglie scales can be shown rigorously via the Schrödinger-Vlasov correspondence [36][37][38], and is well known in the field of quantum hydrodynamics [39].
A kinetic description of the SPEs begins by writing the field using the Wigner distribution (see e.g. Ref. [40]), which describes the occupation probability of modes . This distribution function obeys a collisional Boltzmann equation, with scattering time scale [41]: where is the typical speed in the system (i.e., the virial velocity) and log Λ = log( max / min ) is the Coulomb scattering logarithm for min and max the minimum and maximum length scales in the problem, respectively. This gravitational scattering time scale governs the time over which wavelike effects cause UBDM to depart dynamically from CDM. In addition to the scattering timescale, solution of the SPEs leads to UBDM having distinctive effects on scales of order the de Broglie wavelength. There are three important consequences: 1. Transient "quasi-particle" fluctuations. 2. Formation of long-lived self-bound objects. 3. Interference fringes.
We discuss the first in Section 2.3, and the second in Section 3.1. Interference fringes are observed prominently in numerical simulations of galactic filaments composed of UBDM with ≈ 10 −22 eV [42,43], though the observational consequences are at present unclear.

Galaxies and nonlinear structure
The scale of suppression Eq. (30) can be converted into a DM halo mass by considering the average DM density in a sphere with radius of one half wavelength, = / ,eq : Halos that are significantly more massive than 0 will have the same abundance as in a CDM universe, while halos much lighter than 0 are largely absent. Our estimate for 0 from inspection of the linear equations of motion is within a factor of two of the suppression scale found in -body simulations of nonlinear cosmological structure formation: Ref. [44] finds 0 = 1.9 × 10 10 ( /10 −22 eV) −4/3 , where the different scaling with results from using the half-mode of the transfer function Eq. (32), ( 1/2 ) = 1/2, instead of the Jeans scale. The half-mode is always at < since ( ) decreases below 1/2 , and ( ) = 0. It is possible for structures to form at the half-mode, though they will have suppressed number with respect to CDM. The Jeans scale represents the absolute limit below which no structures form, and corresponds to lower mass halos. Thus using the Jeans scale gives more conservative limits on .
How can we constrain UBDM using our estimate for 0 ? In hierarchical structure formation, low mass halos form first, i.e., at high redshift. Halos with low masses can be identified at high redshift from the light emitted by the galaxies that they host, which is in the form of UV flux from stars, which in turn ionizes hot gas. An approximate relationship between UV flux and halo mass can be derived by so-called abundance matching. One assumes that there is a one-to-one mapping between UV magnitude, AB (the "AB system" for defining magnitude), and halo mass. This can be found assuming that the number of UV sources at some redshift , UV ( ), statistically matches the number of DM halos, ℎ ( ). The matching depends on the observations used to calibrate it, and monotonicity of each function. Current observations (e.g. Ref. [46]) are largely consistent with monotonicity (however, see Ref. [48]), which is consistent with all sources being in halos with mass above 0 . In this case, all the halos we observe are formed on scales far from the Jeans scale, and so the relationship between UV magnitude, , and halo mass, ℎ ( ) is as in Fig The estimate Eq. (44) agrees favourably with complete analyses of similar data [49,44,50].
Another important bound to consider is from the Lyman-forest flux power spectrum. This observable traces the matter power spectrum, ( ), on quasi-linear scales at high redshifts. It can be used to infer the existence of otherwise of the UBDM Jeans scale. Current observations see no evidence for the Jeans scale, and can be used to infer that the UBDM de Broglie wavelength must be correspondingly small. "Abundance matching" between halo mass, ℎ measured in solar masses ( ), and UV magnitude, , assuming CDM, evaluated at different redshifts, . Taken from Ref. [45]. Filled circles show the limiting magnitudes for the Hubble Ultra Deep Field observation [46], while stars are for the future James Webb Space Telescope [47]. The dotted lines represent power law extrapolation from the simulations, while the shaded region denotes the cooling limit below which galaxies cannot form efficiently.
The light from distant quasars is absorbed by neutral hydrogen (HI) along the line of sight. The differing optical thickness of dense clouds of HI leads to a "forest" of absorption features: the optical depth for the absorption traces the HI density, and (since HI clouds lie in gravitational potential wells) the total matter density including DM. A survey of cosmological quasars can then be used to estimate the matter power spectrum by correlation of the absorption feature. For example, HIRES/MIKE covers as large as max ≈ 50 ℎ Mpc −1 [51]. The data are well described by CDM with no evidence for a suppression of power, and so we can derive an approximate bound on the UBDM mass. Using Eq. (30) with the quoted max gives the bound (cf. Fig. 1): which again agrees well with the result derived from more careful analysis [56,57]. Caution is advised with all our estimates on UBDM mass bounds in this section, since they assume that the observations agree perfectly with CDM, and thus that on scales observed UBDM can be treated as such. Strong self-interactions of UBDM also change these bounds, and any other bound based on the suppression of structure formation relative to CDM. A particular example of this is an ALP with large initial field displacement. The ALP potential is ( ) = 2 2 [1 − cos( / )]. An initial displacement = / = − with small leads to large self-interactions at early times, and the field is near an unstable local maximum of the potential. This We convert from Lyman-units for in s km −1 to the more standard Mpc −1 by multiplication with 0 (1 + ). For reviews and discussions of the Lyman-forest as a probe of the matter power spectrum see Refs. [52][53][54][55].
tachyonic instability in the evolution of leads to an increase in the UBDM power spectrum relative to CDM on scales close to the Jeans scale [58]. A displacement ≈ 0.02 is sufficient to evade the bound Eq. (45) and allow ≈ 10 −22 eV to fit the Lyman-power spectrum as well as CDM, while a value ≈ 0.003 leads to a better fit than CDM [59]. The tuned values of require smaller to get the correct relic abundance than in the harmonic approximation, which could make direct detection of this type of tuned UBDM easier by increasing the matter couplings.
UBDM displays dynamics distinct from CDM on scales of order the de Broglie wavelength. A complete description of the effects of sub-de Broglie physics requires numerical simulation. However, analytical understanding is possible in varying degrees of complexity, which has largely been developed in recent years (see e.g. Refs. [54,[60][61][62][63][64]). We will give only the simplest description useful for estimates.
UBDM in a gravitational potential well has a coherence length, ∼ dB = 1/ (ℏ/ in physical units), and coherence time ∼ 1/ 2 , where is the characteristic velocity. The heuristic picture of a wave distribution with these properties is one of quasi-particles of size and lifetime . The quasiparticle mass is : where the Coulomb logarithm in the quasi-particle picture is log Λ = log( / dB ). On timescales longer than relax UBDM departs from the SHM (in the sense that the density distribution is not time-independent) due to heating and cooling. Note the similarity of the relaxation time Eq. (47) to the gravitational scattering timescale Eq. (42) in the kinetic picture if we substitute 2 =¯2. Heating and cooling on the timescale relax can be observed if a tracer population of stars with mass is present in the UBDM halo (when the gravitational potential due to DM is dominant, stars are tracer particles). For qp , heating dominates, while for qp , cooling dominates. Let's estimate qp for some systems of interest. In the solar neighbourhood,¯≈ 0.4 GeV cm −3 = 10 7 kpc −3 and ≈ 100 km s −1 ⇒ = 0.2(10 −22 eV/ ) kpc, which gives qp ≈ 7 × 10 4 (10 −22 eV/ ) 3 . In the solar neighbourhood tracers are stars with ∼ 1 , and the transition from heating to cooling occurs for UBDM mass ≈ 4 × 10 −21 eV, with lighter masses giving rise to heating. The Milky Way in fact possesses a "thick disk" of old stars [65], and this has been argued to provide evidence that in fact DM is composed of UBDM in this so-called fuzzy DM regime [54,64] (for more information, see the "Fuzzy Dark Matter Hints" box below).
On the other hand, if heating is too efficient then the disk will be destroyed completely. Demanding that the relaxation time is shorter than the age of the Universe, i.e. 10 10 years, and applying Eq. (47) we find: which agrees with more accurate modelling [64]. A very strong bound from UBDM heating can be derived by considering the existence of the old, centrally located star cluster in the ultrafaint dwarf galaxy Eridanus II. Observations [66,67] indicate that the DM density is¯= 0.15 pc −3 , and the velocity dispersion is = 6.9 +1.2 −0.9 km s −1 . For UBDM this gives qp = 3(10 −19 eV/ ) 3 , implying that heating dominates for masses 10 −19 eV. The star cluster has a half-light radius of ℎ = 13 pc, an estimated age ∼ 10 10 years, and is close to the centre of Eridanus II. Using Eq. (47), replacing with the half-light radius (since the star cluster is approximately centrally located), with , taking log Λ ∼ O (1), and demanding that the star cluster is stable on the time scale of its age, gives the bound: A potential is said to have a tachyonic region if ( ) < 0, i.e. a local maximum, and negative effective mass squared.
which again agrees very favourably with a more rigorous treatment [63]. Based on the present analysis, the bound from Eridanus II does not, however, apply for 10 −21 eV, where the fluctuation time scale becomes longer than the star cluster orbital period, and potential fluctuations become adiabatic. Another time dependent feature of UBDM halos becomes important at 10 −21 eV: the central soliton (see Section 3.1) undergoes a random walk on scales of order its own radius (which is much larger than the star cluster radius in this case) due to collisions with the quasi-particles in the halo. This again leads to star cluster disruption and could exclude ≈ 10 −22 eV from the Eridanus II star cluster stability. However, the Milky Way tidal potential may lead to sufficient tidal stripping of the quasi-particle atmosphere to quell this random walk, and leave ≈ 10 −22 eV safe from this bound [68].

• ? Problem 3: Relaxation of UBDM
The timescale for gravitational two-body relaxation (diffusion of a body's velocity caused by gravitational interaction in two-body close encounters) of particles with mass moving with velocity in a host of mass with radius can be written as [65]: Use this to derive the relaxation timescale, relax , in Eq. (47).
Solution on page 42.

UBDM Hints: Fuzzy Dark Matter
We have seen a large variety of constraints on UBDM with mass 10 −22 eV from cosmic large scale structure. We have also seen how heating in Eridanus II excludes the range 10 −21 eV 10 −19 eV, and we will see shortly that Black Hole superradiance excludes 10 −19 eV 10 −16 eV. There is only one strong bound in the range just above 10 −22 eV coming from the Lyman-forest flux power spectrum. This bound is sensitive to aspects of astrophysical modelling and, in particular, can be relaxed if the baryon temperature evolves non-monotonically or if significant ionizing photons are produced outside of galactic halos, e.g., in filaments (however, see the recent Ref. [69]). Another possible window is afforded by the Eridanus II bounds around 10 −21 eV, where the statistical modelling is uncertain, and Eridanus II can survive sandwiched between orbital resonances. If either of these bounds (Ly-or Eridanus II) can be relaxed, then there are some hints that DM may in fact be UBDM with masses between about 10 −22 eV and 10 −21 eV, the so-called fuzzy dark matter (FDM) model (cf. Fig. 1). These hints include: • The Milky Way "thick disk": FDM just outside the bound Eq. (48) can help explain the old thick disk in our galaxy [64].
• Suppressed high-galaxy formation: The redshift of reionization is known to be around reion ≈ 8. This relatively low value naturally explained by FDM, which suppresses formation of galaxies at 8. • Solitons and galactic cores: Solitons in FDM halos (see Section 3.1) may help explain cored density profiles in dwarf galaxies without baryonic feedback [42,70].
• Relic density: The relic density is naturally explained by an FDM ALP with close to the GUT scale, as expected in certain string compactifications.
Each hint provides a method to search for FDM. Furthermore, The FDM mass range corresponds to field oscillation frequencies of order one inverse month, making it challenging, but not impossible, to search for via direct detection.

Black hole superradiance
In the following we adopt different units: so-called geometric units where = = 1.
Spinning black holes (BHs) are described by the Kerr metric, which has two parameters: mass, , and dimensionless spin = / ∈ [0, 1]. In "Boyer-Linquist" coordinates the line element is: where we use spherical polar coordinates. The zero solutions of Eq. (53) define the two horizons ± : an inner Cauchy (causal) horizon at − , and the outer physical event horizon at + . The "ergoregion" is defined as radii smaller than ergo , where 00 = 0 (the co-efficient of 2 in the line element). If an object enters the ergoregion between + < < ergo , and ejects some mass which falls into the event horizon, then the object will emerge from the ergoregion with a larger energy than it went in with, and the BH will lose a small amount of energy in the form of mass and spin. This is known as the Penrose process. A wavepacket has a finite extent, and can "eject" part of itself into the BH if it passes through the ergoregion and overlaps with the event horizon. If the wave is trapped near the BH, then this process continually extracts energy from the BH, growing the wavepacket amplitude and becoming "superradiant." The process only ends when the ergoregion has shrunk small enough to remove the overlap (ultimately, the process must stop if = 0, i.e., a Schwarzschild BH with no ergoregion). Such a situation is in fact realized naturally for a massive bosonic field. Gravitational bound states trap the field near the BH, and the hydrogen-like wavefunctions overlap with the superradiant region between the ergosphere and the event horizon. The field in question must be bosonic in order that the wavepacket energy levels can continue to be filled as energy is extracted. "Black hole superradiance" (BHSR) for bosonic fields is discussed in detail in Refs. [72,73]. Consider a scalar field near a Kerr BH. Just like in the tutorial on cosmic structure above, the field obeys the Klein-Gordon equation, Eq. (10), except that now the d'Alembertian ( ) should be evaluated with the metric Eq. (51). Let us write the field as where ℓ ( ) are the spheroidal harmonics (eigenfunctions of the Laplacian on the surface of a spheroid, respecting the axial symmetry of the Kerr spacetime). To avoid confusion, we have labeled the magnetic quantum number and the azimuthal angle . The Klein-Gordon equation can then be reduced to a timeindependent Schrödinger equation for the radial eigenfunctions ℓ , with eigenvalue . The BH provides a background potential ( , ), which possesses a barrier separating the bound states from the horizon, and a potential well with size of order the boson Compton wavelength, 1/ . The system resembles a hydrogen atom with effective fine structure constant eff ≡ , where we temporarily re-instated . The existence of superradiant solutions is determined by the imaginary part of the eigenvalue , which leads to growth of the occupation number of the mode ℓ . The superradiant rate is Γ SR ∝ 4ℓ+4 eff , and numerically it is found to be maximised around eff ∼ 1. This gives an approximate criterion for BHSR: An accessible introduction to general relativity can be found in Ref. [71].
For BHSR to be effective, the superradiant time scale should be longer than any timescale of relevance for the BH, e.g. accretion. If BHSR is effective, then the BH will lose spin. Thus large observed values of will be disfavoured if a boson exists satisfying Eq. (57).
Astrophysical observations indicate the existence of BHs across a wide range of masses, from those formed by collapse of stars at the Chandrasekhar limit ≈ 1.4 , to the supermassive BHs (SMBHs) at the centres of galaxies. The spins of BHs can also be estimated, using X-ray spectroscopy of the accretion disk, or by measurement of the gravitational waveform in the inspiral phase of binary systems. Detectable spins are generally large, 0.5. Assuming that these large values would be disfavoured by a boson satisfying Eq. (57) we can estimate exclusions on UBDM. First consider the stellar BHs, and assume a full spectrum of observations from the Chandrsekhar mass to the LIGO inspiral masses ≈ 30 [74]. This excludes UBDM for: 4 × 10 −12 eV < < 8 × 10 −11 eV (stellar BHs) .
Next, consider SMBHs. The lightest currently known SMBH is in NGC4051, with mass ≈ 1.9×10 6 , while the Event Horizon Telescope has imaged the BH at the centre of M87 and determined the mass ≈ 6.5×10 9 . Again, assuming a continuous spectrum in between we can exclude the range of UBDM masses: 2 × 10 −20 eV < < 7 × 10 −17 eV (supermassive BHs) .
These estimates agree somewhat favourably with more accurate treatments of BHSR modelling and BH population statistics (e.g. Ref. [75]).
To obtain the more accurate picture, the bosonic field equations on the Kerr background should be solved numerically. The oscillation time scale of the field is ∼ 1/ . For real scalar fields the gravitational pressure oscillates with a frequency 2 , sourcing oscillations of the metric potentials on a time scale faster than the superradiant timescale. This makes brute force numerical solution challenging, but many approximation methods are available.
BHSR also works for massive spin-one and spin-two fields (which are also UBDM candidates). The superradiant timescales can be vastly different, and specific treatments are necessary. Reference [76] considers spin-one vectors which have much smaller instability rates, and thus weaker constraints. Reference [77] considers spin two fields, which possess a particular mode mimicking the spin zero case, and thus have similar constraints. A significant difference occurs for complex fields. Due to the underlying U(1) symmetry and conserved particle number, the complex vector field does not source oscillations in the metric potentials with frequency . This greatly simplifies the numerical task, and has allowed direct simulation of superradiance with these so-called Proca fields [78]. The simulations are important because they include nonlinear back-reaction of the superradiant cloud on the Kerr spacetime and demonstrate that BHSR occurs in this more realistic setting.
One known "showstopper" for BHSR is the so-called "Bosenova" caused by attractive quartic selfinteractions, which shut off the instability and prevent growth of the scalar cloud. The self-interaction term in the potential is int = 4 /4!, for some coupling constant . As the cloud grows, this term can become as large as the other terms in the energy budget. At this time, the scalar cloud collapses and superradiance is shut off. This introduces a new timescale into the problem and practically gives rise to a maximum above which the superradiance rate is subdominant to the Bosenova rate and no spin extraction can occur. Numerical simulations [79] determine the maximum cloud occupation number before Bosenova occurs [80]: where is the energy level of the occupied cloud and is the reduced Planck mass. In the last equality we assumed that the scalar potential is of the ALP form ( ) ∝ − cos( / ), giving = 2 / 2 . Using this formula for stellar mass BHs, Ref. [80] finds that BHSR is shut off for 10 13 GeV; for SMBHs this turns out to be 10 16 GeV. Any UBDM interactions can compete with superradiance, and possibly shut it off. Examples include interactions between the cloud and the Standard Model particles in the BH environment or the ALP interaction , which leads to stimulated decay of the cloud [81,82]. Of course, both "showstoppers" (Bosenova and axion-photon interactions) also predict new observables in the form of emission from BH regions for UBDMs of particular masses. Finally, we note that the superradiance phenomenon need not be limited strictly to BHs, and can occur also near stars and neutron stars [83] -even though the astrophysical uncertainties are far greater.

• ? Problem 4: Estimating superradiance properties of UBDM
A simple way to estimate the relevance of BHSR is to inspect terms in the action, where is the Ricci scalar of the Kerr background metric and , , and are the UBDM mass, (dimensionless) self-coupling, and field value, respectively. Note that it is useful to re-instat (or ) for this exercise. Assuming a suitable setup in which superradiance indeed occurs, estimate both the superradiance and Bosenova conditions i.e. Eqs (57) and (60). Note the similarity between your estimate of Bose and Eq. (60) when = 2 / 2 .
Solution on page 43.

UBDM Hints: LIGO and the QCD Axion
The exclusion estimates, Eqs. (58)-(59), assumed continuous BH distributions between the minimum and maximum values. In reality, the distributions are of course incomplete. In fact, this can serve as a discovery tool for UBDM. If light bosons with particular masses exist, then the observed BH mass and spin distribution should contain forbidden regions, and astrophysical BHs should cluster along superradiant "trajectories" in the ( , ) plane. Gravitational wave observations will, over time, provide a very complete survey of this plane.
Furthermore, superradiant clouds emit their own gravitational waves due to level transitions and annihiliation. From these effects, the LIGO observatory provides a discovery channel for UBDM with 10 −13 eV 10 −12 eV [84]. This region is disfavoured by current measurements of BH spins, but the excluded region is determined by the uncertainty on BH masses with a small number of measurements. Thus there is the possibility to make discoveries with precise measurements and greater statistics. The accessible mass region for LIGO corresponds to the QCD axion with ∼ . For proposed GW detectors in lower frequency bands corresponding to higher mass BHs (e.g., Laser Interferometer Space Antenna, LISA), discovery potential moves to lower UBDM masses.

Summary of gravitational constraints
Current constraints on UBDM mass and cosmic density from the CMB, galaxy formation, relaxation, and black hole superradiance, are combined in Fig. 5, along with a selection of forecasts for upcoming surveys. They cover an astonishing 24 orders of magnitude in mass and place sub-percent constraints on the density parameter. We caution that the limits apply strictly only to scalar UBDM with UBDM = −1 in the early Universe and negligible self-interactions, e.g., ALPs and similar cases. However, the limits apply by order-of-magnitude to  Constraints assume a real scalar with potential ( ) = 2 2 /2, see text for clarification on generalizing the bounds. CMB: cosmic microwave background [22,85], PTA: pulsar timing array [32], BHSR: black hole superradiance [75], Ly-a: lyman alpha forest [86,87]. SKA-IM: Square Kilometer Array intensity mapping [25]. Adapted from Ref. [88].
all UBDM, particularly if they come from non-relativistic effects where model dependence is less important. In addition to the effects discussed in detail, we also show projections for the measurement of pulsar timing arrays (PTA) with the Square Kilometer Array [32]. Current bounds from this technique [33] are not yet at the O (1) level for Ω UBDM , and so do not appear.

Axion compact objects
ALP UBDM can form two different types of gravitationally bound objects which are distinct from ordinary DM galactic halos. These objects, miniclusters and axion stars, are interesting phenomenologically since they are far denser than galactic halos. They can thus have observational effects as sources of enhanced DM decay and conversion, gravitational lensing, or in direct detection if they happen to pass through the Earth.

Axion stars
There exist several classes of (pseudo-)solitonic solutions to the Einstein-Klein-Gordon equations. These solutions go by many names, and have been discovered and re-discovered many times. They date back to Wheeler's idea of a "geon": a wave confined to a finite region by gravity, thus mimicking a lump of matter. Ruffinni and Bonnazola [89] found explicit "boson stars" as time-independent fixed particle number state solutions for a complex scalar field coupled to general relativity: these are true solitons, stabilised by the existence of the conserved U(1) scalar field charge. Solutions also exist for a real scalar field. However, in this case there is no conserved charge and instead the solutions have a time-dependent metric, and are known as "oscillatons" [90]. We could continue with the soliton bestiary for some time, but instead we will focus on the most well-motivated class of these objects: axion stars. First, consider the fully relativistic case. We are interested in time-dependent solutions for a scalar field coupled to general relativity. A public code is GRC [91]. Like all stars, axion stars are stabilised by a balance between attraction (gravity, and axion quartic self-interactions) and repulsion (gradient pressure, To continue the bestiary just a little further, solutions are named for all scalar fields: inflaton stars, moduli stars, Higgs stars, etc. http://www.grchombo.org/ and higher order interactions). Initial conditions are found solving the boundary value problem on the initial spacetime volume (hypersurface), and evolved forward in time to investigate their stability. The solutions are a two parameter family in mass, , and axion decay constant, , giving a "phase diagram" that can be explored numerically [92].
The structure of the axion star phase diagram is easy to understand. As the mass of the star increases, the central value of the field 0 also increases. There are two possible instabilities, and which wins depends on . For large , the self-interactions can be neglected. Now the ordinary GR lore applies: collapse to a BH at large mass. At low , the axion has strong self-interactions, and these also drive collapse. Collapse increases 0 further until higher order repulsive interactions take over and expel relativistic axions from the collapsing core in an "axion nova" [93], which occurs at critical mass nova = 10.4 / 4 , where 4 is the coefficient of quartic interactions equal to unity for a cosine potential. For small the restoring interactions become important earlier during collapse, and bring the star back to a stable configuration with only slightly lower mass than before the nova. As increases, it takes more and more of the mass of the star to contract and reach the repulsive core, thus expelling a larger mass in the nova, and reducing the mass of the stable remnant. The two types of instability are divided by a particular value of . As → ∞, oscillatons and boson stars are found to be unstable when 0 ∼ (this defines the "Kaup mass", ∼ 2 / ), while self-interactions become important when 0 ∼ , and so the boundary between the two unstable regions occurs for ∼ . A third phase boundary exists between the nova and BH regions, which simulations have found to be fractal in structure [94]. It is not clear this boundary could be reached by any astrophysical process, and so it is likely only a mathematical curiosity. The "triple point" between all three phases is found numerically to be near where is the "Arnowitt-Deser-Misner" mass [95].
Non-relativistic axion stars are far simpler to study: in the non-relativistic limit the real scalar field possesses an effective conserved particle number. In this case the solutions are simply referred to as solitons and the results apply generically to UBDM in the non-relativistic limit. Solitons are stationary waves of the form ( , ) = ( ) − , where is a dimensionless function giving the radial profile, and is the energy eigenvalue. An important property of the SPEs (see Sec. 2.2) is their scaling symmetry: where is the scale parameter (not to be confused with the quartic interaction strength in Section 2.4). The boundary value problem normalised to (0) = 1, = 1, can be solved numerically and the results are fit by eigenvalue = −0.692 2 and radial density profile: These solutions are the ground state of the SPEs. They are a balance of the nonlinear and non-local gravitational force in the Poisson equation, and the dispersive effect of the gradient energy term in the Schrödinger equation. Soliton dynamics can be studied using the numerical methods already discussed. In the limit of vanishing self-interactions, the soliton solutions are a one-parameter family given by the mass, . Thanks to the scaling symmetry, we only need to find the solution once, and then scale it using (see Problem 5). How might axion stars form in astrophysical environments? Two mechanisms are seen in simulation of the SPEs. Which occurs depends on the scale, , of the gravitational fluctuations compared to the de Broglie wavelength: • Direct collapse: ∼ dB (e.g. Ref. [42]). • Kinetic condensation: dB (e.g. Ref. [41]).
The axion potential is ( ) = 2 2 [1 − cos / ) ]. Taylor expanding this we find that the 4 term is attractive, while higher order terms alternate in sign.
Due to co-ordinate transformations, mass is not a straightforward quantity to define in general relativity (indeed, sometimes it is not defined). The Arnowitt-Deser-Misner mass is defined in the Hamiltonian formulation of general relativity, and is essentially the conserved mass measured in the infinite future.
Direct collapse leads to rapid formation of axion stars on the gravitational free-fall time, and by definition occurs in the smallest objects near to the cut-off scale of gravitational fluctuations, i.e., ∼ 0 , Eq. (43).
This mechanism leads to an axion star in the centre of all DM halos close to the cut-off scale. If all mergers of this first generation are complete up to the largest scale of halos observed, then the numerically determined relationship between the star mass, ★ , and the halo mass, , is: where the constant of proportionality can be found in Ref. [96], and depends on the definition of 0 . This relationship is believed to derive from a combination of the virial theorem, equilibrium between the soliton and its gravitationally bound "atmosphere," and universal mass growth in the merger history of solitons [97]. Slow growth of solitons by accretion leads to significant scatter in the relation. The direct collapse mechanism is particularly relevant to the formation of solitonic cores in dwarf galaxies in the FDM regime (see hint box above), and the formation of axion stars in miniclusters (discussed below). Axion stars formed by this mechanism are in virial equilibrium with their environment for < relax , and do not change appreciably in mass over such scales. The surrounding halo is a hot "atmosphere" for the star. The constant interaction with the halo causes the star to undergo radial oscillations at the normal mode frequencies [100].
Kinetic condensation gives rise to axion star formation in regions much larger than the de Broglie wavelength, for example in the solar neighbourhood for the QCD axion. The scattering timescale thermalizes the distribution function on time scales of order gr , Eq. (42), and at this time the local ground state is found in the form of an axion star which condenses spontaneously. Axion stars formed in this way continue to grow over time as they swallow up matter from the environment, with ∝ ( / gr ) . The index is to be determined numerically, and will evolve slowly in time with the wave distribution function. The growth process will eventually slow down when the star grows a gravitationally bound "atmosphere," at which point it should enter a local virial equilibrium solution close to Eq. (64). Despite progress in our understanding of the formation and growth of axion stars, at the time of writing their abundance and galactic distribution is not fully understood even in benchmark models. The problem is partly one of scale: we do not know the mass above which the relation Eq.  [42,70], or central mass excesses (tracer stars outside the soliton). See Hint box above for more details.
• Direct detection: The passage of axion stars through Earth, though rare, will greatly enhance the signal in a direct search, and could be identified using a co-ordinated network of detectors like the Global Network of Optical Magnetometers to search for Exotic physics (GNOME) and GPS.DM [101,102].
• Indirect detection: The high axion density creates a larger radio signal from decay and conversion of axions into photons (see Section 4.2). Cataclysmic signals could arise if the stars can reach the critical mass for an axion nova or stimulated decay due to interactions. • Relativistic axion stars: If dense enough axion stars can be formed, they may show up as "Exotic Compact Objects" in gravitational wave detectors [103] and multi-messenger astronomy [104].

Miniclusters
A second special class of UBDM compact objects is formed by the process of spontaneous symmetry breaking, if this occurs during the normal course of thermal evolution of the Universe (as opposed to during the initial Very recently some authors have even found a different best-fit exponent [98,99], and numerical convergence may also play a part. The issue is not yet resolved at the time of writing. conditions epoch, inflation or otherwise). The Peccei-Quinn (PQ) phase transition (see Chapter 2) occurs when the temperature of the Universe drops below approximately . Recall that we write the complex PQ field as = , and spontaneous symmetry breaking occurs when the field takes on a vacuum expectation value. The following scenario applies specifically to ALPs where the field is heavy and unstable (such that it decays at late times), while the field is initially massless, but acquires a mass hierarchically smaller than the mass of and at some time much later than the time of PQ symmetry breaking. When PQ symmetry breaking occurs, takes on a non-zero vacuum expectation value, and thus must also be specified. Since the axion field is massless at symmetry breaking, the only terms in the Lagrangian are proportional to , meaning there can be no preferred value for . The axion thus takes on a random value on essentially all scales. Imagine a pencil falling over from its point: in the absence of an external preference, the pencil falls in a random direction specified by an angle, , with the = 0 axis arbitrary.
First, consider the simpler two-dimensional case, illustrated in Fig. 6. Because the PQ field is a continuous function (as all fields must be), for any random configuration of a complex field there will be points in space around which makes a complete wrapping. At the wrapped point, the axion field is undefined (imagine shrinking the circle to a point: at the point the circle must have zero size, and takes every value at once). The only way that this can be possible is if = 0 at the wrapped point. The point in the complex field space where the radial co-ordinate is zero indeed has undefined phase. As long as the complete windings of persist, then at the centre of these windings must remain at the origin, and thus symmetry breaking cannot happen. When the potential is ( ) = (| | 2 − 2 /2) 2 this implies that the potential at the origin is (0) = 4 /4, and this is the value of the potential at the centre of a point around which wraps.
In fact, in three dimensions, cannot wrap just a single point or else the field would be discontinuous. The field must wrap continuous one-dimensional lines (either infinitely long or in closed loops) known as cosmic strings, and in this particular case as axion strings, or global strings (since the symmetry breaking is of a global U(1)). This leads to the existence of PQ strings: continuous one dimensional structures around which makes complete windings and where the radial field is pinned at = 0. The strings are a class of topological soliton: localised field configurations stabilised by the topology of the field space. The formation mechanism is known as the Kibble-Zurek mechanism, and is observed experimentally in condensed matter phase transitions with the same symmetries, for example the transition from normal fluid to superfluid helium [4,5].
What happens to the axion field? The equation of motion for the axion field in Fourier space is: where we are careful to distinguish between the field and the modefunction˜:˜= 0 does not imply = 0 as a preferred value, only that the mode is absent from its spectrum and thus gradients of on the spatial scale 1/ are small. At early times, the QCD axion mass, ( ) is vanishingly small comapred to and can be neglected. Imagine initially that all modes are populated with some amplitude (for example, the inflationary fluctuations of the PQ field), and then the field configuration far from any string is allowed to evolve. Any mode "inside the horizon" has 2 2 2 . These modes will undergo damped oscillation, and decay in amplitude. Modes larger than the horizon, 2 2 2 remain pinned to their initial value by the friction term (co-efficient of ˜) in Eq. (65) given by 3 (Hubble friction). Thus high frequency modes decay and low frequency modes remain static, smoothing the field on scales of order the horizon size. Around any string, the axion field is wound ∈ (− , ], and so we have O (1) variation of the field on horizon size patches around the string. String formation is sketched in Fig. 6. Furthermore, numerical simulations indicate that string dynamics enter into a scaling solution with O (1) strings per horizon volume.
Strings decay when the axion mass becomes relevant to the mode evolution. Recall first that during cosmological expansion temperature, , always decreases as time, , increases. Second, recall that the axion mass and the Hubble rate are both decreasing functions of . The mass term (co-efficient of˜) in Eq. (65) is This is generic for complex fields in three spatial dimensions. It is a topological property. Complex fields have symmetry group U(1) of rotations in the complex plane, i.e. loops. The mapping of U(1) onto R 3 (Euclidean 3-space) is expressed by the first homotopy group 1 (R 3 ). This group is not the empty set, i.e. it is non-trivial, which can be seen by noting that R 3 is the universal  comparable to the friction term when ( ) ≈ ( ). The axion mass term defines a preferred value for the field, = 0, which is exactly why the PQ mechanism solves the strong-CP problem. Eq. (65) is just a damped harmonic oscillator, and so the mode functions for all < (not just the short wavelength modes inside the horizon) will begin to oscillate when ( ) ≈ ( ), defining the special temperature osc . Now, everywhere, the axion field is making harmonic oscillations about zero. There are thus no longer regions around which it makes a complete and continuous winding. The axion field everywhere has average value zero (but importantly of course, nonzero variance and energy density). This means that the radial mode is no longer required to take the value = 0 along the strings. The axion field is everywhere defined, the radial field is not pinned, and it can undergo symmetry breaking at the string locations, i.e. the strings decay (or "unwind").
At this time the axion field has a well specified distribution: in every horizon-size patch, it has O (1) fluctuations, while on larger scales it is uncorrelated. The power spectrum, ( ), is flat (white noise) for ( osc ) ( osc ) = ( osc ) and cut off by the Jeans scale for ( osc ) ( osc ). The normalisation of the power spectrum is fixed by the variance, which should match the variance of the uniform distribution for , 2 = 2 /3. Just prior to ( osc ) the axion equation of state is ≈ −1, and so the fluctuations in do not source any curvature perturbations in the metric: this is what is meant by the term isocurvature. It is this particular power spectrum (white noise isocurvature, with O (1) variance, truncated at the horizon size at osc ), which gives rise to the structures known as axion miniclusters as a remnant of string decay [105]. Figure 7 shows a snapshot from numerical simulation of the axion field after string decay [106], and miniclusters are located using a threshold based on spherical collapse under gravity [107]. cover of 3 , the 3-torus, and 1 ( ) = Z . This last can be seen since one cannot shrink circles on tori to points continuously, and there are distinct circles wrapping . Note that the mode function˜decaying to zero does not imply a preference for the axion field to move to zero: in the massless limit the shift symmetry prevents any such preference.
Intuitively this can be understood because if = −1 exactly then this is equivalent to a cosmological constant, which is constant in space and time, and thus cannot source a spatially varying curvature. z = 3 × 10 4 z = 5 × 10 4 z = 7 × 10 4 z = 9 × 10 4 Fig. 7 Initial conditions for minicluster formation. After string decay, the axion field has large density perturbations, which subsequently collapse into the objects known as miniclusters. The figure shows a small patch of results from the numerical simulations of Ref. [106], which used lattice field theory methods to solve the equations of motion for the complex Peccei-Quinn field in the absence of gravity. The large hierarchies involved necessitate further approximations, and different simulation methods are not currently in precise agreement for the spectrum of perturbations extrapolated to physical values of the particle masses and couplings. Miniclusters are identified at different redshifts using a threshold derived from spherical collapse under gravity [107].
The mass scale of miniclusters is the same as the mass scale of ordinary axion halos: it is fixed by the Jeans scale when field oscillations begin, Eq. (43), but now for very different values of the reference parameters. However, because of the large amplitude of these isocurvature fluctuations (variance of order unity), axion miniclusters begin to collapse earlier than ordinary DM halos and there is significant nonlinear structure formation before matter-radiation equality. The density of any DM halo is related to the background density at the time when it first collapses and leaves the Hubble flow. Thus, miniclusters are denser than ordinary halos, and may survive repeated mergers up to the present day. Let us consider the phenomenology of these low mass, dense objects.
First, we need the minicluster mass which we estimate from the number density of axions within the comoving cosmological horizon at the time when field oscillations begin. We find osc in the usual way, setting Contrast this to the evolution of large scale inhomogeneities seen in the CMB, and the inflationary "adiabatic" mode. These perturbations have small amplitude, and a red slope in the power spectrum leading to smaller amplitude fluctuations on small scales. These fluctuations undergo logarithmic growth at early times in the radiation era, which is important to seed galaxy formation and is one of the pieces of evidence for DM discussed in Chapter 1. However, due to the existence of a free-streaming scale/Jeans scale for particle DM models, there is, in general, no nonlinear structure formation during the radiation era. For WIMPs with 1 TeV, structure formation in the adiabatic mode begins at ≈ 500 with the formation of Earth mass, 10 −6 , halos [108].
( osc ) = 1 ( osc ). Now we need to calculate the horizon volume. Take the volume of rotation of the spherical wave with comoving wavenumber osc = ( osc ) ( osc ) over one-half wavelength: This defines the Hubble volume . Alternative definitions are the cubic volume, ( ) = ( ) −3 , and the spherical volume of of one half wavelength ( ) = 4 3 ( / ) 3 . The expected minicluster mass is simply MC = ( osc osc )Ω crit . We compute MC ( osc ) using the Friedmann equation to fix osc , 3 2 2 = ( 2 /30) ★ 4 , and conservation of entropy to write ( ) ∝ −1 ★ (normalised using the CMB measurement of eq ). The result is well fit by: where S( ) is an activation function that accounts for the behaviour of ★ and ★, in the Standard Model (dominantly the quark-hadron phase transition) and has only been roughly fit here using the results for ★ from Ref. [109]. Note that this expression is valid for any ALP with minicluster-like initial conditions. From the QCD axion ,QCD ( ) dependence based on lattice QCD in Ref. [110], osc ( ,QCD ) is well fit by: over the range of interest (broadly 10 −5 eV ≤ ,QCD ≤ 10 −3 eV) for the relic density in this scenario (see also Now we would like to know the minicluster density profile. Kolb and Tkachev [111] wrote down the equation of motion for spherical collapse of an isolated top-hat density profile, with initial overdensity , in an expanding Universe dominated by radiation. The perturbation first grows in size as the Universe expands. It then turns around, and collapses. Spherical collapse formally leads to a density singularity. However, in real collapses below the threshold for BH formation, aspherical perturbations lead to virialization (equilibrium between average kinetic and gravitational potential energy expressed by the virial theorem, see, e.g., Ref. [112]) and the collapsed object becomes self-supported with a finite average density. virialization occurs when the radius of the perturbation is half of the turn-around radius. Using this information, one can compute the overdensity of the spherical system at the time of virialization. A numerical solution of the ordinary differential equation for spherical collapse gives the final average overdensity: Assuming the minicluster has a constant density, the radius is then calculated to be 1 is simply a constant of proportionality to account for ambiguity defining osc , and its later use in analytical formulae for the evolution of the energy density. Our earlier choices, e.g. Eq. (4), set 1 = 1, physically assuming structure formation begs when the de Broglie wavelength is equal to the Hubble length, −1 . Many authors choose 1 = 3 when estimating the relic density. This ambiguity in the use of osc leads to significant uncertainty in analytical minicluster mass estimates, which can only be resolved by fitting results of numerical simulations.
For the more standard case during matter domination, which applies to ordinary DM halos, see, e.g., Ref. [113]. In the standard case, the equations can be solved analytically, leading to the well-known result that the virial overdensity /¯≈ 200, independent of .
Using¯e q = 2 × 10 14 kpc −3 To make more use of this result, we need to know the minicluster radial profile, ( ). Collapse of isolated density perturbations is self-similar, and leads to a power law density profile with no preferred scale. Miniclusters are not isolated, and the formation proceeds much like ordinary DM halos, leading to Navarro-Frenk-White (NFW) profiles [114][115][116]. It is then natural to associate the radius MC with the NFW scale radius. If the initial distribution of could be measured, one would know the mass and size distribution of miniclusters.
In reality, the problem of miniclusters is far more complex than the simple story given here. Firstly, miniclusters do not have one fixed mass. Structure formation always proceeds hierarchically, and there is a mass function of miniclusters. This can be computed numerically via N-body simulation, or semi-analytically from the initial power spectrum [117,118,116,107]. The mass function takes on a power law spreading over many orders of magnitude around MC . Secondly, the distribution of in initial conditions is not uniquely defined. Numerical thresholding using the spherical collapse results allows some progress to be made [107] but the results still require calibration to N-body simulation. Unfortunately, N-body simulations cannot currently resolve the scale radius on all relevant scales, and are not large enough to capture the rarest, densest, and thus most phenomenologically interesting miniclusters.
Finally, just like other UBDM halos, when the effects of the gradient energy (the UBDM de Broglie wavelength) are included, miniclusters have been shown to form central axion stars [119]. Axion stars in miniclusters follow approximately the same core-halo mass relation, Eq. (64), as an ordinary halo. For the QCD axion, the resulting axion stars are on very different mass scales for the UBDM particle mass and the halo mass than the reference FDM values in Eq. (64).
The minicluster power spectrum, mass function, size function, and central axion stars, can all be used to constrain the QCD axion and ALPs in the post-inflation PQ symmetry breaking scenario. Some examples include: • Microlensing and femtolensing [120][121][122] (see Problem 5 below).
• The CMB isocurvature power spectrum and large scale structure [126,127]. to 10 −6 [128]. The Einstein radius for gravitational microlensing, is = 2[ ★ (1 − ) ] 1/2 , where is the distance from the observer to the lens, is the distance from the observer to the source, and = / . Compare the axion star radius to with = 1/2 and = 770 kpc, the distance to M31. What UBDM particle masses could be probed by HSC lensing due to axion stars? Now consider the mass-radius relation for miniclusters with initial overdensity , Eq. (72) and the minicluster mass relation Eq. (68). What range of the ( osc , ) parameter space can be probed by microlensing?
Solution on page 44.

Stellar and supernova energy loss
In this section we consider only constraints on axion-like couplings, i.e., pseudoscalar, anomalous, or shift symmetric (see Chapter 2). Analogous bounds can of course be derived for scalar and dilaton-like couplings. Axions are pseudoscalars, and their couplings to fermions depend on the orientation of the spin, while couplings of scalar particles do not. The spin dependence can lead to suppression of interactions, since macroscopic bodies are not in general strongly polarized. Being unsuppressed by spin effects, scalar constraints are often stronger. More details on some of the calculations are given in Chapter 5.
First and foremost, it is extremely important to remember that the constraints and effects we discuss in this section apply independently of whether the axion is (a large fraction of) the DM. The axions considered here are produced from Standard Model particles in stellar plasmas. They interact only very weakly and have a long mean-free path inside the plasma. Thus, stars are effectively transparent to axions, and the axions escape, allowing an additional cooling channel for the star. This changes the evolution of stars: in simple terms, it alters the progression of stars along the Hertzsprung-Russell (HR) diagram of stellar luminosity versus temperature. The relationship between the mass, age, and temperature of stars is thus different than in the Standard Model. Stellar physics is generally very well understood in terms of Standard Model physics alone, and can be simulated using a code such as [129], which can be modified to include axion-induced cooling [130]. For more details see Refs [131,10,132].
Stars, including the Sun, can produce axions by the Primakoff process: photons inside the star convert into axions in the ambient magnetic and electric fields of the particles (electrons and nuclear ions) in the plasma. The rate for this process is: where is the photon energy, is the temperature, and is the screening length. In the Debye-Hückel approximation we have where is the free electron density, and is the density of the -th nuclear ion of charge . In a neutral medium with = = 0 the Primakoff rate goes to zero, since there is no background field to facilitate conversion. Photon energies are distributed thermally, and the temperature varies with stellar radius. Kinematically, we require ≥ to produce an axion. Typical stellar interior temperatures are in the keV range, which gives the typical energy of the emitted axions, and approximately the maximum axion mass where this cooling channel is allowed. The luminosity in axions needs to be computed for a given stellar model. Applying this to the Sun gives = 1.85 × 10 −3 ( /10 −10 GeV −1 ) 2 . It is this solar luminosity in axions that helioscope experiments try to detect (see Chapter 5). We can derive a crude bound by demanding that the solar axion luminosity must be less than unity, since the evolution of the Sun is well described by emission dominantly in photons, i.e. ,Sun = 1 by definition. Thus: The bound in Eq. (75) can be improved by considering the statistics of populations of stars. The best understood case is for horizontal branch (HB) stars in globular clusters. Stars in globular clusters are all of a similar age, and differ in their masses. The distribution of the stars gives an HR diagram that can be compared to models. The observable is the ratio of HB stars to red giant branch (RGB) stars, , determined by placing stars on a colour-magnitude diagram. In the Standard Model, this ratio is a function of the primordial helium abundance, He , stellar mass, and metallicity. Globular clusters are old systems, with ages in the range of 10 billion years. This gives a small range of available stellar masses and metallicities, which have a negligible effect on . The value of He can determined observationally by measurement of extragalactic H II regions which gives He = 0.2449 ± 0.0040 [133], which is consistent with the predictions of standard BBN and the CMB measurement of the baryon abundance [3]. Reference [134] reports a measured average value of = 1.38 from 39 globular clusters, consistent with the Standard Model prediction.
With a given stellar evolution model it is possible to compute the effect of the axion-photon coupling, , and the axion-electron coupling, , on . At present there are two different models in the literature for the functional dependence, and each is presented in Ref. [135]. The specific forms are not enlightening, so we simply quote the bounds (derived in Ref. [136]). In both cases the additional cooling channel lowers compared to the Standard Model prediction leading to a degeneracy in the combined constraints, with the maximum value of one coupling only allowed when the other is strictly zero. Setting one coupling to zero, and fixing He = 0.25, the individual bounds are: where the number in brackets refers to the bound using the alternative model for , which we see introduces an O (1) shift in the bound on . For the constraints in the combined parameter space, see Refs. [135,136]. Supernova SN1987A provides an important bound on the axion nuclear couplings, and . During the core collapse process, a proto-neutron star is formed, the gravitational field of which traps neutrinos and causes them to be emitted over a delayed period of time. This model explains the duration of the burst of two dozen observed neutrinos coincident with SN1987A. Axion emission due to nuclear bremsstrahlung: would compete with neutrino emission and cool SN1987A too rapidly, shortening the neutrino burst, unless the total energy loss rate from either axion-nuclear process obeys the bound 10 19 erg g −1 s −1 = 7.2 × 10 −18 eV. For the first process, the cooling rate per unit mass is [132] where is the mass density of the supernova, is the nucleon number density, is the axion angular frequency, and is the spin density structure function, which accounts for the fact that axions couple only to the nuclear spins. The coupling is the nucleon coupling weighted as 2 = 2 + 2 , where and are the neutron and proton abundances, estimated as = 0.3 and = 1 − = 0.7 at the relevant epoch in the supernova. The spin density structure function is non-trivial to compute, but can be estimated in various approximations. The last equality in Eq. (80) defines the dimensionless function from the integral of , which is estimated to be of order unity, and allows a simple estimate of the bound given the supernova internal temperature ≈ 30 MeV. A recent analysis found the more accurate bound [137] 2 < 1.3 × 10 −9 GeV −1 (SN1987A neutrino burst) , a factor of approximately four weaker than the estimate with = 1. The same analysis found an O (1) effect from the modelling of supernova temperature and density profiles. The bound from SN1987A on / is particularly important for the QCD axion, since this couplings is always present, and so the bound can be cast as a model-independent constraint on the QCD axion mass. We use that 2 = 0.066 ⇒ = 0.257, leading to: ,QCD 0.06 0.257 eV (SN1987A neutrino burst) .
For the second process is approximated by: where are the reactant number densities, is the average photon energy, is the supernova mass density, and is the thermally averaged cross section. The nuclear number density and supernova mass density are known, and the other parameters are fixed in terms of the internal temperature, . To estimate the bound on , Ref. [138] approximates the cross section as = 2 2 , leading to:

UBDM Hints: Anomalous White Dwarf Cooling
White dwarfs (WDs) are stellar remnants whose electron-degenerate cores are supported by Fermi pressure against gravitational collapse. Their internal densities are relatively high as their masses are typically comparable to the mass of the Sun (∼ 0.6 ), while their radii are of the order of the Earth's radius [139]. They cannot replenish their internal energy, and therefore continuously cool down over time.
The evolution of WDs can be altered by introducing additional cooling channels. These can be provided by weakly-interacting particles that efficiently carry away energy after being created in, and escaping from, the WD's core. A useful observable to infer the resulting, additional cooling rate is the so-called period increase in variable WDs. These are WDs that periodically change in brightness over time as they pulsate due to non-radial excitations, called "gravity modes" (see e.g. Ref. [140]), with potentially multiple pulsation periods associated with different co-existing sub-modes. For cooling WDs, the periods of their pulsations, Π, tend to increase over time with a rate Π = dΠ/d that can approximately be calculated via where and are the internal temperature and radius of the WD, respectively [141]. The change in radius can usually be neglected for the observed, low-luminosity dwarfs [142]. Axions and ALPs induce an energy loss that is proportional to 2 since axion-electron interactions dominate in the high-density, electron-degenerate interior of the WDs [143,144]. The resulting decrease in temperature, and therefore the additional contribution to the period increase in Eq. (85), is hence also proportional to 2 . Measurements of an anomalous period change can thus be used to estimate the associated axion-electron coupling. From the 250 known variable WDs [145], this has so far only been done for G117-B15A [146][147][148][149][150], R548 [148,151], L19-2 [152], and PG 1351+489 [153]. The reason for this small fraction is that measuring the period change is very difficult: while the periods for the WDs listed here are of the order of a few minutes, their (inherently dimensionless) period changes, Π, are less than about 10 −13 in magnitude.
The left panel of Fig. 8 shows how the measured period increase in G117-B15A (blue line and shading) compared to theoretical prediction from simulating WD evolution with and without axions (red data points). To illustrate the dependence on the axion-electron coupling, we show the theoretical prediction as a function of 2 .
The right panel of Fig. 8 shows the one-dimensional profile likelihoods for the four WDs listed above. Combining these likelihoods hints at an additional cooling channel corresponding to an axion-electron coupling of a few times 10 −13 at more than 3 confidence level [135,155,136]. In addition to difficulties of observing the period change, there are a number of uncertainties involved in the modelling of WDs and their pulsations. Multiple challenges in quantifying the statistical and systematic uncertainties in WD modelling remain, such as the modelling of the transition from the main sequence to the WD phase. More details on WD modelling can be found in Refs. [150][151][152][153]. It is therefore not yet clearly established whether the cooling hints are due to systematics or indicate the presence of new physics-be it in form of a weakly-interacting particle or a completely different astrophysical cooling channel.
A recent, more general review of pulsating WDs can be found in Ref. [140]. Apart from the period increase discussed here, the WD luminosity function can also be used to probe the evolution of WDs and be seen as a hint for ALPs [156,157].

Axion-photon conversion
In the presence of a magnetic field, axions convert into photons, and vice versa, by the Primakoff and inverse Primakoff process (see Chapters 2,4,and 5). This leads to constraints on the axion-photon coupling from any astrophysical environment penetrated by a magnetic field. In the following, we briefly mention some important instances.
Axions produced during supernova SN1987A can escape from the supernova event. These axions are subsequently converted back into visible photons in the form of gamma rays by the magnetic field of the Milky Way. This process would have led to a gamma ray burst coincident with SN1987A, which was not observed [158]. This places constraints on the axion mass and coupling at 95% CL [159]: The bound gets rapidly worse at higher masses due to loss of coherence of the axion field on the scale of the galactic magnetic field and the resulting reduced photon fluence. The bound has an O (1) dependence on the precise model of the galactic magnetic field. Axion-photon conversion also occurs in the intergalactic medium, and leads to modulation of the X-ray spectra of active galactic nuclei (AGN) and quasars (see e.g. Refs [160,161]). The modulation can be modelled statistically with a stochastic model for cluster magnetic fields. The strongest bound arises from the observation of a single source, NGC1275, by the Chandra satellite, which observed no modulations and sets the 3 limit [162,163] < 6-8 × 10 −13 GeV −1 (for 10 −12 eV) .
Still further bounds can be derived from axion-photon conversion on cosmological scales. The conversion of CMB photons by Mpc scale primordial magnetic fields leads to CMB spectral distortions (i.e. departure from a blackbody spectrum) [164,165]. Since the cosmic background explorer (COBE) satellite determined the CMB to be the most perfect black body in the Universe [166], any departures from perfection caused by axion-photon conversion are strongly constrained. On the other hand, the origin and spectrum of large scale, primordial cosmic magnetic fields is highly uncertain (e.g. Ref. [167]). Thus, bounds are given relative to the amplitude of the magnetic field power spectrum averaged on cosmic length scales, = √︁ 2 , as 10 −14 GeV −1 1 nG (for 10 −12 eV) .
These bounds can be improved by up to two orders of magnitude by future CMB spectral measurements.

Solutions to Problems Problem 1: Background evolution of UBDM
The continuity equation, / = −3( + 1) / , can be integrated such that ∝ −3( +1) . Note that this solution is also valid for the case = −1 i.e. when is constant. Substituting this into the Friedmann equation, we find where 1 is some constant. Integrating the equation above for, we obtain Consequently, we find via = / that To change variables from conformal to physical time, note that This yields the following field equation for¯: The solutions of this equation can be expressed in terms of Bessel functions of first ( ) and second kind ( ) with = ( − 1)/2( + 1) and ≡ where 3 and 4 are constants that depend on the initial conditions. While these may be used to derive expressions for the pressure and density, we note that the asymptotic behaviour of Eq. (95) may be derived more generally. Note that, for early times, → 0, we have ∝ 1/ → ∞ and the corresponding term will dominate: More precisely, this is the regime where ( ) .
where 5 and 6 are constants. In typical applications, e.g., a radiation-dominated universe with = 1/3, we have¯= 5 + 6 −1/2 . The second term is divergent for → 0 and we also see that¯ 5 for which are two arguments usually used for ignoring the second term and saying that the background field is simply constant at early times,¯= 5 . On the other hand, for late times ( ) , Eq. (95) is solved by a WKB-like solution, and without loss of generality we take¯∝ −3/2 cos( ) and ¯∝ −3/2 [ cos( ) − sin( )]. Note that also implies that the oscillation periods ∼ 1/ 1/ ∼ such that we can assume that and do not change much over each integration 2 ¯2 = + ⇒ 2 ¯2 = + ≡ ( eff + 1) .
Using sin( ) cos( ) = 0 and sin 2 ( ) = cos 2 ( ) = 1/2, we find that ¯2 ∝ ( / ) 2 /2 + 1/2 , (99) = 2 ¯2 /2 + ∝ 2 /4 + 2 /4 + 2 /4 , For the potential ( ) = 4 , Eq. (95) now becomes For early times, there is no change to the previous argument since the potential is irrelevant. For later times, the shape of Eq. (102) implies that we cannot rely on the WKB-like solutions anymore, and instead follow the general approach presented in Refs [1,2]. We repeat this derivation for a symmetric potential , i.e. (−¯) = (¯), and (¯) =¯( even). At the maximumˆthe total energy density is given by the maximum potential value, (ˆ) = (ˆ) ≡ˆ. Since the oscillations at late times are very rapid, the energy density of the field does not change much over one oscillation period , and ≈ˆ. The definition of then implies that and we find for the oscillation period that which in turn implies that This can be checked by substituting this in Eq. (95) and noting that the non-vanishing terms are all small if . Note that, for the quartic potential (¯) =¯4, the solutions for¯can be expressed in terms of so-called Jacobi elliptic functions times an oscillating function [2], which can be used to check the general solution presented here explicitly.
i.e. eff = 1/3 for = 4, which corresponds to the equation-of-state parameter of radiation. We see that when Hubble friction dominates, i.e., typically at sufficiently early times, any scalar particles generically behave as dark energy with eff = −1. Later on, however, the (dominant terms of) the potential determine the behaviour of the scalar field oscillations, which need not be that of dark matter (pressureless dust) and dark matter bounds hence do not apply.

Problem 2: Derivation of the Schrödinger-Poisson equations for UBDM
First, let us write down the Klein-Gordon equation, where the D'Alembertian is defined as and the potential is given by With Ψ = Φ and = 1, the metric is given by To first order this gives: Thus, the D'Alembertian to first order is given by: We can rewrite the Klein-Gordon equation by multiplying with −(1 + 2Ψ) (since this term goes to −1 in the non-relativistic limit, the result remains unchanged save for an overall minus sign), which reduces the number of terms we need to consider later on: Let us take the ansatz for and write: Since terms for − are the complex conjugate, the terms in front of + need to vanish in order for the Klein-Gordon equation to be fulfilled. Thus, one needs to consider only the terms which go with a + oscillation (or − , respectively). This gives for Eq. (114): Now, let us take the non-relativistic limit and consider either the limits given in the exercise or take → ∞. We calculated in natural units, where = 1. Remember when revoking the 's that Ψ → Ψ/ 2 , that the time derivatives gain a factor of 1/ and that -for each within the brackets -→ . Hence, the terms in red vanish and one is left with: Considering the complex conjugate, this equals Eq. (35) with GP = − . To calculate the energy density, , to leading order we recognize that = − 0 0 . The stress energy tensor is given by: and, thus: Suppose, for simplicity, that we are at the equator ( = /2) and that we just inside the ergosphere at which is also just the Schwarzschild radius of the black hole. We then find For typical field values ∼ in the extreme environment surrounding the black hole, we find that where we ignored the O (1) numerical factor. Note that we might have used dimensional analysis (e.g. that ∼ 1/ S ∼ 1/2 ) instead of Eq. 124 to arrive at a similar result.
To compute the Bosenova condition, the self-coupling of the UBDM particles, , needs to be large enough to play a role. We anticipate this to happen when the corresponding term in the action, 4 /4!, becomes of the same order as the potential term, 2 2 /2. To facilitate this comparison, first note that the total energy density of the UBDM cloud can be equated to its total mass, cloud , divided by its volume, cloud , which are given by where is the number of UBDM particles. By equating 2 2 /2 = cloud / cloud , we find that where the last step made use of Eq. (127). To compare to Eq. (60), we use that = 2 / 2 for an axion and eliminate via the SR condition (127) which is very close to Eq. (60) for the first energy level ( = 1).
The mass of the soliton is then given by The rescaling relates to → 2ˆ. Then, since sol ( ) = | | 2 , Similarly, sol ( ) = 4ˆs ol (ˆ). We can therefore rescale the profile to 4ˆs ol ( ) We can define a scale radius = allowing us to write the density profile as sol ( ) = The Einstein radius for such a lens is given by, Rearranging our mass-radius relation we see For the object to lens like a point-mass, we require that < , therefore which can be rearranged to This upper limit is maximized by being sensitive to masses as small as possible. Setting sol to the smallest mass detectable by HSC < 1.4 × 10 −10 eV .
To calculate the range of the ( osc , ) prameter space, we neglect the activation function in Eq. (68) (since ( ) ∼ O (1)). Requiring again that MC < E and substituting our MC mass into the equation for the Einstein radius we find the region range of ( osc , ) parameter space can be probed by microlensing to be