11institutetext: David J. E. Marsh 22institutetext: Department of Physics,
King’s College London, London, United Kingdom,
22email<EMAIL_ADDRESS>33institutetext: Sebastian Hoof
44institutetext: Institut für Astrophysik, Georg-August-Universität Göttingen,
Germany,
44email<EMAIL_ADDRESS>
# Astrophysical Searches and Constraints on Ultralight Bosonic Dark Matter

David J. E. Marsh and Sebastian Hoof

###### Abstract

Starting from the evidence that dark matter (DM) 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.

##  1 Astrophysical search channels

Astrophysics and cosmology, as outlined in Chapter 1, give convincing evidence
that dark matter (DM) exists in the form of new particles beyond the Standard
Model of particle physics. The space of possible theories in Chapter 2, even
for the subclass of ultralight bosonic DM (UBDM) models considered in this
book, is vast. Beyond the basic fact of the existence of DM, astrophysics can
be used to reign in this vast theoretical parameter space, with a view to
direct detection and measurement of model parameters.

The most basic astrophysical route to constrain UBDM is via the relic density.
There are three channels for UBDM production:

  1. 1.

Coherent field oscillations.

    1. (a)

Vacuum realigment.

    2. (b)

Topological defect decay.

  2. 2.

Thermal production.

  3. 3.

Non-thermal production by direct decay.

Without going into the specifics (see Ref. Marsh:2015xka), 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 2013JCAP...10..020A; Aghanim:2018eyx.

In channel 1a (vacuum realignment) the UBDM relic density is a function of two
parameters, $(m,\phi_{i})$, where $\phi_{i}$ 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}\text{\,}\mathrm{eV}$,${10}^{-1}\text{\,}\mathrm{eV}$)$ for
$\phi_{i}\leq M_{pl}$.111$M_{pl}=1/\sqrt{8\pi G_{N}}$ is the _reduced_ Planck
mass, related to the mass scale given in the “Units and Conversions” section
by the factor of $\sqrt{8\pi}$ coming from Einstein’s equation in general
relativity. For an axion-like particle (ALP) the fundamental parameter from
theory is $f_{a}$: the scale of spontaneous symmetry breaking, also called the
axion decay constant. The parameter $\theta_{i}$, defined via
$\phi_{i}\equiv\theta_{i}f_{a}$, 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, $\phi\rightarrow\phi+\text{constant}$,
and thus $\theta_{i}$ has no preferred value and can be considered a free
random variable (although very small values, or values very close to $\pi$ are
considered fine tuned). Because $\theta_{i}$ is undetermined there is a wide
range of allowed values for the fundamental parameters $(m,f_{a})$ consistent
with the required relic density. In particular, in this channel large values
of the decay constant at the grand unified scale
($\sim${10}^{16}\text{\,}\mathrm{GeV}$$) or the reduced Planck scale
($\sim${10}^{18}\text{\,}\mathrm{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” 1976JPhA....9.1387K; 1985Natur.317..505Z 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
$\mathbb{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
Hiramatsu:2012gg; Klaer:2017ond; Gorghetto:2018myk. 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
$f_{a}\lesssim${10}^{12}\text{\,}\mathrm{GeV}$$ Armengaud:2019uso.

The production mechanism channel 1b works for $f_{a}<T_{\text{max}}$, where
$T_{\text{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. $T_{\text{max}}$ is
bounded from above due to observational constraints on the theory of
inflation. In particular $H_{I}$, the inflationary Hubble rate, is bounded
from above by the fact that tensor-type CMB anisotropies have relative
amplitude $r\lesssim 0.1$ compared to scalar-type perturbations leading to the
constraint $H_{I}\lesssim${10}^{14}\text{\,}\mathrm{GeV}$$. $H_{I}$ sets the
temperature of the Universe during inflation to be the Gibbons-Hawking
temperature, $T_{\text{GH}}=H_{I}/2\pi$. The maximum thermalization
temperature could actually be larger than this, which can easily be seen from
the Friedmann equation during radiation domination,
$3H^{2}M_{pl}^{2}=\pi^{2}g_{\star}T^{4}/30$, where the quantity $g_{\star}$
counts the effective number of relativistic degrees of freedom
1990eaun.book.....K:

$g_{\star}=\frac{7}{8}\sum_{i\in\text{fermions}}g_{i}\left(\frac{T_{i}}{T}\right)^{4}+\sum_{i\in\text{bosons}}g_{i}\left(\frac{T_{i}}{T}\right)^{4}\,,$
(1)

where $g_{i}$ is the degrees of freedom of species $i$ (e.g. two polarizations
for the photon) and $T_{i}$ is the temperature of species $i$, and $T$ is the
photon bath temperature. The value of $g_{\star}$ at very high temperatures is
bounded from below by the Standard Model contribution,
$g_{\star,\text{SM}}=106.75$. $H$ monotonically decreases, and so
$H_{\text{max}}=H_{I}$. If reheating after inflation is instantaneous and 100%
efficient, we find an upper bound for
$T_{\text{max}}\lesssim$8\text{\times}{10}^{15}\text{\,}\mathrm{GeV}$$. ALPs
with values of $f_{a}$ larger than this upper bound on $T_{\text{max}}$ cannot
be produced by mechanism 1b, and must be produced by mechanism 1a. The
observational lower bound on $T_{\text{max}}$ arises from demanding successful
Big Bang nucleosythesis, $T_{\text{max}}\gtrsim$1\text{\,}\mathrm{MeV}$$. For
values of $f_{a}$ in this very large range of allowed $T_{\text{max}}$ values,
it is not determined whether ALPs are produced by mechanism 1a or 1b, either
being possible depending on the model of inflation and reheating.

There are various astrophysical search channels we can use to constrain UBDM:

  1. 1.

Gravitational probes.

  2. 2.

“Indirect detection.”

    1. (a)

Production of UBDM (e.g. in stars or from background radiation).

    2. (b)

Decay/conversion of existing UBDM.

Gravitational probes are the most general form of constraints on UBDM and give
us powerful bounds on the key parameters of mass and density (both cosmic and
local), which are important for the design of direct DM searches. Indirect
detection depends on the UBDM interactions with ordinary matter: null results
provide baseline constraints on couplings to which laboratory searches are
compared, and anomalous results give hints for promising regions of parameter
space to search.

In this Chapter, unless stated otherwise, we use natural units where
$\hbar=c=k_{B}=1$ and express all quantities in electronvolts (eV). We use the
Einstein summation convention for repeated indices. Roman indices $i$, $j$,
etc. run from 1 to 3, while greek indices $\mu$, $\nu$, etc. run from 0 to 3,
with zero labelling the time-like direction. In relativity, we distinguish
covariant (lower) and contravariant (upper) indices, with the metric being
responsible for raising and lowering: $x_{\mu}=g_{\mu\nu}x^{\nu}$.

##  2 Gravitational probes of UBDM

The goal of this section is to assess the validity of UBDM as a model of DM.
Since all current observations are consistent with cold dark matter (CDM,
defined as a pressureless fluid), the bounds we estimate on the UBDM mass $m$
can be thought of as answering the question: “is UBDM observationally
equivalent to CDM?” The answer to this question depends on the observable and
leads to lower bounds on $m$ (and upper bounds on the UBDM density if we allow
for multi-component DM). In order to derive our bounds we must specify the
ways in which UBDM is not equivalent to CDM. These differences further suggest
astrophysical phenomena that could distinguish between UBDM and CDM in the
future, possibly providing evidence for one model over the other.

###  2.1 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
Dodelson:2003ft. Consider a flat, homogeneous, and isotropic spacetime
described by the Friedmann-Robertson-Walker metric:

$\displaystyle g={\rm diag}[-1,a(t)^{2},a(t)^{2},a(t)^{2}]\,.$ (2)

The scale factor is $a(t)$, which obeys Friedmann’s equation for the Hubble
rate $H(t)=\dot{a}/a$:

$\displaystyle H(t)^{2}=\frac{8\pi G_{N}}{3}\bar{\rho}\,,$ (3)

where $\bar{\rho}$ is the total, spatially averaged, energy density. $\rho$ 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 $\vec{v}_{H}=Hr\vec{\hat{r}}$, where
$r$ and $\vec{\hat{r}}$ are the distance and direction from the observer to
the object, respectively. We begin with a Newtonian approximation to cosmology
(see e.g. Ref. mukhanov). Consider an observer at the origin, and a single
particle of UBDM on the Hubble flow. The UBDM de Broglie wavelength is
$\lambda_{H}=1/(mv)=1/(mHr)$, which gives the radial position uncertainty,
$\Delta r$. A net gravitational force in the positive direction along the line
of centres between the observer and the UBDM requires $\Delta r\lesssim
r\Rightarrow r\gtrsim(mH)^{-1/2}$, which defines a critical separation
$r_{\text{crit}}=(mH)^{-1/2}$. On average, UBDM separations larger than
$r_{\text{crit}}$ undergo gravitational clustering, and those smaller than it
do not.

Figure 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
$f_{a}\in[${10}^{6}\text{\,}\mathrm{GeV}$,${10}^{18}\text{\,}\mathrm{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, $g_{\star,S}$. The
quantity $g_{\star,S}$ gives the number of effective relativistic degrees of
freedom contributing to the entropy density; $g_{\star,S}$ takes the same form
as Eq. (1) with the fourth powers replaced by cubes (see, e.g., Ref.
1990eaun.book.....K, Chapter 3).

The cosmological horizon size is approximated by the Hubble length scale
$R_{H}=H^{-1}$. In order for UBDM to have any inhomogeneous gravitational
effect within this radius requires $r_{\text{crit}}<R_{H}$. We show the
cosmological evolution of $H=R_{H}^{-1}$, and the related comoving Hubble
radius $(aH)^{-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.

Evaluating the Hubble length scale today, and using that
$H_{0}=100\,h\,$\mathrm{km}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{Mpc}}^{-1}$=$2.13\text{\times}{10}^{-33}\text{\,}\mathrm{eV}$\times
h$ (where $h$ is the dimensionless Hubble parameter with approximate observed
value $h\sim 0.7$), we arrive at our first bound on the UBDM mass:

$\displaystyle
m>$1.5\text{\times}{10}^{-33}\text{\,}\mathrm{eV}$\times\left(\frac{h}{0.7}\right)\,\quad\text{(size
of the observable universe)}.$ (4)

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.222Very 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 $V(\phi)$, and are known as “quintessence” Copeland:2006wr.
An ultralight bosonic field with a mass less than the bound Eq. (4) is one
such very simple model, with $V(\phi)=m^{2}\phi^{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 Clifton:2011jh.

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 Dodelson:2003ft:

$\displaystyle
H(z)=H_{0}E(z)=H_{0}\,\sqrt{\Omega_{\Lambda}+\Omega_{m}\left[(1+z)^{3}+\frac{(1+z)^{4}}{1+z_{\text{eq}}}\right]}\,,$
(5)

where the second equality defines the energy function $E(z)$. The quantities
$\Omega_{m}$ and $\Omega_{\Lambda}$ are the density parameters of matter and
the cosmological constant, defined as the density divided by the critical
density, i.e. $\Omega_{i}=\bar{\rho}_{i}/\rho_{\text{crit}}$ and
$\rho_{\text{crit}}=3M_{pl}^{2}H_{0}^{2}$. 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, $z_{\text{eq}}$.
The epoch of matter-radiation equality can be found via the relative
redshifting of matter and radiation components:
$\rho_{m}(1+z_{\text{eq}})^{3}=\rho_{r}(1+z_{\text{eq}})^{4}$, with the
density parameters defined today. CMB observations fix $z_{\text{eq}}\approx
3390$, and it is thus slightly earlier in cosmic history than decoupling,
$z_{\rm dec}\approx 1100$.

The _baryon acoustic oscillations_ (BAO, see Sec. 1.1) observed in the CMB and
galaxy surveys like the _Sloan Digital Sky Survey_ 2019ApJS..240...23A 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 $H(z_{\text{eq}})$ we arrive
at the tighter bound (cf. Fig. 1):

$\displaystyle
m>$1.6\text{\times}{10}^{-28}\text{\,}\mathrm{eV}$\times\left(\frac{h}{0.676}\right)\left(\frac{\Omega_{m}}{0.311}\right)^{1/2}\left(\frac{1+z_{\text{eq}}}{3390}\right)^{3/2}$
$\displaystyle\text{(matter-radiation equality)}\,,$ (6)

where we have neglected the small contribution of $\Omega_{\Lambda}$ at
equality, and taken reference parameters from the CMB+BAO combination in Ref.
Aghanim:2018eyx.333Using 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 $\Lambda$CDM).

The matter-radiation equality bound, Eq. (6), is the UBDM equivalent of saying
that DM is not “hot” Primack:2001ib: 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.

###### Example 1

Tutorial: The Growth of Cosmic Structure The challenge in cosmological
perturbation theory bertschinger1995 is to compute the _transfer function_ ,
$T_{X}(t,k)$ for the mode evolution of each cosmological species $X$ (baryons,
photons, neutrinos, dark matter) with Fourier wavenumber $k$, which fully
specifies linear evolution of cosmological fields from Gaussian initial
conditions. That is:

$\displaystyle\zeta_{X}(k,t)=\zeta_{X,i}(k)T_{X}(t,k)\xi_{X}\,,$ (7)

where $\zeta_{X,i}(k)$ is the initial condition of the field and $\xi_{X}$ is
a Gaussian random field defining the initial correlation functions of the
field $\zeta_{X}$.

The codes camb camb and class class are the standards for numerical
computation for CDM (and many other things), while axionCAMB
Hlozek:2014lca444Available at https://github.com/dgrin1/axionCAMB. can be used
for UBDM that is a real scalar field with the self-interaction potential
approximated by $V(\phi)=m^{2}\phi^{2}/2$. This tutorial gives a brief
overview of the most relevant aspects of cosmological perturbation theory for
UBDM constraints.

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 $t$. 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 $\zeta$ of interest are the components of the energy
momentum tensor, written as $T^{0}_{\,\,0}=-(\bar{\rho}+\delta\rho)$,
$T^{i}_{\,\,j}=(\bar{P}+\delta P)\delta^{i}_{j}+\Sigma^{i}_{j}$,
$ik^{i}T^{0}_{\,\,i}=(\bar{\rho}+\bar{P})\theta$, which defines the energy
density, $\rho$, pressure, $P$, and heat flux, $\theta=\nabla\cdot\vec{v}$,
and we assume anisotropic stresses $\Sigma^{i}_{j}$ vanish. This gives the
fields $\delta_{X}=\delta\rho_{X}/\bar{\rho}_{X}$, and $\theta_{X}$, while
pressure is typically described it terms of a sound speed, $c_{s}^{2}=\delta
P/\delta\rho$.

Next, perturb the metric Eq. (2), and switch to conformal time, $\tau$, via
${\rm d}t=a{\rm d}\tau$. The Newtonian gauge considers only scalar metric
perturbations:

$\displaystyle g=a^{2}{\rm diag}[-(1+2\Psi),\,1-2\Phi,\,1-2\Phi,\,1-2\Phi]\,.$
(8)

The potential $\Phi$ is the usual Newtonian potential, and $\Psi$ 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:

$\displaystyle G_{\mu\nu}=8\pi G_{N}T_{\mu\nu}\,,$ (9)

where $G_{\mu\nu}$ 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
$V(\phi)$ is:

$\displaystyle\Box\phi-\partial_{\phi}V=0\,,$ (10)

where the d’Alembertian ($\Box$) is:

$\displaystyle\Box=\frac{1}{\sqrt{-g}}\partial_{\mu}\sqrt{-g}g^{\mu\nu}\partial_{\nu}\,,$
(11)

and where $g$ and $g^{\mu\nu}$ are the metric determinant and the inverse of
the metric, respectively. Setting $V=\frac{1}{2}m^{2}\phi^{2}$ for simplicity,
this leads to the equations of motion for the UBDM background field,
$\bar{\phi}$ and fluctuation mode $\delta\phi_{k}$:

$\displaystyle\bar{\phi}^{\prime\prime}+2\mathcal{H}\bar{\phi}^{\prime}+a^{2}m^{2}\bar{\phi}$
$\displaystyle=0\,,$ (12)
$\displaystyle\delta\phi_{k}^{\prime\prime}+2\mathcal{H}\delta\phi_{k}^{\prime}+(m^{2}a^{2}\delta\phi_{k}+k^{2})\delta\phi_{k}$
$\displaystyle=(\Psi^{\prime}+3\Phi^{\prime})\bar{\phi}^{\prime}-2m^{2}a^{2}\Psi\bar{\phi}\,,$
(13)

where primes denote derivatives with respect to conformal time, and
$\mathcal{H}=a^{\prime}/a=aH$. For the UBDM field, we find $T^{\mu\nu}=\delta
S/(\delta g_{\mu\nu})$ by variation of the action with respect to the metric
tensor, giving:

$T^{\mu\nu}=g^{\mu\alpha}\partial_{\alpha}\phi\partial^{\nu}\phi-g^{\mu\nu}\left[\frac{1}{2}g^{\alpha\beta}\partial_{\alpha}\phi\partial_{\beta}\phi+V(\phi)\right]\,.$
(14)

Working to first order in the metric perturbations and $\delta\phi$, and with
potential $V=m^{2}\phi^{2}/2$ the components are:

$\displaystyle\bar{\rho}$
$\displaystyle=\frac{1}{2}a^{-2}(\bar{\phi}^{\prime})^{2}+\frac{1}{2}m^{2}\bar{\phi}^{2}\,,$
(15) $\displaystyle\bar{P}$
$\displaystyle=\frac{1}{2}a^{-2}(\bar{\phi}^{\prime})^{2}-\frac{1}{2}m^{2}\bar{\phi}^{2}\,,$
(16) $\displaystyle\delta\rho$
$\displaystyle=a^{-2}[\bar{\phi}^{\prime}\delta\phi_{k}^{\prime}-\Psi(\bar{\phi}^{\prime})^{2}]+m^{2}\bar{\phi}\delta\phi_{k}\,,$
(17) $\displaystyle\delta P$
$\displaystyle=a^{-2}[\bar{\phi}^{\prime}\delta\phi_{k}^{\prime}-\Psi(\bar{\phi}^{\prime})^{2}]-m^{2}\bar{\phi}\delta\phi_{k}\,,$
(18) $\displaystyle(\bar{\rho}+\bar{P})\,\theta$
$\displaystyle=a^{-2}ik^{2}\bar{\phi}^{\prime}\delta\phi_{k}\,.$ (19)

###### Question 1

Problem 1: Background evolution of UBDM Assuming a single-fluid universe with
constant equation of state $w$ satisfying $\dot{\rho}=-3H(1+w)\rho$, first
solve Friedmann’s equation Eq. (3) for $a(t)$, and thus $H(t)$. Then change
variables in Eq. (12) to physical time ${\rm d}t=a{\rm d}\tau$. Substituting
your solution for $H(t)$, derive the solution for $\bar{\phi}(t)$ (you may use
exact functions or asymptotic methods). Given that the energy density and
pressure of UBDM are $\bar{\rho}=\frac{1}{2}\dot{\phi}^{2}+V(\phi)$ and
$\bar{P}=\frac{1}{2}\dot{\phi}^{2}-V(\phi)$, derive the behaviour of the
equation of state for UBDM, $w_{\text{UBDM}}=\bar{P}/\bar{\rho}$. What is the
asymptotic value of $w_{\text{UBDM}}$ for $m\ll H$, and $\langle w\rangle$ for
$m\gg H$ (brackets denote period average)? Repeat this exercise for a
$\lambda\phi^{4}$ potential. Comment on the results for $w_{\text{UBDM}}$ in
relation to the approximate UBDM mass bounds above.

Solution on page 6.

CDM is defined as a collisionless and uncoupled fluid, $w_{c}=c_{c}^{2}=0$.
Baryons have a sound speed, $c_{b}^{2}\neq 0$ (computed from the evolution of
the baryon temperature), equation of state $w_{b}=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 ($l$) 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 $l_{\text{max}}$). Truncating this Boltzmann hierarchy at the velocity
moment, the photons resemble a fluid with $w=c_{s}^{2}=1/3$, collisionally
coupled to the baryons. We consider perturbations to the energy density
$\delta_{X}=\delta\rho_{X}/\bar{\rho}_{X}$ and heat flux $\theta_{X}$, defined
via $\bar{\rho}_{X}(1+w_{X})\theta_{X}=ik^{j}\delta(T^{0}{}_{j})_{X}$, where
$(T^{\mu}{}_{\nu})_{X}$ is the $X$ energy momentum tensor.

Let us now consider a number of limits of the full equations of motion, which
can be found in Ref. bertschinger1995. 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
$\theta_{\gamma}=\theta_{b}$. Considering only sub-horizon modes ($k\gg aH$),
and using the Poisson equation and the $ii$ pressure component of the Einstein
equation Eq. (9), the photon fluid at early times obeys the equation of
motion:

$\displaystyle\delta_{\gamma}^{\prime\prime}+\left(c_{s,\gamma}^{2}k^{2}-\frac{16\pi}{3}Ga^{2}\rho_{\gamma}\right)\delta_{\gamma}=4\pi
Ga^{2}\sum_{i}(1+c_{i}^{2})\rho_{i}\delta_{i}\,,$ (20)

where the photon sound speed is $c_{s,\gamma}=1/\sqrt{3}$ (speed of pressure
perturbations in a gas of photons in thermodynamic equilibrium). At very early
times all $\rho_{i}$ in the driving term on the right hand side can be
neglected. Then this equation has sound wave solutions for $k>(16\pi
Ga^{2}\rho_{\gamma})^{1/2}=\sqrt{6}aH$. 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
$z\approx 1100$) 1990eaun.book.....K, 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:

$\displaystyle r_{s}=\int_{0}^{t_{0}}\frac{{\rm
d}t}{a}c_{s,b}\approx\frac{1}{\sqrt{3}}\int^{t_{\rm dec}}_{0}\frac{{\rm
d}t}{a}\,,$ (21)

where $c_{s,b}$ is the baryon sound speed in the plasma, $t_{0}$ is the time
today, and $t_{\rm rec}$ is the time at recombination when $c_{s,b}$ drops
rapidly from $c_{s,\gamma}$ 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:555This 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. Lewis:2006fu; Hlozek:2017zzf.

$\displaystyle\frac{\delta
T}{T}=\int_{0}^{\tau_{0}}\left[\dot{\mu}\left(\Phi+\frac{\delta_{\gamma}}{4}+\vec{\hat{n}}\cdot\vec{v}_{b}+2\dot{\Phi}\right)\right]e^{-\mu}{\rm
d}\tau\,,$ (22)

where $\mu$ is the Thomson scattering opacity, $\vec{v}_{b}$ is the baryon
velocity, $\vec{\hat{n}}$ 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 $z_{\rm dec}\approx 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 $z<z_{\rm eq}\approx 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
$H(a_{\text{osc}})=m$ (this defines the value of the scale factor
$a_{\text{osc}}$ when the background UBDM field, $\bar{\phi}$, begins to
undergo coherent oscillations, see Problem 1). Prior to this time the UBDM is
relativistic and perturbations cannot grow.666For an axion-like potential, the
equation of state is $w=-1$ prior to $a_{\text{osc}}$. For a scalar field with
potential $V=m^{2}\phi^{2}/2+\lambda\phi^{4}$, the equation of state is
$w=1/3$ at early times for large $\phi$ initial conditions. For a complex
scalar, the early time equation of state is $w=1$ due to the conserved charge
and Goldstone mode 2014PhRvD..89h3536L. In each case, perturbations are
suppressed relative to pressureless CDM. For $H\ll m$ the UBDM perturbations
Eq. (13) can be approximated as a fluid with sound speed
2009PhLB..680....1H:777This 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 Hlozek:2014lca.

$\displaystyle
c_{\text{UBDM}}^{2}=\frac{k^{2}/4m^{2}a^{2}}{1+k^{2}/4m^{2}a^{2}}\,.$ (23)

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 $k\gg aH$, super-Compton $k\ll m$ limit, the
CDM+baryon and UBDM fluids obey the coupled equations of motion:

$\displaystyle\ddot{\delta}_{c+b}+2H\dot{\delta}_{c+b}=-\frac{k^{2}}{a^{2}}\Phi\,,$
(24)
$\displaystyle\ddot{\delta}_{\text{UBDM}}+2H\dot{\delta}_{\text{UBDM}}+\frac{k^{4}}{4m^{2}a^{4}}\delta_{\text{UBDM}}=-\frac{k^{2}}{a^{2}}\Phi\,,$
(25) $\displaystyle\frac{k^{2}}{a^{2}}\Phi=4\pi
G_{N}(\bar{\rho}_{c+b}\delta_{c+b}+\bar{\rho}_{\text{UBDM}}\,\delta_{\text{UBDM}})\,.$
(26)

Setting $\bar{\rho}_{\text{UBDM}}=0$ and substituting the Poisson equation Eq.
(26) into Eq. (24) gives the solution
$\delta_{c+b}=A_{+}(k)a+A_{-}(k)a^{-3/2}$. The growing mode initial conditions
sets $A_{-}(k)=0$, and the inflationary initial conditions and matter transfer
function fix $A_{+}(k)$. Due to the zero pressure and sound speed of CDM, all
the $k$-dependence in the solution is fixed by the initial conditions, and the
dynamics are scale invariant.

Now consider a UBDM-dominated universe by taking
$\bar{\rho}_{c+b}=\bar{\rho}_{b}\ll\bar{\rho}_{\text{UBDM}}$ (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
$\delta_{\text{UBDM}}$, which drives growth of $\delta_{\text{UBDM}}$, while
the positive contribution from the sound speed term leads to acoustic
oscillations. The sign of the term proportional to $\delta_{\text{UBDM}}$
depends on $k$ and as such different modes evolve differently. That is, we
find Eq. (25) exhibits a _Jeans scale_ , $k_{J}$, separating growing/decaying
and oscillating modes. The exact solution for pure UBDM is
$\delta_{\text{UBDM}}=A_{+}(k)D_{+}(k,a)+A_{-}(k)D_{-}(k,a)$, where the growth
functions are:

$\displaystyle
D_{+}(k,a)=\frac{3\sqrt{a}}{\tilde{k}^{2}}\sin\left(\frac{\tilde{k}^{2}}{\sqrt{a}}\right)+\left[\frac{3a}{\tilde{k}^{4}}-1\right]\cos\left(\frac{\tilde{k}^{2}}{\sqrt{a}}\right)\,,$
(27) $\displaystyle
D_{-}(k,a)=\left[\frac{3a}{\tilde{k}^{4}}-1\right]\sin\left(\frac{\tilde{k}^{2}}{\sqrt{a}}\right)-\frac{3\sqrt{a}}{\tilde{k}^{2}}\cos\left(\frac{\tilde{k}^{2}}{\sqrt{a}}\right)\,.$
(28) $\displaystyle\tilde{k}=k/\sqrt{mH_{0}}$ (29)

Consider the evolution of three wavenumbers in the pure UBDM case: the horizon
size, $k_{\star}=aH$; the Jeans scale, $k_{J}=a\sqrt{Hm}$; and the Compton
scale, $k_{c}=ma$. The Compton scale defines relativistic modes where
$c_{\text{UBDM}}^{2}=1$; $k_{c}$ increases with time, and more modes become
non-relativistic. If $k_{\star}<k_{c}$, 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 $a_{\text{osc}}$, which
defines the special mode $k_{m}$, the horizon size when the UBDM background
becomes non-relativistic. All $k<k_{m}$ evolve similarly to CDM. All $k>k_{m}$
have suppressed growth.

Figure 2: Evolution of scales for linear perturbations with $m=10^{-26}\text{
eV}$. The Jeans scale, Compton scale, and horizon scale, all intersect at
$a_{\text{osc}}$ when the field begins to oscillate. This defines the scale of
power suppression as the comoving horizon size at this time,
$k_{m}=a_{\text{osc}}H_{\text{osc}}=R_{H}(a_{\rm osc.})^{-1}$. Due to the slow
evolution of $k_{J}$ with $a$, 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. Bauer:2020zsj.

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
$c_{\text{UBDM}}\approx k^{2}/4m^{2}a^{2}$, substituting the Poisson equation,
and solving for $k_{J}$ where the effective mass term in the oscillator
equation for the overdensity vanishes, we find:

$\displaystyle
k_{J,\text{eq}}=9.0\left(\frac{3390}{1+z_{\text{eq}}}\right)^{1/4}\left(\frac{\Omega_{\text{UBDM}}}{0.12}\right)^{1/4}\left(\frac{m}{10^{-22}\text{
eV}}\right)^{1/2}\text{ Mpc}^{-1}\,.$ (30)

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 $m\rightarrow\infty$.
In practice, they are equivalent as long as $k_{J}$ does not play a role in
any observation.

An observable related to the matter clustering is the _matter power spectrum_
defined by
$\langle\delta_{m}(\vec{k}_{1})\delta_{m}(\vec{k}_{2})\rangle=(2\pi)^{3}\delta_{D}(\vec{k}_{1}-\vec{k}_{2})P(k)$,
where $\delta_{m}$ is the total matter (baryon+CDM+UBDM+neutrino) overdensity,
and $\delta_{D}$ is the Dirac delta distribution. The presence of the sound
speed and consequent Jeans scale for UBDM leads to a suppression of $P(k)$
relative to CDM at large wavenumbers. A fit for the relative suppression in
$P(k)$ for UBDM with $V(\phi)=m^{2}\phi^{2}/2$ versus CDM is hu2000:

$\displaystyle P_{\text{UBDM}}(k)=T_{\text{UBDM}}(k)^{2}P_{\rm CDM}\,,$ (31)
$\displaystyle T_{\text{UBDM}}(k)=\frac{\cos x_{J}^{3}(k)}{1+x_{J}^{8}(k)}\,,$
(32) $\displaystyle x_{J}(k)=1.61\left(\frac{m_{a}}{10^{-22}\text{
eV}}\right)^{1/18}\frac{k}{k_{J,\text{eq}}}\,.$ (33)

For the mixed CDM-UBDM system, the behaviour of $P(k)$ can also be derived
amendola2005; 2010PhRvD..82j3528M: perturbations with $k>k_{m}$ experience a
finite amplitude suppression which increases with the ratio
$\Omega_{\text{UBDM}}/\Omega_{m}$.

End of Tutorial

As we have just seen in the above tutorial, two effects distinguish UBDM from
other ingredients in the $\Lambda$CDM model: (1) the background expansion
rate, $H(z)$, driven by the transition in the equation of state
$w_{\text{UBDM}}$ at the epoch $a_{\text{osc}}$, and (2) the growth of
perturbations, driven by the gradient energy in the Klein-Gordon equation and
manifested as an effective sound speed, $c^{2}_{\text{UBDM}}$.

Depending on the value of $a_{\text{osc}}$, the change in $H(z)$ 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 $\ell\approx 200$. This can be compensated by a change in
the value of the Hubble constant, $H_{0}$. After such a compensation there is
a residual _integrated Sachs-Wolfe effect_ which differs from $\Lambda$CDM. If
$w\neq 0$ in the post-recombination Universe, then the gravitational potential
$\dot{\Phi}\neq 0$ into Eq. (22). Due to the fact that the equation of state
$w_{\text{UBDM}}\neq 0,-1$ (the two available equations of state in
$\Lambda$CDM), the evolution of $\Phi$ is different in the presence of a small
contribution of UBDM, and the shape of the $\ell<200$ CMB multipoles is very
sensitive to the value of $\Omega_{\text{UBDM}}$.888This is one of the ways
the CMB is used to constrain the equation of state of dark energy.

Now consider UBDM satisfying the bound Eq. (6). The change in the expansion
rate compared to $\Lambda$CDM now occurs during the radiation dominated epoch.
The horizon size at the time $a_{\text{osc}}$ was smaller than one degree on
the sky, corresponding to multipoles $\ell>200$, i.e. the higher acoustic
peaks. UBDM changes the distance 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 $\Omega_{\text{UBDM}}(a_{\text{osc}})$. However, any fluid
component with $w<1/3$ becomes increasingly sub-dominant to the radiation at
early times (as is the case for axion-like UBDM) and so $\Omega_{\text{UBDM}}$
decreases moving deeper into the radiation era.999A complex scalar with
$w=1,1/3$ prior to $a_{\text{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 2014PhRvD..89h3536L. Because of this decrease in
$\rho_{\text{UBDM}}/\rho_{\gamma}$, the CMB is unable to distinguish between
axion-like UBDM and CDM for $a_{\text{osc}}\lesssim 10^{-5}$ Poulin:2018dzj.
Plugging $z=10^{5}$ in Eq. (5) and requiring $m>H(z_{\text{osc}})$ gives the
bound (see Fig. 1

$m>$2.6\text{\times}{10}^{-25}\text{\,}\mathrm{eV}$\quad\text{(primary CMB
anisotropies)}\,,$ (34)

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 Hlozek:2014lca;
Hlozek:2017zzf; Poulin:2018dzj. We have spent a considerable time deriving
what will turn out to be a rather weak lower bound on $m$. 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.

Figure 3: UBDM effects on the CMB temperature power spectrum. UBDM changes the
expansion rate compared to CDM in the early radiation dominated epoch,
$z\gtrsim 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 $m\leq 10^{-26}\text{ eV}$.
Combining the temperature data with polarization, lensing, and cross-
correlations Hlozek:2017zzf, tightens the bound to be roughly consistent with
our estimate, Eq. (34). On the other hand, UBDM with $m\geq 10^{-24}\text{
eV}$ is indistinguishable from the black best-fit CDM curve. Note that this
plot rescales the $y$-axis in the upper panel by one power of $\ell$ compared
to the usual convention, to enhance the visibility of high-$\ell$ features,
and that the $x$-axis begins at $\ell=50$, since the large scales are not
sensitive to this particular physics.

###### Example 2

UBDM Hints: Precision Cosmology and ALPs from the GUT Scale The realignment
production mechanism of ALPs gives the relic density $\Omega_{a}$ as a
function of mass and initial field value, $\phi_{i}$. Taking $\phi_{i}$ to be
near the GUT scale, $\phi_{i}\in[10^{15},10^{17}]\text{ GeV}$ gives a DM relic
density compatible with the observed value $\Omega_{d}h^{2}\leq 0.12$ for all
masses $m\lesssim 10^{-18}\text{ 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 $\ell_{\text{max}}^{2}$. 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_ SKA_RedBook or HIRAX HIRAX could make significant
inroads into the GUT scale predictions Bauer:2020zsj, as will next generation
Lyman-$\alpha$ forest surveys (see below) and “pulsar timing arrays”
Khmelnitsky:2013lxt; Porayko:2018sfa. These forecasted opportunities are shown
as open regions in Fig. 5.

###  2.2 Schrödinger-Poisson equations

The UBDM condensate101010In the sense that all classical fields can be thought
of as condensates. 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):

$\displaystyle
i\dot{\psi}+\frac{\nabla^{2}}{2m}\psi-m\Phi\psi+\frac{\lambda_{\rm
GP}}{m}|\psi|^{2}\psi=0\,,$ (35) $\displaystyle\nabla^{2}\Phi=4\pi
Gm^{2}\left(|\psi|^{2}-\int{\rm d}^{3}x|\psi|^{2}\right)\,,$ (36)

where we are using the convention that the Newtonian potential is
dimensionless, and the field $\psi$ has canonical mass dimension one such that
the average number density is:

$\displaystyle\bar{n}=m\int{\rm d}^{3}x|\psi|^{2}\,.$ (37)

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 $\lambda_{\rm GP}$, which can
be computed from the relativistic self interaction potential, $V$. The SPEs
fully describe the nonlinear, non-relativistic, structure formation in most
astrophysical environments at low redshifts ($a\gg a_{\text{osc}}$, $L\gg
1/m$, $v\ll 1$, $\Phi\ll 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: $\psi$ 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.

###### Question 2

Problem 2: Derivation of the Schrödinger-Poisson equations for UBDM Take the
metric Eq. (8) in the non-relativistic limit ($\Phi=\Psi$) on a non-expanding
background ($a=1$). Evaluate the d’Alembertian, Eq. (11), to first order in
$\Psi$. Substitute the ansatz:

$\phi=\frac{1}{m\sqrt{2}}\left(\psi e^{imt}+\psi^{*}e^{-imt}\right)\,,$ (38)

into the Klein-Gordon equation with potential
$V(\phi)=m^{2}\phi^{2}/2+\lambda\phi^{4}$. In the Wentzel-Kramers-Brillouin
(WKB) limit, $\dot{\psi}/(m\psi)\ll 1$, and making the non-relativistic
approximation $k/m\ll 1$, $\dot{\Psi}/m\ll 1$, show that the complex field
amplitude $\psi$ obeys the Schrödinger equation Eq. (35). Now take the general
form of the stress energy tensor, Eq. (14), and show that in the same limits
$\rho=|\psi|^{2}$ at leading order, and hence that the Poisson equation Eq.
(36) is obeyed for the overdensity $\delta\rho$.

Solution on page 7.

An instructive change of variables on the SPEs makes use of the _Madelung
transformation_ , $\psi=\sqrt{\rho}e^{i\theta}/m$ to write the wave function
as a fluid with density $\rho$ and velocity $\vec{v}=\nabla\theta$.
Substitution into the SPEs yields the continuity and Euler equations:

$\displaystyle\dot{\delta}_{\text{UBDM}}+a^{-1}\vec{v}_{\text{UBDM}}\cdot\nabla\delta_{\text{UBDM}}=-a^{-1}(1+\delta_{\text{UBDM}})\nabla\cdot\vec{v}_{\text{UBDM}}\,,$
(39)
$\displaystyle\dot{\vec{v}}_{\text{UBDM}}+a^{-1}\left(\vec{v}_{\text{UBDM}}\cdot\nabla\right)\vec{v}_{\text{UBDM}}=-a^{-1}\nabla(\Phi+Q)-H\vec{v}_{\text{UBDM}}\,,$
(40) $\displaystyle\text{where
}Q\equiv-\frac{1}{2m^{2}a^{2}}\frac{\nabla^{2}\sqrt{1+\delta_{\text{UBDM}}}}{\sqrt{1+\delta_{\text{UBDM}}}}\,.$
(41)

The continuity and Euler equations differ from those of CDM by the presence of
the so-called “quantum pressure” $Q$ – a misleading term, as it is neither
quantum, nor a pressure. Expanding these equations to first order in
$\delta_{\text{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 $\psi$ 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 pyultralight
Edwards:2018ccc. 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 $Q$ when $\rho\rightarrow 0$) and
vortex lines, which appear generically in complex fields (the fluid has
$\nabla\times\vec{v}=0$). On scales larger than the UBDM de Broglie
wavelength, standard Newtonian particle mechanics is accurate e.g. the public
code gadget Springel:2005mi. 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 Widrow&Kaiser1993;
Uhlemann:2014npa; Mocz:2018ium, and is well known in the field of quantum
hydrodynamics wyatt_trajectories.

A kinetic description of the SPEs begins by writing the field $\psi$ using the
Wigner distribution (see e.g. Ref. ballentine_book), which describes the
occupation probability of modes $k$. This distribution function obeys a
collisional Boltzmann equation, with scattering time scale Levkov:2018kau:

$\displaystyle\tau_{\rm
gr}\approx\frac{\sqrt{2}}{12\pi^{3}}\frac{mv^{6}}{G^{2}\bar{n}^{2}\log\Lambda}\,,$
(42)

where $v$ is the typical speed in the system (i.e., the virial velocity) and
$\log\Lambda=\log(r_{\text{max}}/r_{\text{min}})$ is the Coulomb scattering
logarithm for $r_{\text{min}}$ and $r_{\text{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. 1.

Transient “quasi-particle” fluctuations.

  2. 2.

Formation of long-lived self-bound objects.

  3. 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 $m\approx 10^{-22}\text{ eV}$
Schive:2014dra; Mocz:2019emo, though the observational consequences are at
present unclear.

###  2.3 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, $R_{J}=\pi/k_{J,\text{eq}}$:

$\frac{M_{0}}{M_{\odot}}=$5.9\text{\times}{10}^{9}$\left(\frac{\Omega_{\text{UBDM}}h^{2}}{0.12}\right)^{1/4}\left(\frac{h}{0.676}\right)\left(\frac{1+z_{\text{eq}}}{3390}\right)^{3/4}\left(\frac{m}{${10}^{-22}\text{\,}\mathrm{eV}$}\right)^{-3/2}.$
(43)

Halos that are significantly more massive than $M_{0}$ will have the same
abundance as in a CDM universe, while halos much lighter than $M_{0}$ are
largely absent. Our estimate for $M_{0}$ from inspection of the linear
equations of motion is within a factor of two of the suppression scale found
in $N$-body simulations of nonlinear cosmological structure formation: Ref.
Schive:2015kza finds
$M_{0}=$1.9\text{\times}{10}^{10}\text{\,}\,$M_{\odot}\,(m/${10}^{-22}\text{\,}\mathrm{eV}$)^{-4/3}$,
where the different scaling with $m$ results from using the half-mode of the
transfer function Eq. (32), $T(k_{1/2})=1/2$, instead of the Jeans scale. The
half-mode is always at $k<k_{J}$ since $T(k)$ decreases below $k_{1/2}$, and
$T(k_{J})=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 $m$.

Figure 4: “Abundance matching” between halo mass, $M_{h}$ measured in solar
masses ($M_{\odot}$), and UV magnitude, $M_{AB}$, assuming CDM, evaluated at
different redshifts, $z$. Taken from Ref. Schultz:2014eia. Filled circles show
the limiting magnitudes for the _Hubble Ultra Deep Field_ observation
Bouwens:2014fua, while stars are for the future _James Webb Space Telescope_
Windhorst:2005as. The dotted lines represent power law extrapolation from the
simulations, while the shaded region denotes the cooling limit below which
galaxies cannot form efficiently.

How can we constrain UBDM using our estimate for $M_{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, $M_{\text{AB}}$ (the “AB system”
for defining magnitude), and halo mass. This can be found assuming that the
number of UV sources at some redshift $z$, $n_{\text{UV}}(z)$, statistically
matches the number of DM halos, $n_{h}(z)$. The matching depends on the
observations used to calibrate it, and monotonicity of each function. Current
observations (e.g. Ref. Bouwens:2014fua) are largely consistent with
monotonicity (however, see Ref. Leung:2018evj), which is consistent with all
sources being in halos with mass above $M_{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, $M_{AB}$, and halo mass, $M_{h}(M_{AB})$ is as in Fig. 4
(computed from simulations of CDM with no Jeans scale). The limiting magnitude
of the _Hubble Ultra Deep Field_ UV luminosity function at $z=8$ is
$M_{\text{AB,lim}}\approx-18$ which we read off from Fig. 4 as giving a
limiting halo mass of $M_{h}\approx 10^{10}M_{\odot}$. Demanding
$M_{0}>M_{h}(M_{\text{AB,lim}})$ we find the bound:

$\displaystyle
m>$0.7\text{\times}{10}^{-22}$\left(\frac{\Omega_{\text{UBDM}}h^{2}}{0.12}\right)^{1/6}\left(\frac{h}{0.676}\right)^{2/3}\left(\frac{1+z_{\text{eq}}}{3390}\right)^{1/2}$\mathrm{eV}$$
$\displaystyle\hfill\text{(high redshift galaxies)}\,.$ (44)

The estimate Eq. (44) agrees favourably with complete analyses of similar data
Bozek:2014uqa; Schive:2015kza; Corasaniti:2016epp.

Another important bound to consider is from the Lyman-$\alpha$ forest flux
power spectrum. This observable traces the matter power spectrum, $P(k)$, 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.

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 $k$ as large as
$k_{\text{max}}\approx 50\,h\,$\mathrm{M}\mathrm{pc}^{-1}$$
Viel:2013apy.111111We convert from Lyman-$\alpha$ units for $k$ in s km-1 to
the more standard Mpc-1 by multiplication with $H_{0}(1+z)$. For reviews and
discussions of the Lyman-$\alpha$ forest as a probe of the matter power
spectrum see Refs. Gnedin:2001wg; Tegmark:2002cy; Hui:2016ltb;
Chabanier:2019eai. 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 $k_{\text{max}}$ gives the bound (cf.
Fig. 1):

$\displaystyle
m>$1.5\text{\times}{10}^{-21}\text{\,}\mathrm{eV}$\quad\text{(Lyman-}\alpha\text{
forest)}\,.$ (45)

which again agrees well with the result derived from more careful analysis
Armengaud:2017nkf; Irsic:2017yje.

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
$V(\phi)=f_{a}^{2}m^{2}[1-\cos(\phi/f_{a})]$. An initial displacement
$\theta=\phi/f_{a}=\pi-\delta\theta$ with small $\delta\theta$ leads to large
self-interactions at early times, and the field is near an unstable local
maximum of the potential. This _tachyonic instability_ 121212A potential is
said to have a tachyonic region if $V^{\prime\prime}(\phi)<0$, i.e. a local
maximum, and negative effective mass squared. in the evolution of $\delta$
leads to an _increase_ in the UBDM power spectrum relative to CDM on scales
close to the Jeans scale Zhang:2017dpp. A displacement $\delta\theta\approx
0.02$ is sufficient to evade the bound Eq. (45) and allow $m\approx
10^{-22}\text{ eV}$ to fit the Lyman-$\alpha$ power spectrum as well as CDM,
while a value $\delta\theta\approx 0.003$ leads to a _better_ fit than CDM
Leong:2018opi. The tuned values of $\delta\theta$ require smaller $f_{a}$ 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. Hui:2016ltb; 2019ApJ...871...28B; El-Zant:2019ios;
Lancaster:2019mde; Marsh:2018zyw; 2019MNRAS.485.2861C). We will give only the
simplest description useful for estimates.

UBDM in a gravitational potential well has a coherence length,
$L\sim\lambdabar_{\text{dB}}=1/mv$ ($\hbar/mv$ in physical units), and
coherence time $\tau\sim 1/mv^{2}$, where $v$ is the characteristic velocity.
The heuristic picture of a wave distribution with these properties is one of
quasi-particles of size $L$ and lifetime $\tau$. The quasiparticle mass is :

$\displaystyle M_{\rm qp}\sim\lambdabar_{\text{dB}}^{3}\bar{\rho}\,,$ (46)

where $\bar{\rho}$ is the average local density in a volume encompassing a
large number of quasi-particles (i.e., in the solar neighbourhood, 0.4 GeV cm
-3). Two body relaxation between quasi-particles leads to the relaxation time
(see Problem 3) Hui:2016ltb:

$\displaystyle t_{\rm
relax}\sim\frac{10^{10}}{\log\Lambda}\left(\frac{m}{10^{-22}\text{
eV}}\right)^{3}\left(\frac{v}{100\text{ km
s}^{-1}}\right)^{2}\left(\frac{R}{5\text{ kpc}}\right)^{4}\text{ yr}\,,$ (47)

where the Coulomb logarithm in the quasi-particle picture is
$\log\Lambda=\log(R/\lambdabar_{\text{dB}})$. On timescales longer than
$t_{\rm 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
$v^{2}=Gm\bar{n}R^{2}$.

Heating and cooling on the timescale $t_{\rm relax}$ can be observed if a
tracer population of stars with mass $m_{t}$ is present in the UBDM halo (when
the gravitational potential due to DM is dominant, stars are tracer
particles). For $m_{t}\ll M_{\rm qp}$, heating dominates, while for $m_{t}\gg
M_{\rm qp}$, cooling dominates. Let’s estimate $M_{\rm qp}$ for some systems
of interest. In the solar neighbourhood, $\bar{\rho}\approx 0.4\text{ GeV
cm}^{-3}=10^{7}\,M_{\odot}\text{ kpc}^{-3}$ and $v\approx 100\text{ km
s}^{-1}\Rightarrow\lambdabar=0.2(10^{-22}\text{ eV}/m)\text{ kpc}$, which
gives $M_{\rm qp}\approx 7\times 10^{4}(10^{-22}\text{ eV}/m)^{3}M_{\odot}$.
In the solar neighbourhood tracers are stars with $m_{t}\sim 1M_{\odot}$, and
the transition from heating to cooling occurs for UBDM mass $m\approx 4\times
10^{-21}\text{ eV}$, with lighter masses giving rise to heating. The Milky Way
in fact possesses a “thick disk” of old stars BinneyTremaine2008, and this has
been argued to provide evidence that in fact DM is composed of UBDM in this
so-called fuzzy DM regime Hui:2016ltb; 2019MNRAS.485.2861C (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}\text{ years}$, and applying Eq. (47) we find:

$\displaystyle m\gtrsim 10^{-22}\text{ eV}\,\quad\text{(Milky Way disk
heating)}\,,$ (48)

which agrees with more accurate modelling 2019MNRAS.485.2861C.

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 Li:2016utv; 2016ApJ...824L..14C indicate that
the DM density is $\bar{\rho}=0.15M_{\odot}\text{ pc}^{-3}$, and the velocity
dispersion is $\sigma_{v}=6.9^{+1.2}_{-0.9}\text{ km s}^{-1}$. For UBDM this
gives $M_{\rm qp}=3(10^{-19}\text{ eV}/m)^{3}M_{\odot}$, implying that heating
dominates for masses $m\lesssim 10^{-19}\text{ eV}$. The star cluster has a
half-light radius of $r_{h}=13\text{ pc}$, an estimated age $t\sim
10^{10}\text{ years}$, and is close to the centre of Eridanus II. Using Eq.
(47), replacing $R$ with the half-light radius (since the star cluster is
approximately centrally located), $v$ with $\sigma_{v}$, taking
$\log\Lambda\sim\mathcal{O}(1)$, and demanding that the star cluster is stable
on the time scale of its age, gives the bound:

$\displaystyle m\gtrsim 10^{-19}\text{ eV}\quad\text{(Eridanus~{}II)}\,,$ (49)

which again agrees very favourably with a more rigorous treatment
Marsh:2018zyw.

Based on the present analysis, the bound from Eridanus II does _not_ ,
however, apply for $m\lesssim 10^{-21}\text{ 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 $m\lesssim 10^{-21}\text{ 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 $m\approx 10^{-22}\text{ 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 $m\approx 10^{-22}\text{ eV}$ safe from this bound
Schive:2019rrw.

###### Question 3

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 $m$ moving with velocity
$v$ in a host of mass $M$ with radius $R$ can be written as
BinneyTremaine2008:

$\displaystyle t_{\text{relax}}=0.1\,\frac{R}{v}\,\frac{M}{m\log\Lambda}\,,$
(50)

Use this to derive the relaxation timescale, $t_{\rm relax}$, in Eq. (47).

Solution on page 8.

###### Example 3

UBDM Hints: Fuzzy Dark Matter We have seen a large variety of constraints on
UBDM with mass $m\lesssim${10}^{-22}\text{\,}\mathrm{eV}$$ from cosmic large
scale structure. We have also seen how heating in Eridanus II excludes the
range $${10}^{-21}\text{\,}\mathrm{eV}$\lesssim
m\lesssim${10}^{-19}\text{\,}\mathrm{eV}$$, and we will see shortly that Black
Hole superradiance excludes $${10}^{-19}\text{\,}\mathrm{eV}$\lesssim
m\lesssim${10}^{-16}\text{\,}\mathrm{eV}$$. There is only one strong bound in
the range just above ${10}^{-22}\text{\,}\mathrm{eV}$ coming from the
Lyman-$\alpha$ 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. Rogers:2020ltq). Another possible window is afforded by the
Eridanus II bounds around ${10}^{-21}\text{\,}\mathrm{eV}$, where the
statistical modelling is uncertain, and Eridanus II can survive sandwiched
between orbital resonances. If either of these bounds (Ly-$\alpha$ or Eridanus
II) can be relaxed, then there are some hints that DM may in fact be UBDM with
masses between about ${10}^{-22}\text{\,}\mathrm{eV}$ and
${10}^{-21}\text{\,}\mathrm{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 2019MNRAS.485.2861C.

  * •

Suppressed high-$z$ galaxy formation: The redshift of reionization is known to
be around $z_{\rm reion}\approx 8$. This relatively low value naturally
explained by FDM, which suppresses formation of galaxies at $z\gtrsim 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
Schive:2014dra; Marsh:2015wka.

  * •

Relic density: The relic density is naturally explained by an FDM ALP with
$f_{a}$ 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.

###  2.4 Black hole superradiance

In the following we adopt different units: so-called geometric units where
$G_{N}=c=1$.

Spinning black holes (BHs) are described by the Kerr metric, which has two
parameters: mass, $M$, and dimensionless spin $a_{J}=J/M\in[0,1]$. In “Boyer-
Linquist” coordinates the line element is:131313An accessible introduction to
general relativity can be found in Ref. 2004sgig.book.....C.

$\displaystyle{\rm d}s^{2}_{\rm Kerr}$ $\displaystyle=$
$\displaystyle-\left(1-\frac{2Mr}{\Sigma}\right){\rm
d}t^{2}-\frac{4Ma_{J}r\sin^{2}\theta}{\Sigma}{\rm d}t{\rm
d}\phi+\frac{\Sigma}{\Delta}{\rm d}r^{2}+\Sigma{\rm d}\theta^{2}$ (51)
$\displaystyle+$
$\displaystyle\frac{(r^{2}+a_{J}^{2})^{2}-a_{J}^{2}\Delta\sin^{2}\theta}{\Sigma}\sin^{2}\theta{\rm
d}\phi^{2}\,,$ $\displaystyle\Sigma$ $\displaystyle\equiv$ $\displaystyle
r^{2}+a_{J}^{2}\cos^{2}\theta\,,$ (52) $\displaystyle\Delta$
$\displaystyle\equiv$ $\displaystyle r^{2}+a_{J}^{2}-2Mr\,,$ (53)
$\displaystyle r_{\pm}$ $\displaystyle\equiv$ $\displaystyle
M\pm\sqrt{M^{2}-a_{J}^{2}}\,.$ (54) $\displaystyle r_{\text{ergo}}$
$\displaystyle\equiv$ $\displaystyle
M+\sqrt{M^{2}-a_{J}^{2}\cos^{2}\theta}\,,$ (55)

where we use spherical polar coordinates. The zero solutions of Eq. (53)
define the two horizons $r_{\pm}$: an inner Cauchy (causal) horizon at
$r_{-}$, and the outer physical event horizon at $r_{+}$. The “ergoregion” is
defined as radii smaller than $r_{\text{ergo}}$, where $g_{00}=0$ (the co-
efficient of $dt^{2}$ in the line element). If an object enters the ergoregion
between $r_{+}<r<r_{\text{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 $a_{J}=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. 2011PhRvD..83d4026A; 2015LNP...906.....B.

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 ($\Box$) should be evaluated with the metric Eq.
(51). Let us write the field as

$\displaystyle\phi=\sum_{\ell,\alpha}e^{-i\omega
t+i\mu\varphi}S_{\ell\mu}(\theta)\psi_{\ell\mu}(r)+\text{h.c.}\,,$ (56)

where $S_{\ell\mu}(\theta)$ 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 $\mu$ and the azimuthal angle $\varphi$. The Klein-Gordon equation can
then be reduced to a time-independent Schrödinger equation for the radial
eigenfunctions $\psi_{\ell\mu}$, with eigenvalue $\omega$. The BH provides a
background potential $V(r,\omega)$, which possesses a barrier separating the
bound states from the horizon, and a potential well with size of order the
boson Compton wavelength, $1/m$. The system resembles a hydrogen atom with
effective fine structure constant $\alpha_{\text{eff}}\equiv G_{N}Mm$, where
we temporarily re-instated $G_{N}$.

The existence of superradiant solutions is determined by the imaginary part of
the eigenvalue $\omega$, which leads to growth of the occupation number of the
mode $\psi_{\ell\mu}$. The superradiant rate is $\Gamma_{\rm
SR}\propto\alpha_{\text{eff}}^{4\ell+4}m$, and numerically it is found to be
maximised around $\alpha_{\text{eff}}\sim 1$. This gives an approximate
criterion for BHSR:

$\displaystyle m\sim\frac{8\pi M_{pl}^{2}}{M}=1.33\times 10^{-10}\text{
eV}\left(\frac{1M_{\odot}}{M}\right)\,.$ (57)

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 $a_{J}$ 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
$M\approx 1.4M_{\odot}$, 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,
$a_{J}\gtrsim 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 $M\approx 30M_{\odot}$
Abbott:2016blz. This excludes UBDM for:

$\displaystyle$4\text{\times}{10}^{-12}\text{\,}\mathrm{eV}$<m<$8\text{\times}{10}^{-11}\text{\,}\mathrm{eV}$\quad\text{(stellar
BHs)}\,.$ (58)

Next, consider SMBHs. The lightest currently known SMBH is in NGC4051, with
mass $M\approx 1.9\times 10^{6}M_{\odot}$, while the _Event Horizon Telescope_
has imaged the BH at the centre of M87 and determined the mass $M\approx
6.5\times 10^{9}M_{\odot}$. Again, assuming a continuous spectrum in between
we can exclude the range of UBDM masses:

$\displaystyle$2\text{\times}{10}^{-20}\text{\,}\mathrm{eV}$<m<$7\text{\times}{10}^{-17}\text{\,}\mathrm{eV}$\quad\text{(supermassive
BHs)}\,.$ (59)

These estimates agree somewhat favourably with more accurate treatments of
BHSR modelling and BH population statistics (e.g. Ref. Stott:2018opm).

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 $\tau\sim 1/m$. For real scalar fields the gravitational pressure
oscillates with a frequency $2m$, 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 Dolan:2018dqv considers spin-one vectors
which have much smaller instability rates, and thus weaker constraints.
Reference Brito:2013wya 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
$\mathbb{U}(1)$ symmetry and conserved particle number, the complex vector
$A_{\mu}$ field does not source oscillations in the metric potentials with
frequency $m$. This greatly simplifies the numerical task, and has allowed
direct simulation of superradiance with these so-called Proca fields
East:2017ovw. 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 self-interactions, which shut off the instability and
prevent growth of the scalar cloud. The self-interaction term in the potential
is $V_{\text{int}}=\lambda\phi^{4}/4!$, for some coupling constant $\lambda$.
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 $\lambda$ above which the superradiance rate is
subdominant to the Bosenova rate and no spin extraction can occur. Numerical
simulations Yoshino:2012kn determine the maximum cloud occupation number
before Bosenova occurs Arvanitaki:2014wva:

$\displaystyle N_{\text{Bose}}\sim
150\frac{n^{4}}{\alpha_{\text{eff}}\lambda}=$5\text{\times}{10}^{78}$\frac{n^{4}}{\alpha_{\text{eff}}^{3}}\left(\frac{M}{10M_{\odot}}\right)^{2}\left(\frac{f_{a}}{M_{pl}}\right)^{2}\,,$
(60)

where $n$ is the energy level of the occupied cloud and $M_{pl}$ is the
reduced Planck mass. In the last equality we assumed that the scalar potential
is of the ALP form $V(\phi)\propto-\cos(\phi/f_{a})$, giving
$\lambda=m^{2}/f_{a}^{2}$. Using this formula for stellar mass BHs, Ref.
Arvanitaki:2014wva finds that BHSR is shut off for
$f_{a}\lesssim${10}^{13}\text{\,}\mathrm{GeV}$$; for SMBHs this turns out to
be $f_{a}\lesssim${10}^{16}\text{\,}\mathrm{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 $g_{a\gamma\gamma}$,
which leads to stimulated decay of the cloud Rosa:2017ury; Ikeda:2018nhb. 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
Day:2019bbh – even though the astrophysical uncertainties are far greater.

###### Question 4

Problem 4: Estimating superradiance properties of UBDM A simple way to
estimate the relevance of BHSR is to inspect terms in the action,

$\displaystyle S=\int{\rm
d}^{4}x\,\sqrt{-g}\left[\frac{M_{pl}^{2}}{2}R-\frac{1}{2}(\partial\phi)^{2}-\frac{1}{2}m^{2}\phi^{2}+\frac{\lambda}{4!}\phi^{4}\right]\,,$
(61)

where $R$ is the Ricci scalar of the Kerr background metric $g$ and $m$,
$\lambda$, and $\phi$ are the UBDM mass, (dimensionless) self-coupling, and
field value, respectively. Note that it is useful to re-instat $M_{pl}$ (or
$G_{N}$) 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
$N_{\text{Bose}}$ and Eq. (60) when $\lambda=m^{2}/f_{a}^{2}$.

Solution on page 9.

###### Example 4

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 $(m,a_{J})$ 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}\text{ eV}\lesssim
m\lesssim 10^{-12}\text{ eV}$ Arvanitaki:2016qwi. 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 $f_{a}\sim M_{pl}$. 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.

###  2.5 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 $w_{\text{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 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_
Khmelnitsky:2013lxt. Current bounds from this technique Porayko:2018sfa are
not yet at the $\mathcal{O}(1)$ level for $\Omega_{\text{UBDM}}$, and so do
not appear.

Figure 5: Summary of gravitational constraints (shaded) on UBDM, and forecasts
(open) for upcoming surveys. Constraints assume a real scalar with potential
$V(\phi)=m^{2}\phi^{2}/2$, see text for clarification on generalizing the
bounds. CMB: cosmic microwave background Hlozek:2017zzf; Hlozek:2016lzm, PTA:
pulsar timing array Khmelnitsky:2013lxt, BHSR: black hole superradiance
Stott:2018opm, Ly-a: lyman alpha forest Kobayashi:2017jcf; Aghamousa:2016zmz.
SKA-IM: Square Kilometer Array intensity mapping Bauer:2020zsj. Adapted from
Ref. Grin:2019mub.

##  3 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.

###  3.1 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 1969PhRv..187.1767R 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 $\mathbb{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” 1991PhRvL..66.1659S. 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.141414To continue the
bestiary just a little further, solutions are named for all scalar fields:
inflaton stars, moduli stars, Higgs stars, etc.

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 GRChombo Clough:2015sqa.151515http://www.grchombo.org/ Like all stars,
axion stars are stabilised by a balance between attraction (gravity, and axion
quartic self-interactions) and repulsion (gradient pressure, and higher order
interactions).161616The axion potential is
$V(\phi)=m^{2}f_{a}^{2}[1-\cos\phi/f_{a})]$. Taylor expanding this we find
that the $\phi^{4}$ term is attractive, while higher order terms alternate in
sign. 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,
$M$, and axion decay constant, $f_{a}$, giving a “phase diagram” that can be
explored numerically Helfer:2016ljl.

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 $\phi_{0}$ also
increases. There are two possible instabilities, and which wins depends on
$f_{a}$. For large $f_{a}$, the self-interactions can be neglected. Now the
ordinary GR lore applies: collapse to a BH at large mass. At low $f_{a}$, the
axion has strong self-interactions, and these also drive collapse. Collapse
increases $\phi_{0}$ further until higher order repulsive interactions take
over and expel relativistic axions from the collapsing core in an “axion nova”
Levkov:2016rkk, which occurs at critical mass
$M_{\text{nova}}=10.4\,M_{pl}f_{a}/mg_{4}$, where $g_{4}$ is the coefficient
of quartic interactions equal to unity for a cosine potential. For small
$f_{a}$ 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 $f_{a}$ 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 $f_{a}$. As
$f_{a}\rightarrow\infty$, oscillatons and boson stars are found to be unstable
when $\phi_{0}\sim M_{pl}$ (this defines the “Kaup mass”, $M\sim
M_{pl}^{2}/m$), while self-interactions become important when $\phi_{0}\sim
f_{a}$, and so the boundary between the two unstable regions occurs for
$f_{a}\sim M_{pl}$. A third phase boundary exists between the nova and BH
regions, which simulations have found to be fractal in structure
Michel:2018nzt. 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
$(M,f_{a})=(2.4M_{pl}^{2}/m,0.3M_{pl})$, where $M$ is the “Arnowitt-Deser-
Misner” mass 1959PhRv..116.1322A.171717Due 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.

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 $\psi(r,t)=M_{pl}\chi(mr)e^{-i\gamma mt}$, where
$\chi$ is a dimensionless function giving the radial profile, and $\gamma$ is
the energy eigenvalue. An important property of the SPEs (see Sec. 2.2) is
their _scaling symmetry_ :

$\displaystyle(t,x,\psi,\Phi)\rightarrow(\lambda^{-2}t,\lambda^{-1}x,\lambda^{2}\psi,\lambda^{2}\Phi)\,,$
(62)

where $\lambda$ is the scale parameter (not to be confused with the quartic
interaction strength in Section 2.4). The boundary value problem normalised to
$\chi(0)=1$, $\lambda=1$, can be solved numerically and the results are fit by
eigenvalue $\gamma=-0.692\lambda^{2}$ and radial density profile:

$\displaystyle\frac{\rho_{\rm
sol}(r)}{m^{2}M_{pl}^{2}}=\chi^{2}(mr)=\frac{1}{[1+(0.230mr)^{2}]^{8}}\,.$
(63)

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, $M$. Thanks to the scaling symmetry, we
only need to find the solution once, and then scale it using $\lambda$ (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, $R$, of the
gravitational fluctuations compared to the de Broglie wavelength:

  * •

Direct collapse: $R\sim\lambda_{\text{dB}}$ (e.g. Ref. Schive:2014dra).

  * •

Kinetic condensation: $R\gg\lambda_{\text{dB}}$ (e.g. Ref. Levkov:2018kau).

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., $M\sim M_{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, $M_{\star}$, and the halo mass, $M$, is:

$M_{\star}\propto\left(\frac{M}{M_{0}}\right)^{1/3}M_{0}\,,$ (64)

where the constant of proportionality can be found in Ref. Schive:2014hza, and
depends on the definition of $M_{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 Du:2016aik. Slow growth of solitons by accretion
leads to significant scatter in the relation.181818Very recently some authors
have even found a different best-fit exponent Mina:2020eik; Nori:2020jzx, and
numerical convergence may also play a part. The issue is not yet resolved at
the time of writing.

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
$t<t_{\rm 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 Veltmaat:2018dfz.

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 $\tau_{\rm 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 $M\propto(t/\tau_{\rm
gr})^{p}$. The index $p$ 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. (64) breaks down and
halos have no central soliton, but instead grow many small solitons in the
kinetic regime.

Axion stars have a host of possible phenomenological consequences:

  * •

Galactic cores: Solitons composed of Fuzzy DM with $m\sim 10^{-22}\text{ eV}$
may help explain flat central densities in Milky Way dwarf satellites (tracer
stars reside within the soliton) Schive:2014dra; Marsh:2015wka, 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
JacksonKimball:2017qgk; Derevianko:2013oaa.

  * •

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
Giudice:2016zpa and multi-messenger astronomy Dietrich:2018jov.

###  3.2 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 conditions
epoch, inflation or otherwise). The Peccei-Quinn (PQ) phase transition (see
Chapter 2) occurs when the temperature of the Universe drops below
approximately $f_{a}$. Recall that we write the complex PQ field as
$\varphi=Re^{i\theta}$, and spontaneous symmetry breaking occurs when the
field $R$ takes on a vacuum expectation value. The following scenario applies
specifically to ALPs where the field $R$ is heavy and unstable (such that it
decays at late times), while the field $\theta$ is initially massless, but
acquires a mass hierarchically smaller than the mass of $R$ and at some time
much later than the time of PQ symmetry breaking.

When PQ symmetry breaking occurs, $R$ takes on a non-zero vacuum expectation
value, and thus $\theta$ must also be specified. Since the axion field is
massless at symmetry breaking, the only terms in the Lagrangian are
proportional to $\partial\theta$, meaning there can be no preferred value for
$\theta$. 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,
$\theta$, with the $\theta=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 $\theta$ makes a complete wrapping. At the wrapped point, the axion
field $\theta$ is undefined (imagine shrinking the circle to a point: at the
point the circle must have zero size, and $\theta$ takes every value at once).
The only way that this can be possible is if $R=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 $\theta$ persist,
then at the centre of these windings $R$ must remain at the origin, and thus
symmetry breaking cannot happen. When the potential is
$V(\varphi)=\lambda(|\varphi|^{2}-f_{a}^{2}/2)^{2}$ this implies that the
potential at the origin is $V(0)=\lambda f_{a}^{4}/4$, and this is the value
of the potential at the centre of a point around which $\theta$ wraps.

In fact, in three dimensions, $\theta$ 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 $\mathbb{U}(1)$).191919This is
generic for complex fields in three spatial dimensions. It is a topological
property. Complex fields have symmetry group $\mathbb{U}(1)$ of rotations in
the complex plane, i.e. loops. The mapping of $\mathbb{U}(1)$ onto
$\mathbb{R}^{3}$ (Euclidean 3-space) is expressed by the first homotopy group
$\pi_{1}(\mathbb{R}^{3})$. This group is not the empty set, i.e. it is non-
trivial, which can be seen by noting that $\mathbb{R}^{3}$ is the universal
cover of $T^{3}$, the 3-torus, and $\pi_{1}(T^{n})=\mathbb{Z}^{n}$. This last
can be seen since one cannot shrink circles on tori to points continuously,
and there are $n$ distinct circles wrapping $T^{n}$. This leads to the
existence of PQ strings: continuous one dimensional structures around which
$\theta$ makes complete windings and where the radial field is pinned at
$R=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
1976JPhA....9.1387K; 1985Natur.317..505Z.

What happens to the axion field? The equation of motion for the axion field in
Fourier space is:

$\displaystyle\ddot{\tilde{\theta}}_{k}+3H\dot{\tilde{\theta}}_{k}+(k^{2}/a^{2}+m(T)^{2})\tilde{\theta}_{k}=0\,,$
(65)

where we are careful to distinguish between the field $\theta$ and the
modefunction $\tilde{\theta}_{k}$: $\tilde{\theta}_{k}=0$ does not imply
$\theta=0$ as a preferred value, only that the mode $k$ is absent from its
spectrum and thus gradients of $\theta$ on the spatial scale $1/k$ are small.
At early times, the QCD axion mass, $m(T)$ is vanishingly small comapred to
$H$ 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 $k^{2}\gg a^{2}H^{2}$. These modes will undergo damped
oscillation, and decay in amplitude. Modes larger than the horizon, $k^{2}\ll
a^{2}H^{2}$ remain pinned to their initial value by the friction term (co-
efficient of $\dot{\tilde{\theta}}_{k}$) in Eq. (65) given by $3H$ (Hubble
friction). Thus high frequency modes decay and low frequency modes remain
static, smoothing the field on scales of order the horizon size.202020Note
that the mode function $\tilde{\theta}_{k}$ decaying to zero does not imply a
preference for the axion field $\theta$ to move to zero: in the massless limit
the shift symmetry prevents any such preference. Around any string, the axion
field is wound $\theta\in(-\pi,\pi]$, and so we have $\mathcal{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 $\mathcal{O}(1)$
strings per horizon volume.

Figure 6: Formation of axion strings from spontaneous symmetry breaking. The
complex plane of the PQ field $\varphi$ is illustrated in the upper right
corner with the radius marked where $R$ takes the vacuum expectation value,
and the circle along which the axion field is defined. When spontaneous
symmetry breaking occurs, this circle is mapped onto the co-ordinates in the
real plane, $(x,y)$. The complex phase of the field (i.e. the value of the
axion field) is shown by the orientation of the arrow. Complete wrappings of
the field lead to defects where the phase is not defined. In this two-
dimensional case, the defects are points, in thee dimensions they are lines,
i.e. strings. In the centre of the string, the PQ field has
$\langle\varphi\rangle=0$ and the potential takes the value indicated by the
colour bar. Damping of sub-horizon modes $k>aH$ in Eq. (65) smooths
fluctuations on length scales $L\sim H^{-1}$.

Strings decay when the axion mass becomes relevant to the mode evolution.
Recall first that during cosmological expansion temperature, $T$, always
decreases as time, $t$, increases. Second, recall that the axion mass and the
Hubble rate $H$ are both decreasing functions of $T$. The mass term (co-
efficient of $\tilde{\theta}_{k}$) in Eq. (65) is comparable to the friction
term when $m(T)\approx H(T)$. The axion mass term defines a preferred value
for the field, $\theta=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 $k<aH$ (not just the short wavelength modes inside the
horizon) will begin to oscillate when $m(T)\approx H(T)$, defining the special
temperature $T_{\text{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 $R=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 $\mathcal{O}(1)$ fluctuations, while on larger
scales it is uncorrelated. The power spectrum, $P(k)$, is flat (white noise)
for $k\ll a(T_{\text{osc}})H(T_{\text{osc}})=k_{J}(T_{\text{osc}})$ and cut
off by the Jeans scale for $k\gg a(T_{\text{osc}})H(T_{\text{osc}})$. The
normalisation of the power spectrum is fixed by the variance, which should
match the variance of the uniform distribution for $\theta$,
$\langle\theta^{2}\rangle=\pi^{2}/3$. Just prior to $a(T_{\text{osc}})$ the
axion equation of state is $w\approx-1$, and so the fluctuations in $\theta$
do not source any curvature perturbations in the metric:212121Intuitively this
can be understood because if $w=-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. this is what is meant by the term
_isocurvature_. It is this particular power spectrum (white noise
isocurvature, with $\mathcal{O}(1)$ variance, truncated at the horizon size at
$T_{\text{osc}}$), which gives rise to the structures known as _axion
miniclusters_ as a remnant of string decay 1988PhLB..205..228H. Figure 7 shows
a snapshot from numerical simulation of the axion field after string decay
Vaquero:2018tib, and miniclusters are located using a threshold based on
spherical collapse under gravity Ellis:2020gtq.

Figure 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. Vaquero:2018tib, 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 Ellis:2020gtq.

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.222222Contrast 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 $m\lesssim
1\text{ TeV}$, structure formation in the adiabatic mode begins at $z\approx
500$ with the formation of Earth mass, $10^{-6}M_{\odot}$, halos Green:2003un.
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 $T_{\text{osc}}$ in the usual way, setting
$H(T_{\text{osc}})=A_{1}m_{a}(T_{\text{osc}})$.232323$A_{1}$ is simply a
constant of proportionality to account for ambiguity defining
$T_{\text{osc}}$, and its later use in analytical formulae for the evolution
of the energy density. Our earlier choices, e.g. Eq. (4), set $A_{1}=1$,
physically assuming structure formation begs when the de Broglie wavelength is
equal to the Hubble length, $H^{-1}$. Many authors choose $A_{1}=3$ when
estimating the relic density. This ambiguity in the use of $T_{\text{osc}}$
leads to significant uncertainty in analytical minicluster mass estimates,
which can only be resolved by fitting results of numerical simulations. Now we
need to calculate the horizon volume. Take the volume of rotation of the
spherical wave with comoving wavenumber $k_{\text{osc}}=a(T_{\rm
osc})H(T_{\text{osc}})$ over one-half wavelength:

$\displaystyle V(k)=\int{\rm
d}^{3}xj_{0}(kr)=\frac{4\pi}{(ak)^{3}}\int_{0}^{\pi}{\rm
d}y\,y\sin(y)=\frac{4\pi^{2}}{(ak)^{3}}\equiv\frac{V_{H}}{(ak)^{3}}\,.$ (66)

This defines the Hubble volume $V_{H}$. Alternative definitions are the cubic
volume, $V(k)=(ak)^{-3}$, and the spherical volume of of one half wavelength
$V(k)=\frac{4}{3}\pi(\pi/ak)^{3}$. The expected minicluster mass is simply
$M_{\text{MC}}=V(a_{\text{osc}}k_{\text{osc}})\Omega_{a}\rho_{\text{crit}}$.
We compute $M_{\text{MC}}(T_{\text{osc}})$ using the Friedmann equation to fix
$k_{\text{osc}}$, $3H^{2}M_{pl}^{2}=(\pi^{2}/30)g_{\star}T^{4}$, and
conservation of entropy to write $a(T)\propto T^{-1}g_{\star S}$ (normalised
using the CMB measurement of $z_{\text{eq}}$). The result is well fit by:

$\displaystyle
M_{\text{MC}}(T_{\text{osc}})=$9.2\text{\times}{10}^{-13}$\,V_{H}\,\left(\frac{T_{\text{osc}}}{2\text{
GeV}}\right)^{-3}\mathcal{S}\left(\log_{10}\left(\frac{T_{\text{osc}}}{$\mathrm{GeV}$}\right)\right)\,M_{\odot}\,,$
(67) $\displaystyle\mathcal{S}(x)\equiv 0.5[1+\tanh 4(x-8.2)]+1.3[1+\tanh
4(8.2-x)]\,,$ (68)

where $\mathcal{S}(x)$ is an activation function that accounts for the
behaviour of $g_{\star}$ and $g_{\star,S}$ in the Standard Model (dominantly
the quark-hadron phase transition) and has only been roughly fit here using
the results for $g_{\star}$ from Ref. Wantz:2009it. Note that this expression
is valid for any ALP with minicluster-like initial conditions.

From the QCD axion $m_{a,\text{QCD}}(T)$ dependence based on lattice QCD in
Ref. Borsanyi:2016ksw, $T_{\text{osc}}(m_{a,\text{QCD}})$ is well fit
by:242424For the present purposes a simple power law $m\propto T^{-4}$,
matched to the zero-temperature result $m_{a,\text{QCD}}\equiv
m_{a,\text{QCD}}(0)=$5.72\text{\,}\mathrm{\SIUnitSymbolMicro
eV}$(${10}^{12}\text{\,}\mathrm{GeV}$/f_{a})$ at $T_{\rm
QCD}=$140\text{\,}\mathrm{MeV}$$ is accurate enough.

$\displaystyle
T_{\text{osc}}(m_{a,\text{QCD}})=2\,\left(\frac{m_{a,\text{QCD}}}{100A_{1}\,$\mathrm{\SIUnitSymbolMicro
eV}$}\right)^{0.165}\text{ GeV}\,,$ (69)

over the range of interest (broadly $${10}^{-5}\text{\,}\mathrm{eV}$\leq
m_{a,\text{QCD}}\leq${10}^{-3}\text{\,}\mathrm{eV}$$) for the relic density in
this scenario (see also Fig. 1). For smaller axion masses there there is an
$\mathcal{O}(1)$ change in the constant at the front of Eq. (69), while the
power law remains approximately the same.

Now we would like to know the minicluster density profile. Kolb and Tkachev
Kolb:1994fi wrote down the equation of motion for spherical collapse of an
isolated top-hat density profile, with initial overdensity $\delta$, 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. sigl) 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:252525For the more standard case during matter domination, which
applies to ordinary DM halos, see, e.g., Ref. peacock. In the standard case,
the equations can be solved analytically, leading to the well-known result
that the virial overdensity $\langle\rho_{f}\rangle/\bar{\rho}\approx 200$,
independent of $\delta$.

$\displaystyle\langle\rho_{f}\rangle=140\bar{\rho}_{\text{eq}}\delta^{3}(1+\delta).$
(70)

Assuming the minicluster has a constant density, the radius is then calculated
to be

$\displaystyle
R_{\text{MC}}=\left(\frac{3M_{\text{MC}}}{4\langle\rho_{f}\rangle}\right)^{1/3}\,.$
(71)

Using $\bar{\rho}_{\text{eq}}=2\times 10^{14}M_{\odot}\text{kpc}^{-3}$

$\displaystyle R_{\text{MC}}=3\times
10^{-10}\,\text{kpc}\,\frac{1}{\delta(1+\delta)^{1/3}}\,\left(\frac{T_{\mathrm{osc}}}{2\text{GeV}}\right)^{-1}\,.$
(72)

To make more use of this result, we need to know the minicluster radial
profile, $\rho(r)$. 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 Zurek:2006sy;
Gosenca:2017ybi; Eggemeier:2019khm. It is then natural to associate the radius
$R_{\mathrm{MC}}$ with the NFW scale radius. If the initial distribution of
$\delta$ 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 Enander:2017ogx;
Fairbairn:2017sil; Eggemeier:2019khm; Ellis:2020gtq. The mass function takes
on a power law spreading over many orders of magnitude around $M_{\text{MC}}$.
Secondly, the distribution of $\delta$ in initial conditions is not uniquely
defined. Numerical thresholding using the spherical collapse results allows
some progress to be made Ellis:2020gtq 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 Eggemeier:2019jsu. 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 Kolb:1995bu; Katz:2018zrn; Fairbairn:2017dmf
(see Problem 5 below).

  * •

Radio signals from minicluster-neutron star collisions Tkachev:2014dpa;
Iwazaki:2014wta; Bai:2017feq.

  * •

The CMB isocurvature power spectrum and large scale structure Hardy:2016mns;
Feix:2020txt.

###### Question 5

Problem 5: Microlensing constraints on UBDM Show that the scaling symmetry,
Eq. (62), is a symmetry of the SPEs. Use this relationship and the profile,
Eq. (63), to write down the mass-radius relation in units of solar masses
$(M_{\odot})$ and kiloparsecs. The _Subaru Hyper Suprime Cam_ (HSC)
microlensing survey of M31 probes PBHs in the range $10^{-11}M_{\odot}$ to
$10^{-6}M_{\odot}$ Niikura:2017zjd. The Einstein radius for gravitational
microlensing, is $R_{E}=2[GM_{\star}x(1-x)d_{s}]^{1/2}$, where $d$ is the
distance from the observer to the lens, $d_{s}$ is the distance from the
observer to the source, and $x=d/d_{s}$. Compare the axion star radius to
$R_{E}$ with $x=1/2$ and $d_{s}=770\text{ 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 $\delta$, Eq. (72) and the minicluster mass relation Eq. (68).
What range of the $(T_{\text{osc}},\delta)$ parameter space can be probed by
microlensing?

Solution on page 10.

##  4 Indirect detection of UBDM

###  4.1 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 mesa mesa, which can be modified to include axion-induced cooling
Friedland:2012hj. For more details see Refs Raffelt:1990yz;
1990eaun.book.....K; 2008LNP...741...51R.

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:

$\displaystyle\Gamma_{\gamma\rightarrow
a}=\frac{g_{a\gamma\gamma}^{2}T\kappa_{s}^{2}}{32\pi}\left[\left(1+\frac{\kappa_{s}^{2}}{4E^{2}}\right)\ln\left(1+\frac{4E^{2}}{\kappa_{s}^{2}}\right)-1\right]\,,$
(73)

where $E$ is the photon energy, $T$ is the temperature, and $\kappa_{s}$ is
the screening length. In the Debye-Hückel approximation we have

$\displaystyle\kappa_{s}^{2}=\frac{4\pi\alpha}{T}\left(n_{e}+\sum_{j}Z_{j}^{2}n_{j}\right)\,,$
(74)

where $n_{e}$ is the free electron density, and $n_{j}$ is the density of the
$j$-th nuclear ion of charge $Z_{j}$. In a neutral medium with $n_{e}=n_{j}=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 $E\geq m$ 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
$L_{a}=1.85\times 10^{-3}(g_{a\gamma\gamma}/10^{-10}\text{
GeV}^{-1})^{2}L_{\odot}$. 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. $L_{\gamma,{\rm Sun}}=1L_{\odot}$ by definition. Thus:

$\displaystyle g_{a\gamma\gamma}<2.3\times 10^{-9}\text{
GeV}^{-1}\quad\text{(luminosity of the Sun)}\,.$ (75)

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, $R$, determined by placing stars on a
colour-magnitude diagram. In the Standard Model, this ratio is a function of
the primordial helium abundance, $Y_{\rm 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 $R$. The value of $Y_{\rm He}$ can determined
observationally by measurement of extragalactic H II regions which gives
$Y_{\rm He}=0.2449\pm 0.0040$ Aver:2015iza, which is consistent with the
predictions of standard BBN and the CMB measurement of the baryon abundance
Aghanim:2018eyx. Reference Ayala:2014pea reports a measured average value of
$R=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, $g_{a\gamma\gamma}$, and the axion-electron
coupling, $g_{aee}$, on $R$. At present there are two different models in the
literature for the functional dependence, and each is presented in Ref.
Giannotti:2015kwo. The specific forms are not enlightening, so we simply quote
the bounds (derived in Ref. Hoof:2018ieb). In both cases the additional
cooling channel lowers $R$ 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 $Y_{\rm He}=0.25$, the individual bounds are:

$\displaystyle g_{a\gamma\gamma}$ $\displaystyle<4.95\;(9.56)\times
10^{-11}\,\text{GeV}^{-1}\quad$ $\displaystyle\text{(95\% CL, HB/RGB
stars)}\,,$ (76) $\displaystyle g_{aee}$ $\displaystyle<2.95\;(3.53)\times
10^{-13}\quad$ $\displaystyle\text{(95\% CL, HB/RGB stars)}\,,$ (77)

where the number in brackets refers to the bound using the alternative model
for $R$, which we see introduces an $\mathcal{O}(1)$ shift in the bound on
$g_{a\gamma\gamma}$. For the constraints in the combined parameter space, see
Refs. Giannotti:2015kwo; Hoof:2018ieb.

Supernova SN1987A provides an important bound on the axion nuclear couplings,
$g_{aNN}$ and $g_{d}$. 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:

$\displaystyle N+N$ $\displaystyle\rightarrow N+N+a\quad$
$\displaystyle(\text{$g_{aNN}$ coupling})\,,$ (78) $\displaystyle N+\gamma$
$\displaystyle\rightarrow N+a\quad$ $\displaystyle(\text{$g_{d}$
coupling})\,.$ (79)

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 $\varepsilon_{a}\lesssim 10^{19}\,\text{erg
g}^{-1}\text{s}^{-1}=$7.2\text{\times}{10}^{-18}\text{\,}\mathrm{eV}$$.

For the first process, the cooling rate per unit mass is 2008LNP...741...51R

$\displaystyle\varepsilon_{a}=\frac{1}{\rho}\left(\frac{C_{N}}{2f_{a}}\right)^{2}\frac{n_{N}}{4\pi^{2}}\int_{0}^{\infty}{\rm
d}\omega\,\,\omega^{4}S_{\sigma}=\left(\frac{C_{N}}{2f_{a}}\right)^{2}\frac{T^{4}}{\pi^{2}m_{N}}F\,,$
(80)

where $\rho$ is the mass density of the supernova, $n_{N}$ is the nucleon
number density, $\omega$ is the axion angular frequency, and $S_{\sigma}$ is
the spin density structure function, which accounts for the fact that axions
couple only to the nuclear spins. The coupling $C_{N}$ is the nucleon coupling
weighted as $C_{N}^{2}=Y_{n}C_{n}^{2}+Y_{p}C_{p}^{2}$, where $Y_{n}$ and
$Y_{p}$ are the neutron and proton abundances, estimated as $Y_{p}=0.3$ and
$Y_{n}=1-Y_{p}=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 $F$ from the integral of $S_{\sigma}$, which is
estimated to be of order unity, and allows a simple estimate of the bound
given the supernova internal temperature $T\approx 30\text{ MeV}$. A recent
analysis found the more accurate bound Chang:2018rso

$\displaystyle\frac{C_{N}}{2f_{a}}<$1.3\text{\times}{10}^{-9}\text{\,}\mathrm{G}\mathrm{eV}^{-1}$\quad\text{(SN1987A
neutrino burst)}\,,$ (81)

a factor of approximately four weaker than the estimate with $F=1$. The same
analysis found an $\mathcal{O}(1)$ effect from the modelling of supernova
temperature and density profiles. The bound from SN1987A on $C_{N}/f_{a}$ 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 $C_{KSVZ}^{2}=0.066\Rightarrow C_{KSVZ}=0.257$,
leading to:

$\displaystyle m_{a,\text{QCD}}\lesssim
0.06\left(\frac{0.257}{C_{N}}\right)\text{\,eV}\quad\text{(SN1987A neutrino
burst)}\,.$ (82)

For the second process $\varepsilon_{a}$ is approximated by:

$\displaystyle\varepsilon_{a}=\frac{\Gamma}{V\rho}\approx\frac{\langle
E_{\gamma}\rangle n_{N}n_{\gamma}\langle\sigma v\rangle}{\rho}\,,$ (83)

where $n_{i}$ are the reactant number densities, $\langle E_{\gamma}\rangle$
is the average photon energy, $\rho$ is the supernova mass density, and
$\langle\sigma v\rangle$ 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, $T$. To estimate the bound on
$g_{d}$, Ref. Graham:2013gfa approximates the cross section as $\langle\sigma
v\rangle=g_{d}^{2}T^{2}$, leading to:

$\displaystyle g_{d}<4\times 10^{-9}\text{ GeV}^{-2}\,.$ (84)

###### Example 5

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 ($\sim
0.6\,M_{\odot}$), while their radii are of the order of the Earth’s radius
10.1086/156841. 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. Corsico:2019nmr), with
potentially multiple pulsation periods associated with different co-existing
sub-modes. For cooling WDs, the periods of their pulsations, $\Pi$, tend to
increase over time with a rate $\dot{\Pi}=\mathrm{d}\Pi/\mathrm{d}t$ that can
approximately be calculated via

$\frac{\dot{\Pi}}{\Pi}\approx-\frac{1}{2}\,\frac{\dot{T}}{T}+\frac{\dot{R}}{R}\approx-\frac{1}{2}\,\frac{\dot{T}}{T}\,,$
(85)

where $T$ and $R$ are the internal temperature and radius of the WD,
respectively 10.1038/303781a0. The change in radius can usually be neglected
for the observed, low-luminosity dwarfs 10.1086/163398. Axions and ALPs induce
an energy loss that is proportional to $g_{aee}^{2}$ since axion-electron
interactions dominate in the high-density, electron-degenerate interior of the
WDs Nakagawa:1987pga; Nakagawa:1988rhp. The resulting decrease in temperature,
and therefore the additional contribution to the period increase in Eq. (85),
is hence also proportional to $g_{aee}^{2}$.

Measurements of an anomalous period change can thus be used to estimate the
associated axion-electron coupling. From the 250 known variable WDs
10.1093/mnras/stz2571, this has so far only been done for G117-B15A
Isern:1992gia; Altherr:1993zd; BischoffKim:2007ve; Isern:2008nt;
Corsico:2012ki, R548 BischoffKim:2007ve; Corsico:2012sh, L19-2
Corsico:2016okh, and PG 1351+489 Battich:2016htm. 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, $\dot{\Pi}$, are less than about
$10^{-13}$ in magnitude.

Figure 8: Cooling hints in white dwarfs. Left: Comparison of the predicted and
measured period change as a function of axion-electron coupling squared for
the WD variable G117-B15A (data from Refs. 2012Kepler; Corsico:2012ki). Right:
Overview of the likelihood functions for different WD variables (reproduced
from Ref. Hoof:2018ieb).

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 $g_{aee}^{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\sigma$ confidence level Giannotti:2015kwo;
Giannotti:2017hny; Hoof:2018ieb.

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. Corsico:2012ki; Corsico:2012sh; Corsico:2016okh;
Battich:2016htm. 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.
Corsico:2019nmr. 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 Bertolami:2014wua; Isern:2018uce.

###  4.2 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
Chupp:1989kx. This places constraints on the axion mass and coupling at 95% CL
Payez:2014xsa:

$g_{a\gamma\gamma}<$5\text{\times}{10}^{-12}\text{\,}\mathrm{G}\mathrm{eV}^{-1}$\quad(\text{for
}m_{a}\lesssim${10}^{-9}\text{\,}\mathrm{eV}$)\,.$ (86)

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 $\mathcal{O}(1)$ dependence on the
precise model of the galactic magnetic field.

Axion-photon conversion also occurs in the intergalactic medium, and leads to