FERMILAB-PUB-24-0295-T

# Axion Stars: Mass Functions and Constraints

Jae Hyeok Chang  Patrick J. Fox  and Huangyu Xiao

###### Abstract

The QCD axion and axion-like particles, as leading dark matter candidates, can
also have interesting implications for dark matter substructures if the
Peccei-Quinn symmetry is broken after inflation. In such a scenario, axion
perturbations on small scales will lead to the formation of axion miniclusters
at matter-radiation equality, and subsequently the formation of axion stars.
Such compact objects open new windows for indirect searches for axions. We
compute the axion star mass function based on recent axion minicluster studies
and Bose star simulations. Applying this mass function, we find post-inflation
axion-like particles with masses $m_{a}<3.3\times 10^{-17}$ eV are constrained
by the lack of dynamical heating of stars in ultrafaint dwarfs. We also find
that current microlensing surveys are insensitive to QCD axion stars. While we
focus on the gravitational detectability of axion stars, our result can be
directly applied to other interesting signatures of axion stars, e.g. their
decay to photons, that require as input the abundance, mass, and density
distribution of axion stars.

##  1 Introduction

The QCD axion, currently a leading solution to the strong CP problem in
quantum chromodynamics (QCD) that addresses the experimental nearly zero value
of the neutron electric dipole moment, is also a viable dark matter candidate
[1, 2, 3, 4, 5, 6, 7] and has gained significant attention in recent years,
stimulating many experimental efforts on the direct detection of axion dark
matter [8, 9, 10, 11, 12]. In the post-inflationary scenario of axion, Peccei-
Quinn symmetry is broken after inflation and axions develop different field
values across different horizons, leading to a predictable relic abundance
that is independent of initial conditions [13, 14], which sets an important
target for QCD axion experiments. For those axion masses, $m_{a}\lesssim\rm
meV$, the occupation number of axion dark matter is so high that it can be
considered as wave-like dark matter which acts as a coherently-oscillating
classical field. On the other hand, the current spectrum of axion
possibilities is broad, which motivates us to consider axion-like particles
(ALPs) that do not solve the strong CP problem but have a broader range of
axion parameters [15].

The post-inflationary scenario of axions or axion-like particles involves
interesting cosmological consequences such as the formation of axion
miniclusters, also referred to as minihalos, and axion stars as early as
matter-radiation equality [16, 17, 18], which opens a new window for indirect
searches for axions that are complimentary to existing efforts on direct
detection. The search for signatures of axion minihalos is not as mature as
the search for the diffuse axion background. Specifically, they have not
imposed any constraints on the parameter space of QCD axions. This is
partially due to the computational challenge of simulating the evolution of
axion miniclusters but more because of the observational challenge of
detecting them gravitationally. Even though axion miniclusters form at very
early times, they are still too diffuse to be detectable with femto-, pico-,
and microlensing surveys [18, 19, 20]. There are a few future proposals that
might eventually detect QCD axion miniclusters such as pulsar timing arrays
[21, 22, 23], microlensing of highly magnified extragalactic stars [24], and
fast radio burst (FRB) timing [25]. Axion stars, generally expected to form at
the center of axion miniclusters, are more compact objects and can potentially
increase sensitivity to the QCD axion or axion-like particles. This motivates
us to study the properties of stars forming in these minihalos, specifically
the timescale for their formation and the mass and density of the stars
formed; for both the QCD axion and axion-like particles. With this axion star
mass function in hand, it is possible to determine the prospects for their
detection through various astrophysical observables.

Enabled by recent progress in our understanding of the condensation of wave-
like dark matter and the evolution of axion miniclusters, from both analytical
and numerical studies, we can quantitatively study the mass distribution of
axion stars and their compactness. There are still uncertainties that remain
to be sorted out by future numerical studies on the exact growth model of
axion stars or generic Bose stars. In the spirit of motivating indirect
searches of axion stars as a powerful probe to axion dark matter in the post-
inflationary scenario and stimulating new search strategies, we calculate the
mass function of axion stars using a semi-analytical method for both QCD axion
and axion-like particles, which can be further used to compare with other
observations.

In this work, we focus on the gravitational interactions of axion stars, which
are independent of the details of the axion self-couplings or couplings to SM
fields. One could also think of other observables such as radio signals,
coming from the axion-photon coupling, which we leave for future work. Our
proposal does not rely on the couplings between axions and standard model
particles but only on the gravitational interaction. Note that only post-
inflation axions may lead to the formation of axion miniclusters and axion
stars. Meanwhile, the post-inflationary scenario reduces the axion diffuse
flux that is the target of direct detection because most axions are bounded
inside axion miniclusters or axion stars, which requires better sensitivity
for direct detection [26, 27]. Therefore, our work opens an exciting new
window for indirect searches for the axion dark matter complementary to
existing direct detection efforts.

This work is organized as follows: In Sec. 2, a brief overview on the current
status of axion minicluster and axion star study is presented. In Sec. 3, we
perform our new calculation on the axion star mass function for post-inflation
axion models. In Sec. 4, we apply our axion star mass function to current
observables and obtain new limits on axion masses. Our bound on the axion mass
is $m_{a}>3.3\times 10^{-17}$ eV, which is currently the leading constraint
for post-inflation axions. In Sec. 5, we present the conclusion and discussion
of this work.

##  2 Overview of Axion Stars and Axion Minihalos

Axion stars, gravitational bound states of axion dark matter, have attracted
much interest and have been extensively studied in recent years. Generically,
there are two types of compact objects formed from axions: axion miniclusters
and axion stars. Axion miniclusters are gravitationally virialized objects
formed from density perturbations on small scales, which naturally arise from
isocurvature axion perturbations if the Peccei-Quinn symmetry is broken after
inflation. The formation of axion miniclusters starts as early as the time of
matter-radiation equality, after which the minihalos will subsequently merge
and grow in mass as time evolves. In the literature, the term “axion minihalo”
has been used interchangeably with “axion minicluster”, to emphasize the
similarity to the usual dark matter halos, especially regarding the density
profile [28]. We will refer to these virialized halo objects as axion
minihalos in the rest of this work.

On the other hand, axion stars are more compact than axion minihalos. Axion
stars are the ground state of axion field configurations sustained by a
balance between gravitational attraction and the repulsion of gradient
pressure; axion stars do not necessarily emit light like normal stars do.
Therefore, axion stars are expected to form through Bose-Einstein condensation
at the center of axion minihalos due to their large central density. There are
two branches of axion star solutions: the dilute branch has a balance between
gravity and axion kinetic pressure, and the dense branch has an unstable
equilibrium between the attractive forces of gravity and axion self-
interactions and the repulsion of kinetic pressure. Those two branches are
separated by a critical mass, above which dilute axion stars will collapse
into dense axion stars and subsequently explode by emitting relativistic
axions. Such axion star explosions can continue to occur recurrently if axion
accretion is fast enough, which are called recurrent axinovae, leading to the
depletion of axion dark matter [29]. In the scenario of QCD axion stars,
however, one can show the axion star mass is one order of magnitude below the
critical mass even considering the heaviest axion star formed in the current
Universe [30]. Therefore QCD axion stars are naturally stable dark compact
objects supported by the axion kinetic pressure. For ALPs, axion stars can
also be stabilized if their self-couplings are sufficiently small and do not
trigger axinovae.

There are extensive studies on the stable solutions of axion stars (or generic
Bose stars) [31, 32], the formation [33], and the collapse of axion stars [34,
35, 36]. However, the population analysis and detectable signatures are still
lacking. Recently, there have been calculations to determine the axion star
explosion rate [37, 29] and their encounter rates with neutron stars [38]. The
mass function and density of axion minihalos, as the environments for axion
star formation, are crucial to determining the formation rate and population
of axion stars. In this work, we adopt the state-of-art axion minihalo
population and axion star formation models to systematically study the axion
star mass function. In what follows, we will give an overview of the current
literature on axion star solutions, the evolution of axion minihalos, and the
formation of axion stars. Applying these theoretical tools, we will derive a
mass function of axion stars, which will be important inputs for indirect
searches of axion dark matter.

###  2.1 Stable configurations of axion stars

The stable axion-field configuration for the self-gravitating bound-state of
non-relativistic axions can be found by solving the Schrödinger-Poisson
equations (see [39, 40] for reviews). One could obtain generic solutions for
the normal modes of axion stars (or Bose stars) by solving the wave function
equations [41]. The ground state of a dense gas of axions is an axion star
surrounded by a less dense gas of free axions. Thus, one expects that, given
enough time, these objects will form in dense environments. The axion field or
ALP field $a$, is a pseudo-Goldstone boson of a global $U(1)_{PQ}$ symmetry
broken at a scale $f_{a}$. The axion mass term and self-couplings are
generated through the instanton potential. Since the $U(1)_{PQ}$ is anomalous
under the color gauge group (or a QCD-like group) after the confinement scale,
$\Lambda$, the axion potential will be generated through instanton effects,
leading to a light mass term for axions.

The axion potential takes the form [42]

$V(a)=-\frac{\Lambda^{4}}{c_{ud}}\sqrt{1-4c_{ud}\,\sin^{2}\left(\frac{a}{2f_{a}}\right)}~{}.$
(2.1)

For the QCD axion, $\Lambda\approx 80\,\rm MeV$ and $c_{ud}\approx
m_{u}m_{d}/(m_{u}+m_{d})^{2}\approx 0.2$. For ALPs, the scale associated with
the axion potential is a free parameter and $\Lambda$ can take on a wide range
of values. Even though the above expression specifically describes the QCD
axion potential, we expect ALPs to have a similar sine potential and very
similar phenomenology for the self-couplings. In what follows we will only be
interested in the mass and quartic interactions contained in the full
potential. Expanding Eq. 2.1 to quartic order one can find the mass term and
an attractive self-interaction term

$V(a)=\frac{1}{2}m_{a}^{2}a^{2}-\frac{\lambda}{4!}a^{4}~{}.$ (2.2)

For the QCD axion $m_{a}=\Lambda^{2}/f_{a}$ and
$\lambda=(1-3c_{ud})m_{a}^{2}/f_{a}^{2}$, that is the mass and self-coupling
are determined by the strong coupling scale of QCD ($\Lambda$) and the
breaking scale of PQ symmetry ($f_{a}$). For ALPs, where the $\Lambda$ is a
free parameter $m_{a}$ and $f_{a}$ can be treated as independent, we still
expect $\lambda\ll 1$.

If we expand the axion potential to higher-order terms, the next term to
appear is a repulsive self-coupling, at dimension-6. The analysis of the
initial collapse of axion stars is always governed by the quartic term and the
repulsive term will only become important during the further contraction of
axion stars once they pass the point of instability. These higher order terms
have effects on the remnants of axion star explosions. In this work, we focus
on the study of gravitational couplings of axions, which will dominate the
axion star configuration and the formation in most cases. The self-couplings
cause the instability of axion stars when they continuously accrete axion dark
matter and reach the critical mass. For the QCD-axion, this critical mass is
not reached in the age of the universe. For the case of an ALP it is always
possible to change the PQ breaking scale, $f_{a}$, such that the self-
couplings are too small to lead to critical star formation at the star masses
we are interested in.

In the following, we sketch how to determine the stable axion field
configurations. In the most generic case, the Klein-Gordon equation of axion
field in curved spacetime can be written as

$\displaystyle\frac{1}{\sqrt{-g}}\partial_{\mu}(\sqrt{-g}g^{\mu\nu}\partial_{\nu}a)=-V^{\prime}(a).$
(2.3)

It is often not required to work with the full theory of General Relativity.
The axion stars, especially the stable ones on the dilute branch are
nonrelativistic objects as long as $f_{a}\ll M_{\rm pl}$. Therefore we can
always simplify the calculation by studying the self-gravitating axion stars
in the following metric

$\displaystyle
ds^{2}=g_{\mu\nu}dx^{\mu}dx^{\nu}=(1+2\Phi)dt^{2}-(1-2\Phi)dr^{2}-r^{2}d\Omega^{2},$
(2.4)

where $\Phi$ is the gravitational potential given by the self-gravity of axion
stars and $d\Omega$ is the differential solid angle. In the non-relativistic
limit, we can also use single harmonics to approximate the axion field, which
gives the following expression

$\displaystyle a(r,t)=f_{a}\,\Theta(r)\cos(\omega_{a}t),$ (2.5)

where $\Theta(r)$ is the spatial-dependent part of the axion field and
$\omega_{a}$ is the energy of the axions in the star. For the non-relativistic
stars $\omega_{a}=m_{a}+p^{2}/2m_{a}\approx m_{a}$. If we neglect the axion
self-interaction and higher order terms, Eq. 2.3 becomes

$\displaystyle\boldsymbol{\nabla}^{2}\Theta=2\left(m^{2}_{a}\Phi+\frac{m^{2}_{a}-\omega^{2}_{a}}{2}\right)\Theta,$
(2.6)

where the spatial Laplace operator can be written as
$\boldsymbol{\nabla}^{2}=\frac{\partial^{2}}{\partial
r^{2}}+\frac{2}{r}\frac{\partial}{\partial r}$ since we study the ground state
configuration of axions that is spherically symmetric. In the Newtonian limit,
the gravitational potential can be simply determined from the energy density
via the Poisson equation

$\displaystyle\boldsymbol{\nabla}^{2}\Phi=4\pi G\rho_{\star},$ (2.7)

with $\rho_{\star}$ being the density of axion stars that contains the rest
mass, kinetic energy, gradient energy, and self-energy, which takes the form

$\displaystyle\rho_{\star}=\frac{m_{a}^{2}f_{a}^{2}\Theta^{2}}{2}+\frac{f_{a}^{2}}{4}\left(\frac{d\Theta}{dr}\right)^{2}-\frac{\lambda
f_{a}^{4}\Theta^{4}}{64}.$ (2.8)

Where we have taken the time average of the axion field. Without fully solving
these equations numerically, one could alternatively use a Gaussian ansatz to
approximate the field profile [43, 44, 36], which can demonstrate the
competing effects driving the physics of stabilizing the axion stars [32].
With this method, the energy component of an axion star of mass $M_{\star}$
and radius $R_{\star}$ can be determined quantitatively as

$E_{\star}=-\frac{G_{N}M_{\star}^{2}}{R_{\star}}+c_{1}\frac{M_{\star}}{2\,m_{a}^{2}R_{\star}^{2}}-c_{2}\frac{\lambda
M_{\star}^{2}}{12\,m_{a}^{4}R_{\star}^{3}}~{}.$ (2.9)

The sum of these terms gives the total energy, with individual components
corresponding to the self-energy term from gravity, the gradient energy, and
the self-energy from the axion potential. The numerical coefficients, $c_{i}$,
are solved numerically for the full system [45, 46, 32] to be $c_{1}=9.9$,
$c_{2}=0.85$. As expected, the formalism of axion star energy shows that the
axion star is balanced by gravity and kinetic pressure (or gradient energy).
By adding mass to the axion star, one increases the gravitational self-
interaction, which requires more kinetic pressure to balance the gravity,
leading to a smaller axion star radius since the gradient energy goes like
$1/R_{\star}^{2}$ while gravitational energy goes like $1/R_{\star}$. However,
such balance will eventually be broken down by the quartic self-coupling, with
an energy component that scales as $1/R_{\star}^{3}$. At the critical point,
the density of axion stars is so high that self-couplings become important,
and kinetic pressure cannot balance it by further decreasing the radius. The
mass-radius relation for axion stars can be found by minimizing $E_{\star}$,
which corresponds to the relaxation to the ground state of axion
configurations and gives the following radius

$R_{\star}^{\pm}=\frac{c_{1}}{2\,G_{N}M_{\star}m_{a}^{2}}\left(1\pm\sqrt{1-\frac{c_{2}}{c_{1}^{2}}\lambda
G_{N}M_{\star}^{2}}\right)~{}.$ (2.10)

The $R_{\star}^{+}$ root corresponds to the dilute branch of axion stars,
which is the axion field configuration stabilized by kinetic pressure and
gravity while the self-interaction is negligible. Since the quartic self-
couplings can be ignored in the dilute branch, the term related to $\lambda$
can be dropped, leading to a mass-radius relation

$R_{\star}^{+}=9.9\frac{M_{\mathrm{pl}}^{2}}{m_{a}^{2}M_{\star}}\approx
2.6\times 10^{5}\,{\rm km}\left(\frac{10^{-6}\rm
eV}{m_{a}}\right)^{2}\left(\frac{10^{-12}M_{\odot}}{M_{\star}}\right)~{}.$
(2.11)

In general, a lighter axion mass $m_{a}$ corresponds to a larger axion star
due to the large de Broglie wavelength. On the other hand, a heavier axion
star has a smaller radius, because it requires more kinetic pressure to
balance the gravity. Overall, from Eq. 2.11 we see that the density of a star
grows as the fourth power of its mass, $\rho_{\star}\sim
M_{\star}^{4}(G_{N}m_{a}^{2})^{3}$. This will have important implications when
we study the gravitational detectability of QCD axion stars. Since more
compact objects can be probed with observations such as microlensing, one
needs to form a sufficiently heavy axion star to reach a high density for
detection. Therefore, the evolution and growth of axion stars will be a
crucial part of the story (Unfortunately, as we will discuss later in Fig. 1,
the QCD axion stars will not grow to a mass that will be compact enough for
the current sensitivity of microlensing). As the mass of the axion star grows,
Eq. 2.10 shows that at some point a sensible solution can no longer be
obtained. This is due to the breakdown of the balance between gravity and
gradient pressure caused by the non-negligible addition of attractive self-
couplings, leading to the collapse of axion stars. The maximal mass that still
corresponds to a stabilized axion star is

$M_{\star}^{\textrm{max}}=\frac{10.7}{\sqrt{\lambda}}M_{\mathrm{pl}}\approx
1.16\times 10^{-10}M_{\odot}\left(\frac{\rm
10^{-6}eV}{m_{a}}\right)\left(\frac{f_{a}}{10^{12}\rm GeV}\right)~{}.$ (2.12)

The two solutions (2.10) merge at this critical mass (The $R_{\star}^{-}$ root
is sometimes called the critical branch, which is metastable and not relevant
to our study). After the collapse of critical mass stars, we will obtain the
dense branch of axion stars which will soon emit relativistic axions [47, 48,
35] and are thus cosmologically unstable, with a lifetime $\sim
10^{3}m_{a}^{-1}$ [32]. In this work, we will not consider strong self-
couplings that will destabilize axion stars. This corresponds to a large
$f_{a}$ that is not excluded by Ref. [29, 30] that considers the collapse of
axion stars and their observational constraints. In this work, we will focus
on the stable axion stars and study their gravitational detectabilities.

###  2.2 An overview on axion minihalos

Axion minihalos are natural products of axion cosmology in the post-
inflationary scenario where Peccei-Quinn symmetry is broken after inflation
[16, 17, 18]. The axion field will randomly take different values at different
horizon patches when its potential is flat. During the QCD phase transition,
axions acquire their mass, converting the inhomogeneity of field values to
matter density fluctuations. This leads to order unity matter perturbations at
the scale corresponding to the size of horizon patches when axions obtain
their mass. These large perturbations subsequently and form axion minihalos
at, or shortly after, matter-radiation equality. The axion isocurvature
fluctuations at superhorizon scales follow a white noise power spectrum
because different horizon patches were not in casual contact so there is no
correlation. Furthermore, because many horizon patches contribute to
perturbations at superhorizon scales, the spectrum of fluctuations is Gaussian
and the power spectrum fully describes the fluctuations. Note that these
superhorizon modes will enter the horizon before matter-radiation equality and
grow linearly with the scale factor. From the white-noise power spectrum, one
can study the population of axion minihalos at different cosmic times.

The cosmological evolution of axion minihalos after its initial formation has
been studied using semi-analytic techniques [49, 50, 19, 51, 52] as well as
full numerical simulations [53, 54, 28, 26, 55], which suggest $\sim 75\%$ of
the axion dark matter will be bound into axion minihalos in the late Universe
due to the large axion perturbations. Even though they are abundant, it is
very challenging to observe these axion minihalos. Some ideas might involve
probing these structures via their photon signals such as the stimulated decay
of axions [56, 57, 58, 59, 60, 61, 62, 63, 64, 65], and axion-photon
conversion in neutron stars [66, 67, 68, 69]. Other proposals considered the
lensing or photon signatures from axion stars [70, 71, 72, 73, 74, 75, 76, 77,
78, 79]. Our work will focus on the calculation of the axion star mass
function and its gravitational effects, which will have important implications
for the rate of detectable axion star events.

According to simulations of string dynamics [80, 81, 82], the dimensionless
matter power spectrum $\Delta^{2}(k)$, after the axion string network has
collapsed, follows an approximate (white noise) power law at large scales
before becoming approximately constant at short distance scales and then
ultimately decreasing at very small scales. The exact form of the power
spectrum is still under active study. Here we choose to concentrate on the
power law part of the power spectrum and to remain agnostic to the exact form
once this power law breaks. We do this by modeling the spectrum as a pure
white noise power law at all scales but place a lower bound on halo masses so
that we do not probe the non-white noise region of the power spectrum. Thus,
the dimensionless matter power spectrum $\Delta^{2}(k)$ takes the form

$\Delta^{2}(k)\equiv\frac{k^{3}}{2\pi^{2}}P(k)=A_{0}\left(\frac{k}{k_{0}}\right)^{3}\,,$
(2.13)

where $k_{0}$, which is the comoving wavenumber where the matter power
spectrum in full simulations becomes constant, and $A_{0}$ is the amplitude at
this scale, which is also the peak amplitude of $\Delta^{2}(k)$. Here
$P(k)\equiv V^{-1}|\tilde{\delta}(\mathbf{k})|^{2}$ is the matter power
spectrum, $V$ is the volume, $\tilde{\delta}$ is the Fourier transform of the
linear dark matter overdensity. For the QCD axion $A_{0}\approx 0.4$, while
for ALPs $A_{0}\approx 0.3$ [82]. The scale where the power law breaks,
$k_{0}$, is related to the comoving horizon scale at the time of axion
oscillation, $k_{\rm osc}\equiv a_{\rm osc}{H}_{\rm osc}$. For the case of the
QCD axion $k_{0}\approx 3k_{\rm osc}$, while for ALPs $k_{0}\approx k_{\rm
osc}$. Associated with the scale $k_{0}$ is a mass, which we denote the
characteristic halo mass,

$M_{h,c}=\frac{4\pi}{3}A_{\rm 0}\left(\frac{1}{k_{\rm
0}}\right)^{3}\bar{\rho}_{a,0}~{}.$ (2.14)

We focus only on halos whose mass is greater than $M_{h,c}$. Halos with
$M_{h}<M_{h,c}$ will correspond to a higher wavenumber, $k>k_{0}$, where the
spectrum of perturbation is less well understood. These high wavenumbers do
not contribute much mass to the total minihalo distribution.

The mass function of axion minihalos, especially at the high mass end, which
constitutes the majority of the mass, is well studied with N-body simulations
that also agree with a calibrated analytical formula [28]. A semi-analytical
model by Press $\&$ Schechter [83] gives the following halo mass function:

$\frac{df_{h}}{d{\rm ln}M_{h}}=\sqrt{\frac{\nu}{2\pi}}{\rm
exp}(-\nu/2)\frac{d\,{\rm ln}\nu}{d\,{\rm ln}M_{h}}\,,$ (2.15)

where $M_{h}$ is the axion minicluster(or minihalo) mass and $df_{h}/d{\rm
ln}M_{h}$ is the mass fraction of axion minihalos per logarithmic mass i.e.
$f_{h}$ is the mass fraction of dark matter bound in minihalos. We choose
$df_{h}/d{\rm ln}M_{h}$ to present the mass distribution of the population of
axion minihalos, as it can characterize the most significant mass ranges of
these objects that contain the majority of the total mass when considering
logarithmic mass bins. $f_{h}=M_{h}n_{h}/\bar{\rho}_{a,0}$ is the mass density
fraction of axion in minihalos. $\nu$ is defined as:

$\nu\equiv\frac{\delta_{c}^{2}}{\sigma^{2}(M_{h})D(z)^{2}}\,,$ (2.16)

where $\delta_{c}=1.686$ is the critical density required for spherical
collapse and $D(z)$ is the growth function which is normalized to be
$D(z_{i})=1$ where $z_{i}$ is the initial redshift in radiation dominated era.
For subhorizon modes, the growth function is described by the
M$\rm\acute{e}$sz$\rm\acute{a}$ros equation which gives a solution

$D(z)=1+\frac{3}{2}\left(\frac{1+z_{\rm eq}}{1+z}\right)~{}.$ (2.17)

We are ignoring the effects of dark energy since axion star formation occurs
early when dark energy is a subdominant fraction of the universe. The variance
of the axion perturbations, $\sigma^{2}(M_{h})$, in the initial density
fluctuations with a tophat filter of scale $R=(3M_{h}/4\pi\bar{\rho})^{1/3}$,
can be calculated as:

$\sigma^{2}(M_{h})\equiv\int\frac{dk}{k}\frac{k^{3}P(k)}{2\pi^{2}}\left|\widetilde{W}(kR)\right|^{2},$
(2.18)

where $\widetilde{W}(x)=(3/x^{3})[{\rm sin}(x)-x{\rm cos}(x)]$ is the
spherical top-hat window function. There are other reasonable choices for a
window function, e.g. sharp $k$-space filter, which give qualitatively similar
results. The variance of the white-noise power spectrum from axion
perturbations is

$\sigma(M_{h})=\sqrt{\frac{3\pi}{2}\frac{M_{h,c}}{M_{h}}}\,.$ (2.19)

As we see in the above expression, the only relevant parameter for
$\sigma(M_{h})$ is the dimensionless quantity $x_{h}\equiv M_{h}/M_{h,c}$,
which will give the full description of the dimensionless mass function in Eq.
2.15. At $\nu<1$, the exponential part of the mass function can be ignored and
the mass function exhibits a power law $df_{h}/d{\rm ln}M_{h}\propto
x_{h}^{1/2}$. As discussed before, we cut off the mass ranges $M_{h}<M_{h,c}$
(or $x_{h}<1$) to filter out the contributions from high wave numbers; we only
consider $x_{h}\geq 1$.

Therefore the mass distribution of axion minihalos is peaked at larger masses,
which suggests that the observable is sensitive to the axion perturbations at
relatively large scales. Note that the above mass function includes the white-
noise power spectrum only, which only applies at $z\gtrsim 20$ when large-
scale structures have not formed yet. Afterward, the infall of axion minihalos
into large-scale structures will stop the evolution of the mass function, as
extensively discussed in Ref. [28]. In this work, we focus on studying axion
minihalos around matter-radiation equality, which maximizes the axion star
formation rate. Therefore, the above mass function that describes the
minihalos at early times should suffice for our purpose.

The density of axion minihalos is also needed to compute the formation
timescale of axion stars. Axion minihalos are highly concentrated
substructures, with a Navarro-Frenk-White (NFW) [84] density profile as tested
in N-body simulations [28]. The size and density of an NFW halo can be
characterized by the scale radius and the scale density, which will ultimately
determine the axion star formation rate in axion minihalos. As computed in
N-body simulations [28], the scale radius of axion minihalos formed from a
white-noise power spectrum is given by:

$r_{s}\approx 3125{\rm
AU}\left(\frac{M_{h,c}}{10^{-12}M_{\odot}}\right)^{-1/2}\left(\frac{M_{h}}{10^{-6}M_{\odot}}\right)^{5/6},$
(2.20)

where $M_{h}$ is the axion minihalo mass. Given the scale radius and the
minihalo mass, one can determine $\rho_{s}$ by:

$\rho_{s}=\frac{M_{h}}{4\pi r_{s}^{3}(\log(1+c)-c/(1+c))}\,,$ (2.21)

where $c$ is the halo concentration number, which is defined by the ratio of
virial radius to scale radius, $c=r_{\rm vir}/r_{s}$. Simulations show that
the concentration of axion minihalos can be expressed as $c\approx 1.4\times
10^{4}\,(1+z)^{-1}\sqrt{M_{h,c}/M_{h}}$ [28], which agrees well with an
analytical model that the concentration of axion minihalos is roughly 4 at
birth ($\sigma(M_{h})=\delta_{c}$). At matter-radiation equality, most axion
minihalos born at that time should have a concentration close to 4. When we
calculate the condensation timescale of axion stars in realistic axion
minihalos later, the central density and velocity will be needed. For
conservative estimates, we can use the scale density as a benchmark for the
central density. The characteristic velocity can be estimated as the circular
velocity at scale radius

$v_{\rm c}=\sqrt{2\pi Gr_{s}^{2}\rho_{s}\left(\log 4-1\right)}~{}.$ (2.22)

###  2.3 Formation of axion stars in axion minihalos

Since the axion star configuration, as discussed earlier in this section, is
an energetically favored state of axion fields, the axion waves prefer to
reach the axion star configuration as long as there is enough time for
relaxation. Such processes can be understood as Bose-Einstein condensation of
Bosonic dark matter. Therefore, the relevant timescale of axion star formation
is the scattering timescale between axions as wave-like particles, which
approximately gives the condensation timescale of axion stars in simulations
of axion fields [33, 85, 86, 87]

$\tau\sim\left(f_{\mathrm{BE}}n\sigma v\right)^{-1}\,,$ (2.23)

where $f_{\mathrm{BE}}=6\pi^{2}n(m_{a}v)^{-3}$ is the phase density of axions,
with $n$ the typical axion number density and $v$ is the axion speed. The
axion wave scattering cross-section is $\sigma$, this scattering can be due to
axion self-couplings or gravitational interactions.

This scattering of axion waves is Bose-enhanced due to the large phase space
density. Therefore, a light axion mass or a small velocity will greatly
enhance the formation rate of axion stars. Quartic self-interactions of axions
can further enhance the rate but also lead to the collapse of axion stars and
interesting constraints on the attractive self-couplings [29]. In this work,
we focus on the gravitational interactions of axions, which gives the minimal
rate of axion star formation. All axion stars formed with negligible self-
couplings will be cosmologically stable objects. The cross-section of
gravitational scattering is
$\sigma_{\mathrm{gr}}=8\pi(G_{N}m_{a}v^{-2})^{2}\log(m_{a}vr_{\rm vir})$,
where the Coulomb logarithm applies a cutoff at the size of the axion
minihalo, $r_{\rm vir}$. Therefore, the gravitational condensation timescale
is

$\tau=\frac{b}{48\pi^{3}}\frac{m_{a}v^{6}}{G_{N}^{2}n^{2}\log\left(m_{a}vr_{\rm
vir}\right)}\,.$ (2.24)

The parameter $b\sim\mathcal{O}(1)$ is a numerical coefficient determined by
numerical simulations [88]. $r_{\rm vir}=c\,r_{s}$ is the virial radius of the
minihalo that provides the cutoff scale. We focus on the axion star formation
around matter-radiation equality, and we choose $c=4$ similar to the choice in
the previous subsection for a simple estimation of the logarithmic factor.
These formulae can calculate the relaxation scale of axion gas with a given
velocity distribution and agree well with the results found in numerical
simulations [33, 85, 86, 89].

The gravitational condensation time, which is proportional to
$v^{6}/\rho^{2}$, is dramatically shorter at a smaller radius in the halo
center. In this work, we calculate the characteristic density and speed at the
scale radius for conservative estimations. Since the star formation rate peaks
at the center, we generically expect the axion star to form at the halo
center. While multiple stars may form in a halo when the self-interaction is
strong enough and axion gas fragments to stars in different regions of one
minihalo, we do not consider that possibility and focus on the scenario of
only one axion star per axion minihalo. The story of the axion star formation
rate is not yet finished since we have not determined the relevant mass scales
of axion stars formed after the relaxation. It was recently found the growth
of Bose stars in a halo environment is self-similar [90], and the growth
behavior can be described by

$\frac{(1+x_{\star}^{3}/\epsilon^{2})^{3}}{(1-x_{\star})^{5}}\approx\frac{t+0.1\tau}{1.1\tau},$
(2.25)

where we define $x_{\star}(t)\equiv M_{\star}(t)/M_{h}$ with the axion star
mass $M_{\star}$ and $x_{\star}=0$ at $t=\tau$, with $\tau$ the timescale for
the minihalo system to become thermalized. The quantity $\epsilon\approx
3M_{\star,c}/M_{h}\approx 10vM_{\rm pl}^{2}/(m_{a}M_{h})$ is a combination of
total mass and energy in the minihalo system, which is a measure of the star-
halo mass relation. $M_{\star,c}$ is the characteristic mass scale of the
axion star and can be well approximated by equating the virial velocity of the
axion star and the halo velocity, which gives [33, 85]

$M_{\star,c}\approx\frac{3\,v}{G_{N}m_{a}}\approx
3\rho_{a}^{1/6}G_{N}^{-1/2}m_{a}^{-1}M_{h}^{1/3}~{},$ (2.26)

where $v$ is the velocity of axion waves in axion minihalos, $\rho_{a}$ is the
scale density of axion minihalos, and $M_{h}$ is the axion minihalo mass.
Sometimes the characteristic mass of axion stars $M_{\star,c}$ is also called
the “saturation” mass since the growth slows down after this mass is reached
according to numerical simulations. However, the mass growth will continue
even after $x_{\star}\sim\epsilon$ in this self-similar model discussed above.
The power law is indeed different once $x_{\star}^{3}/\epsilon^{2}\gg 1$ is
satisfied. As $x_{\star}$ approaches 1, it takes more and more time for the
star growth, leading to the slow down of axion star formation. The above
expressions can be used to calculate the mass of the axion stars once we have
the axion minihalo properties.

##  3 Axion Star Mass Functions

While the formation of axion minihalos, or the formation of axion stars within
certain backgrounds, has been studied, the current time mass distribution of
axion stars is not well-known. Determining an exact axion star mass function
requires extensive numerical study taking into account physics on widely
separated scales and wave effects, but we provide an analytical estimate of
the axion star mass distribution using several conservative assumptions. As a
simple estimate, we suggest taking the axion minihalo distribution at the time
of matter-radiation equality ($z\approx 3400$) and considering axion stars
formed within each minihalo over one Hubble time from the time of equality.

The first minihalo starts to form roughly at matter-radiation equality because
the Hubble time at equality is approximately the same as the minihalo
dynamical timescale needed to form a virialized minihalo. After the equality,
axion stars can form and grow, but axion minihalos also grow and merge. As
axion minihalos become heavier, the axion star formation time increases faster
than the Hubble time, causing the growth of axion stars to cease. Also, we
only consider one axion star per minihalo because they typically form in the
center of minihalos, and if more than one stars are formed they quickly merge.
Thus, we focus on the formation around matter-radiation equality, and the late
evolution of axion stars such as mergers, which requires more extensive
numerical studies, is not considered in this work. Our calculation
underestimates the mass of individual stars, but we expect the total mass
fraction in axion stars is well captured. This simple estimate already
suggests that an $\mathcal{O}(1)$ fraction of axion dark matter forms axion
stars, which can lead to interesting observational constraints and potential
signals.

Once we have the axion minihalo properties, the axion star mass formed within
the minihalo can be determined. We can calculate the axion star mass for a
given axion mass $m_{a}$ and minihalo mass $M_{h}$, at a particular time, by
solving Eq. 2.25. Which gives

$M_{\star}(m_{a},M_{h})=M_{h}\,x_{\star}\left(\frac{1}{H(z_{\rm
eq})}\right)~{}.$ (3.1)

Recall that $x_{\star}$ depends upon the typical density and size of the
minihalo as well as the speeds of axions in the minihalo. All of these are
determined, through Eqs. 2.20, 2.21, and 2.22, by the minihalo mass (along
with the parameters associated with the power spectrum, $k_{0}$ and $A_{0}$).
Since we limit ourselves to the case where there is at most one axion star per
minihalo the axion star mass fraction is related to the axion minihalo mass
fraction $df_{\star}/dM_{h}=(M_{\star}/M_{h})df_{h}/dM_{h}$. Furthermore, from
the star-minihalo mass relation, the axion minihalo population described in
Eq. 2.15 can be used to determine the axion star mass function by the
following expression

$\frac{df_{\star}}{dM_{\star}}=\frac{M_{\star}}{M_{h}}\frac{dM_{h}}{dM_{\star}}\frac{df_{h}}{dM_{h}}\,,$
(3.2)

where $f_{\star}$ is the mass fraction of axions in axion stars to total
ambient axions. We will show the mass distributions of axion stars for both
QCD axions and ALPs, in the next two subsections.

###  3.1 QCD Axion Stars

Figure 1:  This plot presents the axion star mass and density as a function of
axion minihalo mass $M_{h}$ and axion mass $m_{a}$. We scale the $y$-axis to
use the mass ratio $x_{h}=M_{h}/M_{h,c}\geq 1$ to represent the minihalo mass.
The black contours show different axion star masses
($10^{-14}M_{\odot}\lesssim M_{\star}\lesssim 10^{-11}M_{\odot}$), and the red
contours show different densities of axion stars
($10^{4}\,M_{\odot}/\textrm{pc}^{-3}\lesssim\rho_{\star}\lesssim
10^{6}\,M_{\odot}/\textrm{pc}^{-3}$). The blue dashed line indicates the pivot
minihalo mass $x_{h,p}$, which gives the maximum axion star mass for a given
axion mass. The grey region does not lead to a sufficient axion star formation
rate because the condensation timescale is too long compared to the Hubble
time at matter-radiation equality. The ratio of the star to minihalo mass are
shown as blue-shaded regions. A smaller $x_{h}$ corresponds to a higher halo
density and thus a larger axion star formation rate. Therefore, the axion star
to halo mass ratio $M_{\star}/M_{h}$ increases with a smaller $x_{h}$. We
focus on the mass range $m_{a}<10^{-5}\,\rm eV$ since we mainly study halos
with $x_{h}>1$ that will dominate the population of axion minihalos. It is
possible to have $x_{h}<1$, especially considering the effects from axion
strings [91], which will further enhance the axion star formation for larger
$m_{a}$.  Figure 2: This plot presents the mass function of QCD axion stars
for various axion masses $m_{a}$, which correspond to different typical values
of axion star masses. The solid curves represent the contribution from
$x_{h}<x_{h,p}$ while the dashed curves come from larger axion minihalos with
$x_{h}>x_{h,p}$. Note that axion stars formed from smaller axion minihalos
($x_{h}<x_{h,p}$) will dominate the mass function. A smaller $m_{a}$ will lead
to a larger dark matter fraction in axion stars, $f_{\star}$, due to the
shorter gravitational condensation timescale.

It is crucial to determine the QCD axion minihalo mass before obtaining the
mass function of QCD axion stars, which will be discussed in this subsection.
The time when axions start to oscillate is determined from $m_{a}(T_{\rm
osc})=1.6H(T_{\rm osc})$, which will give the characteristic mass, $M_{h,c}$,
of axion minihalos. Since the mass term of the QCD axion is induced by the
temperature-dependent QCD instanton potential, it is suppressed at
temperatures higher than the QCD confinement scale. Therefore, the oscillation
of the QCD axion is generally delayed compared to the naive expectation from
its zero temperature mass. The temperature-dependent QCD axion mass is given
by $m_{a}(T)\sim m_{a}(0)(T/\Lambda_{\rm QCD})^{-4}$ for $T\gtrsim\Lambda_{\rm
QCD}$ [92]. For the QCD axion mass range we will consider, $m_{a}>10^{-8}$ eV,
the oscillation temperature is higher than the QCD phase transition scale,
$T_{\rm osc}\gtrsim\Lambda_{\rm QCD}$ and the temperature dependent mass must
be considered. For smaller axion masses, oscillation will occur at the zero
temperature mass. We consider this possibility, over a broader range of masses
in Sec. 3.2.

The mass range for the QCD axion is motivated by the studies on the relic
abundance, with lighter axion masses predicting a larger relic density since
the oscillation starts later and the axion density is thus less diluted. The
vacuum misalignment mechanism predicts a mass of $m_{a}\sim 10^{-5}$ eV for
the correct relic density. Extra contributions to the axion density from
string emissions will predict a relatively larger relic density and thus
prefer larger masses [93, 14, 13, 94, 95, 96]. However, non-standard thermal
history such as an early-matter-dominated epoch that reheats after the QCD era
can inject extra entropy into the Universe, diluting the axion density. We
consider a broad QCD axion mass ranges
$10^{-8}\,\mathrm{eV}<m_{a}<10^{-5}\,\mathrm{eV}$ in this work.

Using the full temperature-dependent mass of the QCD axion, the oscillation
scale $k_{\rm osc}$ is [97, 80]

$k_{\rm osc}\approx
14~{}\textrm{pc}^{-1}\left(\frac{m_{a}}{\mu\mathrm{eV}}\right)^{0.17}\,.$
(3.3)

Given that $k_{0}\sim 3k_{\rm osc}$ and $A_{0}\sim 0.4$, as discussed earlier,
the characteristic mass of axion minihalos Eq. 2.14 is

$M_{h,c}\approx 7.0\times
10^{-13}M_{\odot}\left(\frac{\mu\mathrm{eV}}{m_{a}}\right)^{0.51}\,.$ (3.4)

As for the halo mass function (see Sec. 2.2) it is convenient to introduce the
dimensionless parameter $x_{h}\equiv M_{h}/M_{h,c}>1$ to represent the axion
minihalo mass in a $M_{h,c}$ independent way. Then combining Eq. 2.20, Eq.
2.21, and Eq. 2.26, we can get the characteristic axion star mass

$M_{\star,c}\approx 1.4\times
10^{-13}M_{\odot}\,x_{h}^{1/12}\left(\frac{\mu\mathrm{eV}}{m_{a}}\right)^{1.17}\,.$
(3.5)

With the determination of the axion minihalo mass, we can obtain the formation
timescale (Eq. 2.24) for QCD axion stars:

$\tau H(z_{\rm
eq})\approx\frac{0.083x_{h}^{7/2}(m_{a}/\mu\mathrm{eV})^{1.98}}{\log[29x_{h}^{11/12}(m_{a}/\mu\mathrm{eV})^{0.66}]}$
(3.6)

With these expressions, the QCD axion star mass present in a halo shortly
after matter-radiation equality can be computed, from Eq. 3.1.

In Fig. 1, we present the properties of QCD axion stars formed in a given
axion minihalo with a mass of $M_{h}$ within the Hubble time at matter-
radiation equality. We also show the contours of different masses and
densities of axion stars as a function of $m_{a}$ and $M_{h}/M_{h,c}$, where
$M_{h,c}$ is the characteristic mass of initial axion minihalos. For larger
minihalos, the axion star mass decreases due to a more suppressed formation
rate and the corresponding contribution to the total mass is also suppressed.
Meanwhile, lighter axion minihalos have a larger axion star formation rate due
to the higher density, leading to a larger axion star to minihalo mass ratio
$M_{\star}/M_{h}$. The diagonal dashed blue line indicates the pivot minihalo
mass, $x_{h,p}$ which is the $x_{p}$, or equivalently minihalo mass, that
maximizes the axion star mass formed in the minihalo for a given axion mass.
The grey region presents the parameter space without any axion star formation
around matter-radiation equality, which happens for larger axion minihalos or
heavier axion masses. An important note is that we focus on the white-noise
part of the axion power spectrum, which allows us to derive all the halo
properties analytically. However, axion minihalos with $x_{h}<1$ are predicted
by simulations of axion string dynamics, which can potentially enhance the
axion star formation rate [91]. Note that such axion minihalos usually do not
contribute a significant mass fraction of the total axion dark matter
population, which is dominated by the axion minihalos at the high mass end
[28]. Therefore, as mentioned earlier, we focus on the white-noise part of the
power spectrum which captures the most important physics in a simple analytic
treatment.

The mass function of QCD axion stars is plotted in Fig. 2 for various axion
masses. We labeled two contributions to the axion star mass function
separately with solid and dashed curves, which come from lighter axion
minihalos ($x_{h}<x_{h,p}$) and heavier axion minihalos ($x_{h}>x_{h,p}$)
respectively. The actual mass function is the sum of these, but the lighter
axion minihalo part always dominates, if it exists. Note that there is no
solid line for $m_{a}=10^{-6}~{}\mathrm{eV}$ because $x_{h,p}<1$ for this mass
as shown in Fig. 1. At the pivot minihalo mass, the mass function diverges
because $dM_{h}/dM_{\star}$ diverges, as can be seen from the shape of the
black $M_{\star}$ contours in Fig. 1. However, physical observables depend on
the integration of the mass function, which is finite. As we see in the plot,
the total fraction of axion dark matter in axion stars, $f_{\star}$,
monotonically decreases with axion mass. This behavior is related to the
temperature-dependent mass of the QCD axion. The axion minihalo mass is
determined by the critical time when axions start to oscillate and behave as
nonrelativistic particles at $m_{a}(T)\approx H$. Usually, the QCD axion
oscillates before the QCD epoch, where the axion mass is smaller than its zero
temperature value after the QCD phase transition. This relatively late
oscillation leads to a larger axion minihalo mass and a smaller axion star
formation rate. As we see in Fig. 4, the mass fraction in ALP stars is
invariant with different axion masses since they start the oscillation at
zero-temperature mass. As we discussed before, the mass function of axion
stars is typically peaked at heavier masses as shown in Fig. 2.

###  3.2 ALP Stars

Figure 3: This figure has the same information as in Fig. 1 but focuses on
axion-like particles here. Again there is no axion star formation in the grey
region because axion minihalos with lower masses have a larger density and
thus a higher axion star formation rate. The characteristic axion minihalo
mass $M_{h,c}$ is determined by the mass contained within the horizon when
$m_{a}=H$, which is the canonical scenario for the vacuum misalignment
mechanism. The red curves are the contours of axion star densities, which are
nearly horizontal up to changes on $g_{\star}$. Note that the critical density
of the Universe today is about $10^{-7}M_{\odot}/\rm pc^{3}$ while the axion
stars have densities larger than $10^{3}M_{\odot}/\rm pc^{3}$. The black
curves are the contours of axion star masses, which are approximately vertical
and mostly determined by axion masses. A smaller $m_{a}$, corresponding to a
longer de Broglie wavelength, will give a larger axion star mass. Therefore,
axion stars with large masses and high densities can lead to interesting
observational constraints. Figure 4: The mass function of ALP stars with the
same scheme in Fig. 2. The mass function is self-similar across different
axion masses. This self-similarity differs from that in QCD axion stars where
the axion mass is temperature-dependent around the initial axion oscillation
time $m\sim H$. Still, the solid curves dominate the total mass of axion
stars, suggesting that the low mass end of axion minihalos contributes most to
the axion star population. A large fraction of axion dark matter,
$f_{\star}\gtrsim 0.4$, is in the form of axion stars even though we only
include the formation around matter-radiation equality.

In this subsection, we consider an ALP that starts to oscillate at its zero
temperature mass. This determines the oscillation temperature through [98]

$m_{a}=q_{0}H(T_{\rm
osc})=q_{0}\sqrt{\frac{\pi^{2}}{90}}\frac{g_{*}^{1/2}(T_{\rm osc})T_{\rm
osc}^{2}}{M_{\rm pl}},$ (3.7)

with $q_{0}=1.6$ and $T_{\rm osc}$ is the temperature of the Universe at the
onset of axion oscillation. Note that we assume a standard cosmology where the
Universe was dominated by radiation at early times, one could also consider
alternative cosmological scenarios such as early matter domination [99, 98].
Since entropy is conserved after axion oscillation, we can relate the scale
factor at oscillation to the current temperature via $g_{*S}(T_{\rm
osc})a_{\rm osc}^{3}T_{\rm osc}^{3}=g_{*S}(T_{0})T_{0}^{3}$, where $T_{0}$ is
the current CMB temperature. Therefore we can obtain the follow expression for
$k_{\rm osc}$

$k_{\rm osc}=\left(\frac{g_{*S}(T_{0})}{g_{*S}(T_{\rm
osc})}\right)^{1/3}\frac{T_{0}}{T_{\rm osc}}\frac{m_{a}}{q_{0}}.$ (3.8)

With this value of $k_{\text{osc}}$, and recalling that for an ALP
$k_{0}=k_{\text{osc}}$ and $A_{0}\approx 0.3$, we can now repeat the process
in the previous subsection to calculate $M_{h,c}$, $M_{\star,c}$, and $\tau$
for an ALP. This results in

$\displaystyle M_{h,c}$ $\displaystyle\approx$ $\displaystyle
14~{}M_{\odot}\left(\frac{10^{-16}~{}\mathrm{eV}}{m_{a}}\right)^{3/2}\left(\frac{g_{*}(T_{\rm
osc})}{g_{*}(T_{0})}\right)^{1/4}~{},$ (3.9) $\displaystyle M_{\star,c}$
$\displaystyle\approx$ $\displaystyle
37~{}M_{\odot}\,x_{h}^{1/12}\left(\frac{10^{-16}~{}\mathrm{eV}}{m_{a}}\right)^{3/2}\left(\frac{g_{*}(T_{\rm
osc})}{g_{*}(T_{0})}\right)^{1/12}~{},$ (3.10) $\displaystyle\tau H(z_{\rm
eq})$ $\displaystyle\approx$ $\displaystyle\frac{3.4\times
10^{-5}~{}x_{h}^{7/2}\left(\frac{g_{*}(T_{\rm
osc})}{g_{*}(T_{0})}\right)^{1/2}}{\log\left[2.2x_{h}^{11/12}\left(\frac{g_{*}(T_{\rm
osc})}{g_{*}(T_{0})}\right)^{1/6}\right]}~{}.$ (3.11)

Again, $T_{\rm osc}$ is determined by $m_{a}$ from Eq. 3.7. Here, we have
ignored the small difference between $g_{*}(T)$ and $g_{*,S}(T)$. With the
above expressions, we obtained all the equations needed to compute the axion
star mass function for ALPs. Note that $\tau$ is independent of $m_{a}$ except
for a minor dependence from $g_{*}$. This gives an almost universal mass
distribution shape for all $m_{a}$, although as $m_{a}$ is increased the peak
in the distribution of star masses also increases.

In Fig. 3, we present the same axion star properties as Fig. 1 for axion-like
particles. The black and red curves are the mass and density contours of ALP
stars, respectively. We also show the pivot axion minihalo mass $x_{h,p}$ with
the blue dashed line. Similar to the QCD axion scenario, axion minihalos with
smaller masses still have a larger axion star formation rate. Therefore
smaller $x_{h}$ can lead to axion star formation while large $x_{h}$
correspond to the grey region without any axion stars forming during matter-
radiation equality. Unlike QCD axion stars, the ALP star formation rate is
independent of $m_{a}$ except for the minor $g_{*}$ dependence. Thus, the red
curves and the boundary of the grey region are nearly horizontal. The mass
contours are almost vertical because the de Broglie wavelength of axion waves
sets the size of axion stars while the axion star densities are approximately
independent of $m_{a}$. A smaller $m_{a}$ corresponds to a larger $M_{\star}$,
which implies constraints on massive compact objects can place a lower bound
on the ALP mass.

As shown in Fig. 4, the shape of the mass function of ALP stars is invariant
for different $m_{a}$. The shape slightly changes between $m_{a}\lesssim
10^{-17}~{}\mathrm{eV}$ and $m_{a}\gtrsim 10^{-17}~{}\mathrm{eV}$ due to the
sudden change of $g_{*}$ from the electron-positron annihilations, but the
total mass fraction does not change much. Note that each mass function has a
divergent peak at $x_{h,p}$ as discussed in Sec. 3.1. As shown in Fig. 4, the
mass fraction of axions that have formed axion stars is $f_{\star}\sim 0.42$.
Such a large axion star fraction can lead to interesting observational
constraints such as dynamical heating effects in ultrafaint dwarfs [100, 101],
which we discuss in the next section.

The axion star mass function shown in Fig. 4 shows that the distribution of
axion star masses span roughly one order of magnitude, for a given axion mass.
We have thus far only considered the formation process around matter-radiation
equality. We expect that this dominates the overall rate because the axion
star formation rate is peaked around this time due to the large density of
axion minihalos formed at equality. At later times the halos, as the halo
masses grow, become more dilute, and the star formation time increases. We
have ignored the merger of axion stars at later times, which may increase the
average star mass and extend the star mass function at the high mass end.

##  4 Constraints

Figure 5: Current status of constraints on ultralight dark matter. Here we
compare our results with the existing constraints in the literature that apply
to various dark matter models. For generic fuzzy dark matter, there are
Lyman-$\alpha$ constraints ($m_{a}>2\times 10^{-21}$ eV [102]) and ultrafaint
dwarf galaxies constraints ($m_{a}>3\times 10^{-19}$ eV [103]). Lyman-$\alpha$
also constrains post-inflationary axions ($m_{a}>3\times 10^{-18}$ eV [104]),
which is the same scenario as discussed here, but the bound in this work is
stronger by one-order-of-magnitude ($m_{a}>3.3\times 10^{-17}\rm eV$). There
is also a generic bound on post-inflation dark matter ($m_{a}>10^{-19}$ eV)
considering the axion isocurvature perturbations and the free-streaming
effect, which will be again constrained by Lyman-$\alpha$ [105]. The
constraint from ALP stars is derived by considering their dynamical heating
effects in ultrafaint dwarfs [100, 101].

Using the axion star mass functions determined above, we now turn to possible
constraints on axion dark matter coming from the effects of these exotic
compact objects. As mentioned earlier, we focus only on axion stars
interacting through gravitational interactions, to be as model independent as
possible. For ALPs, we considered a wider, and generally smaller, range of
axion mass than for QCD axions which leads to axion stars which are more
massive than in the case of the QCD axion. We did consider axions of both
types in the mass range $10^{-8}~{}\mathrm{eV}<m_{a}<10^{-6}~{}\mathrm{eV}$,
and here although the range of axion star masses are comparable their
distribution and the total mass fraction are different.

Generally, over the QCD axion mass range the star masses are very light,
$10^{-14}M_{\odot}\lesssim M_{\star}\lesssim 10^{-11}M_{\odot}$. This means
that QCD axion stars are not massive or dense enough to lead to meaningful
constraints from galactic microlensing [106]. For the ALP scenario the star
mass is much heavier, for the range of axions we considered we find stars in
the range $0.1M_{\odot}\lesssim M_{\star}\lesssim 10^{8}M_{\odot}$. These
masses are beyond the reach of microlensing searches, due to the low number
density of the axion stars. However, this mass range can instead be
constrained by observations of ultra faint dwarfs, which we now discuss.

It has recently been pointed out that the motion of pointlike [107, 108, 109,
110, 100] and diffuse dark matter substructures [101] inside dwarf galaxies
can lead to heating of the stars in the dwarf, which in turn can place a
constraint on the number of these structures. The ensemble of dark matter
substructures is typically hotter than the gas of stars in these systems and
thus, through gravitational interactions, heat will flow from one system to
the other. This heating will increase the average speed, and thus the scale
radius, of the stars in the dwarf. From observations of the distribution of
stars in Segue-I, an ultra-faint dwarf (UFD) galaxy, Refs. [100, 101] placed
strong constraints on the fraction of dark matter that can be in these
substructures, as a function of the clump mass and density. The UFD
constraints were derived assuming that all clumps had the same mass and for
each mass they can be interpreted as a bound on the fraction of dark matter
that can be in substructures of that mass.

For a given axion model our analysis gives a prediction for the axion star
mass function, shown in Figs. 4. Since the contribution of each substructure
to heating UFD is linearly adding up, we can translate the constraints on a
single mass distribution to the constraints on a specific mass distribution
with

$\int
f_{\textrm{UFD}}^{-1}(M_{\star},\rho_{\star})\frac{df_{\star}}{dM_{\star}}dM_{\star}\leq
1\,,$ (4.1)

where $f_{\textrm{UFD}}$ is the constraint on the mass fraction of compact
objects assuming a single mass $M_{\star}$ and density $\rho_{\star}$ in Ref.
[101]. Any axion star mass function that does not satisfy Eq. 4.1 will be
ruled out. We can easily test if certain ALP masses are excluded by applying
the axion star mass function presented in Fig. 4 to the above integral Eq.
4.1. Lighter axion masses are more constrained because there are stronger
bounds for more massive objects. As a result, we exclude the ALP mass
$m_{a}<3.3\times 10^{-17}~{}\mathrm{eV}$. We present our result in Fig. 5 in
comparison to other constraints on fuzzy DM, post-inflation DM, and post-
inflation axions. Our work focuses on the observational constraints of post-
inflation axions, which feature an enhancement on the matter power spectrum as
well as the formation of axion minihalos and axion stars. The constraint on
axion minihalos only from Lyman-$\alpha$ has already been considered and leads
to a weaker constraint [104]. There are also generic bounds on the mass of
post-inflation dark matter [105] and fuzzy dark matter [102, 111, 103, 112],
which are also weaker than our bound obtained from the gravitational effects
of axion stars.

##  5 Discussion

In the post-inflationary axion scenario, where the PQ symmetry is broken after
inflation has ended, the axion field has a white-noise spectrum of
isocurvature perturbations at short scales. Once structure can grow, after
matter-radiation equality, these perturbations can quickly form axion
minihalos. The mass distribution of these minihalos can be calculated based on
the spectrum of perturbations. Furthermore, numerical simulations, and
analytic insights, indicate that axion stars, compact bound states of axions
held together by gravity, will form in the center of these minihalos.

Using the growth rate of the stars inside minihalos and the mass function of
the minihalos we determine the axion star mass function,
$df_{\star}/dM_{\star}$, for both QCD axion and axion-like particles. These
results are useful for understanding the rate of possible signals coming from
axion stars, such as the radio signals of axion stars, due to an axion-photon
coupling, from their collisions with other objects. We used the star mass
function to determine the signal coming from the axion’s gravitational
interactions _alone_. In particular, the effect of axion stars on the heating
of normal stars in ultra faint dwarfs.

On the dilute branch axion stars are stabilized by kinetic pressure and
gravity, the mass-radius relation is simply determined by the axion mass.
Therefore, the axion star mass function also contains information about the
density and radius of axion stars. We used this to translate recent bounds on
heating in ultra faint dwarfs, by dark compact objects of a single mass, into
bounds on axions. For QCD axion stars, we did not find any current
constraints. ALP stars, on the other hand, can be much more massive and this
approach places a constraint of $m_{a}>3.3\times 10^{-17}$ eV, see Fig. 5 for
a comparison of constraints. In addition, the ALP star mass function suggests
that almost half of the axion dark matter is in the form of ALP stars for a
broad range of ALP masses.

There are several ways in which our approach can be improved. First, our
calculation is based on a pure white-noise power spectrum of axion
isocurvature perturbations but axion strings will contribute to the power
spectrum at higher wavenumber, altering the spectrum. We did not look at
minihalos outside of the white noise regime but including them would lead to
more light minihalos. For the ALP case, this enhancement is not very relevant
because the lighter axion minihalos are not a significant fraction of the
total axion dark matter mass. For QCD axions with $m_{a}\gtrsim 10^{-6}$ eV,
this enhancement may be more important. Second, we focused on the axion star
population formed around the time of matter-radiation equality, which is
motivated by the fact that the axion star formation rate is on average peaked
at earlier times due to the high density of the average minihalos. We have not
considered late evolution such as the merger of axion stars, which leads to
more massive and thus more detectable axion stars. However, to include late
time evolution one has to be careful not to double count mass in axion stars,
a complete analysis would require careful tracking of stars in minihalos. We
note that tidal disruption is only a problem for minihalos close to galactic
centers and is unlikely to be an issue for denser axion stars.

#### Note added:

During the preparation of this work, Ref. [91] appeared, studying the
enhancements in QCD axion star formation based on the power spectrum from
string simulations, which extends to smaller scales. This allows the
consideration of larger QCD axion masses that can form axion stars. The range
of QCD axion masses considered is complementary to those considered here.

## Acknowledgement

We thank Peter Graham and Harikrishnan Ramani for the helpful discussions and
for providing us with the data associated with their UFD constraint. We also
thank Han Gil Choi, Wayne Hu, Xuheng Luo, and Dennis Maseizik for the useful
discussions. JHC, PJF, and HX are supported by Fermi Research Alliance, LLC
under Contract DE-AC02-07CH11359 with the U.S. Department of Energy.

## References

  * [1] R. D. Peccei and H. R. Quinn, _CP Conservation in the Presence of Instantons_ , _Phys. Rev. Lett._ 38 (1977) 1440–1443.
  * [2] F. Wilczek, _Problem of strong $p$ and $t$ invariance in the presence of instantons_, _Phys. Rev. Lett._ 40 (Jan, 1978) 279–282.
  * [3] J. E. Kim, _Weak-interaction singlet and strong $\mathrm{CP}$ invariance_, _Phys. Rev. Lett._ 43 (Jul, 1979) 103–107.
  * [4] L. Abbott and P. Sikivie, _A Cosmological Bound on the Invisible Axion_ , _Phys. Lett. B_ 120 (1983) 133–136.
  * [5] M. Dine and W. Fischler, _The Not So Harmless Axion_ , _Phys. Lett. B_ 120 (1983) 137–141.
  * [6] J. Preskill, M. B. Wise and F. Wilczek, _Cosmology of the Invisible Axion_ , _Phys. Lett._ B120 (1983) 127–132.
  * [7] R. D. Peccei, _The Strong CP problem and axions_ , _Lect. Notes Phys._ 741 (2008) 3–17, [hep-ph/0607268].
  * [8] ADMX collaboration, N. Du et al., _A Search for Invisible Axion Dark Matter with the Axion Dark Matter Experiment_ , _Phys. Rev. Lett._ 120 (2018) 151301, [1804.05750].
  * [9] ADMX collaboration, T. Braine et al., _Extended Search for the Invisible Axion with the Axion Dark Matter Experiment_ , _Phys. Rev. Lett._ 124 (2020) 101303, [1910.08638].
  * [10] D. Alesini et al., _Galactic axions search with a superconducting resonant cavity_ , _Phys. Rev. D_ 99 (2019) 101101, [1903.06547].
  * [11] S. Lee, S. Ahn, J. Choi, B. R. Ko and Y. K. Semertzidis, _Axion Dark Matter Search around 6.7 $\mu$eV_, _Phys. Rev. Lett._ 124 (2020) 101802, [2001.05102].
  * [12] D. Alesini et al., _Search for Galactic axions with a high-Q dielectric cavity_ , _Phys. Rev. D_ 106 (2022) 052007, [2208.12670].
  * [13] M. Buschmann, J. W. Foster, A. Hook, A. Peterson, D. E. Willcox, W. Zhang et al., _Dark matter from axion strings with adaptive mesh refinement_ , _Nature Commun._ 13 (2022) 1049, [2108.05368].
  * [14] M. Gorghetto, E. Hardy and G. Villadoro, _More Axions from Strings_ , _SciPost Phys._ 10 (2021) 050, [2007.04990].
  * [15] A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper and J. March-Russell, _String Axiverse_ , _Phys. Rev. D_ 81 (2010) 123530, [0905.4720].
  * [16] C. J. Hogan and M. J. Rees, _Axion Miniclusters_ , _Phys. Lett. B_ 205 (1988) 228–230.
  * [17] E. W. Kolb and I. I. Tkachev, _Axion miniclusters and Bose stars_ , _Phys. Rev. Lett._ 71 (1993) 3051–3054, [hep-ph/9303313].
  * [18] E. W. Kolb and I. I. Tkachev, _Femtolensing and picolensing by axion miniclusters_ , _Astrophys. J. Lett._ 460 (1996) L25–L28, [astro-ph/9510043].
  * [19] M. Fairbairn, D. J. E. Marsh, J. Quevillon and S. Rozier, _Structure formation and microlensing with axion miniclusters_ , _Phys. Rev. D_ 97 (2018) 083502, [1707.03310].
  * [20] A. Katz, J. Kopp, S. Sibiryakov and W. Xue, _Femtolensing by Dark Matter Revisited_ , _JCAP_ 12 (2018) 005, [1807.11495].
  * [21] J. A. Dror, H. Ramani, T. Trickle and K. M. Zurek, _Pulsar Timing Probes of Primordial Black Holes and Subhalos_ , _Phys. Rev. D_ 100 (2019) 023003, [1901.04490].
  * [22] H. Ramani, T. Trickle and K. M. Zurek, _Observability of Dark Matter Substructure with Pulsar Timing Correlations_ , 2005.03030.
  * [23] V. S. H. Lee, A. Mitridate, T. Trickle and K. M. Zurek, _Probing Small-Scale Power Spectra with Pulsar Timing Arrays_ , _JHEP_ 06 (2021) 028, [2012.09857].
  * [24] L. Dai and J. Miralda-Escudé, _Gravitational Lensing Signatures of Axion Dark Matter Minihalos in Highly Magnified Stars_ , _Astron. J._ 159 (2020) 49, [1908.01773].
  * [25] H. Xiao, L. Dai and M. McQuinn, _Detecting Dark Matter Substructures on Small Scales with Fast Radio Bursts_ , 2401.08862.
  * [26] B. Eggemeier, C. A. J. O’Hare, G. Pierobon, J. Redondo and Y. Y. Y. Wong, _Axion minivoids and implications for direct detection_ , _Phys. Rev. D_ 107 (2023) 083510, [2212.00560].
  * [27] C. A. J. O’Hare, G. Pierobon and J. Redondo, _Axion minicluster streams in the solar neighbourhood_ , 2311.17367.
  * [28] H. Xiao, I. Williams and M. McQuinn, _Simulations of axion minihalos_ , _Phys. Rev. D_ 104 (2021) 023515, [2101.04177].
  * [29] P. J. Fox, N. Weiner and H. Xiao, _Recurrent axion stars collapse with dark radiation emission and their cosmological constraints_ , _Phys. Rev. D_ 108 (2023) 095043, [2302.00685].
  * [30] P. J. Fox, N. Weiner and H. Xiao. To appear, 2024.
  * [31] E. Braaten, A. Mohapatra and H. Zhang, _Dense Axion Stars_ , _Phys. Rev. Lett._ 117 (2016) 121801, [1512.00108].
  * [32] L. Visinelli, S. Baum, J. Redondo, K. Freese and F. Wilczek, _Dilute and Dense Axion Stars_ , _Phys. Lett._ B777 (2018) 64–72, [1710.08910].
  * [33] D. G. Levkov, A. G. Panin and I. I. Tkachev, _Gravitational Bose-Einstein Condensation in the Kinetic Regime_ , _Phys. Rev. Lett._ 121 (2018) 151301, [1804.05857].
  * [34] E. Braaten, A. Mohapatra and H. Zhang, _Emission of Photons and Relativistic Axions from Axion Stars_ , _Phys. Rev. D_ 96 (2017) 031901, [1609.05182].
  * [35] D. G. Levkov, A. G. Panin and I. I. Tkachev, _Relativistic Axions from Collapsing Bose Stars_ , _Phys. Rev. Lett._ 118 (2017) 011301, [1609.03611].
  * [36] P.-H. Chavanis, _Collapse of a Self-Gravitating Bose-Einstein Condensate with Attractive Self-Interaction_ , _Phys. Rev. D_ 94 (2016) 083007, [1604.05904].
  * [37] X. Du, D. J. E. Marsh, M. Escudero, A. Benson, D. Blas, C. K. Pooni et al., _Soliton Merger Rates and Enhanced Axion Dark Matter Decay_ , 2301.09769.
  * [38] D. Maseizik and G. Sigl, _Distributions and Collision Rates of ALP Stars in the Milky Way_ , 2404.07908.
  * [39] H. Zhang, _Axion Stars_ , _Symmetry_ 12 (2019) 25, [1810.11473].
  * [40] J. Eby, M. Leembruggen, L. Street, P. Suranyi and L. C. R. Wijewardhana, _Global view of QCD axion stars_ , _Phys. Rev. D_ 100 (2019) 063002, [1905.00981].
  * [41] J. H.-H. Chan, S. Sibiryakov and W. Xue, _Boson star normal modes_ , _JHEP_ 08 (2023) 045, [2304.13054].
  * [42] G. Grilli di Cortona, E. Hardy, J. Pardo Vega and G. Villadoro, _The QCD axion, precisely_ , _JHEP_ 01 (2016) 034, [1511.02867].
  * [43] P.-H. Chavanis, _Mass-Radius Relation of Newtonian Self-Gravitating Bose-Einstein Condensates with Short-Range Interactions: I. Analytical Results_ , _Phys. Rev. D_ 84 (2011) 043531, [1103.2050].
  * [44] P. H. Chavanis and L. Delfini, _Mass-Radius Relation of Newtonian Self-Gravitating Bose-Einstein Condensates with Short-Range Interactions: II. Numerical Results_ , _Phys. Rev. D_ 84 (2011) 043532, [1103.2054].
  * [45] R. Ruffini and S. Bonazzola, _Systems of selfgravitating particles in general relativity and the concept of an equation of state_ , _Phys. Rev._ 187 (1969) 1767–1783.
  * [46] M. Membrado, J. Abad, A. F. Pacheco and J. Sanudo, _Newtonian Boson Spheres_ , _Phys. Rev. D_ 40 (1989) 2736–2738.
  * [47] J. Eby, P. Suranyi and L. C. R. Wijewardhana, _The Lifetime of Axion Stars_ , _Mod. Phys. Lett. A_ 31 (2016) 1650090, [1512.01709].
  * [48] J. Eby, M. Leembruggen, P. Suranyi and L. C. R. Wijewardhana, _Collapse of Axion Stars_ , _JHEP_ 12 (2016) 066, [1608.06911].
  * [49] J. Enander, A. Pargner and T. Schwetz, _Axion minicluster power spectrum and mass function_ , _JCAP_ 12 (2017) 038, [1708.04466].
  * [50] N. Blinov, M. J. Dolan and P. Draper, _Imprints of the Early Universe on Axion Dark Matter Substructure_ , _Phys. Rev. D_ 101 (2020) 035002, [1911.07853].
  * [51] D. Ellis, D. J. E. Marsh and C. Behrens, _Axion Miniclusters Made Easy_ , _Phys. Rev. D_ 103 (2021) 083525, [2006.08637].
  * [52] D. Ellis, D. J. E. Marsh, B. Eggemeier, J. Niemeyer, J. Redondo and K. Dolag, _Structure of axion miniclusters_ , _Phys. Rev. D_ 106 (2022) 103514, [2204.13187].
  * [53] K. M. Zurek, C. J. Hogan and T. R. Quinn, _Astrophysical Effects of Scalar Dark Matter Miniclusters_ , _Phys. Rev. D_ 75 (2007) 043511, [astro-ph/0607341].
  * [54] B. Eggemeier, J. Redondo, K. Dolag, J. C. Niemeyer and A. Vaquero, _First simulations of axion minicluster halos_ , 1911.09417.
  * [55] X. Shen, H. Xiao, P. F. Hopkins and K. M. Zurek, _Disruption of Dark Matter Minihaloes in the Milky Way environment: Implications for Axion Miniclusters and Early Matter Domination_ , 2207.11276.
  * [56] I. Tkachev, _Fast Radio Bursts and Axion Miniclusters_ , _JETP Lett._ 101 (2015) 1–6, [1411.3900].
  * [57] A. Iwazaki, _Axion stars and fast radio bursts_ , _Phys. Rev. D_ 91 (2015) 023008, [1410.4323].
  * [58] L. van Waerbeke and A. Zhitnitsky, _Fast Radio Bursts and the Axion Quark Nugget Dark Matter Model_ , _Phys. Rev. D_ 99 (2019) 043535, [1806.02352].
  * [59] T. Dietrich, F. Day, K. Clough, M. Coughlin and J. Niemeyer, _Neutron star–axion star collisions in the light of multimessenger astronomy_ , _Mon. Not. Roy. Astron. Soc._ 483 (2019) 908–914, [1808.04746].
  * [60] S. Sun and Y.-L. Zhang, _Gravitational Waves and Possible Fast Radio Bursts from Axion Clumps_ , 2003.10527.
  * [61] J. H. Buckley, P. B. Dev, F. Ferrer and F. P. Huang, _Fast radio bursts from axion stars moving through pulsar magnetospheres_ , 2004.06486.
  * [62] A. Caputo, M. Regis, M. Taoso and S. J. Witte, _Detecting the Stimulated Decay of Axions at RadioFrequencies_ , _JCAP_ 03 (2019) 027, [1811.08436].
  * [63] Y. Sun, K. Schutz, A. Nambrath, C. Leung and K. Masui, _Axion dark matter-induced echo of supernova remnants_ , _Phys. Rev. D_ 105 (2022) 063007, [2110.13920].
  * [64] P. S. B. Dev, F. Ferrer and T. Okawa, _On the Galactic radio signal from stimulated decay of axion dark matter_ , 2311.13653.
  * [65] Y. Sun, K. Schutz, H. Sewalls, C. Leung and K. W. Masui, _Looking in the axion mirror: An all-sky analysis of stimulated decay_ , 2310.03788.
  * [66] B. R. Safdi, Z. Sun and A. Y. Chen, _Detecting Axion Dark Matter with Radio Lines from Neutron Star Populations_ , _Phys. Rev. D_ 99 (2019) 123021, [1811.01020].
  * [67] R. A. Battye, B. Garbrecht, J. I. McDonald, F. Pace and S. Srinivasan, _Dark matter axion detection in the radio/mm-waveband_ , _Phys. Rev. D_ 102 (2020) 023504, [1910.11907].
  * [68] J. W. Foster, Y. Kahn, O. Macias, Z. Sun, R. P. Eatough, V. I. Kondratiev et al., _Green Bank and Effelsberg Radio Telescope Searches for Axion Dark Matter Conversion in Neutron Star Magnetospheres_ , _Phys. Rev. Lett._ 125 (2020) 171301, [2004.00011].
  * [69] A. Hook, Y. Kahn, B. R. Safdi and Z. Sun, _Radio Signals from Axion Dark Matter Conversion in Neutron Star Magnetospheres_ , _Phys. Rev. Lett._ 121 (2018) 241102, [1804.03145].
  * [70] Y. Bai and Y. Hamada, _Detecting Axion Stars with Radio Telescopes_ , _Phys. Lett. B_ 781 (2018) 187–194, [1709.10516].
  * [71] M. P. Hertzberg and E. D. Schiappacasse, _Dark Matter Axion Clump Resonance of Photons_ , _JCAP_ 11 (2018) 004, [1805.00430].
  * [72] M. P. Hertzberg, Y. Li and E. D. Schiappacasse, _Merger of Dark Matter Axion Clumps and Resonant Photon Emission_ , _JCAP_ 07 (2020) 067, [2005.02405].
  * [73] A. Prabhu, _Optical Lensing by Axion Stars: Observational Prospects with Radio Astrometry_ , 2006.10231.
  * [74] A. Prabhu and N. M. Rapidis, _Resonant Conversion of Dark Matter Oscillons in Pulsar Magnetospheres_ , 2005.03700.
  * [75] P. Carenza, A. Mirizzi and G. Sigl, _Dynamical evolution of axion condensates under stimulated decays into photons_ , _Phys. Rev. D_ 101 (2020) 103016, [1911.07838].
  * [76] K. Fujikura, M. P. Hertzberg, E. D. Schiappacasse and M. Yamaguchi, _Microlensing constraints on axion stars including finite lens and source size effects_ , _Phys. Rev. D_ 104 (2021) 123012, [2109.04283].
  * [77] M. Escudero, C. K. Pooni, M. Fairbairn, D. Blas, X. Du and D. J. E. Marsh, _Axion Star Explosions: A New Source for Axion Indirect Detection_ , 2302.10206.
  * [78] X. Gan, L.-T. Wang and H. Xiao, _Detecting Axion Dark Matter with Black Hole Polarimetry_ , 2311.02149.
  * [79] Z. Yin and L. Visinelli, _Axion star condensation around primordial black holes and microlensing limits_ , 2404.10340.
  * [80] A. Vaquero, J. Redondo and J. Stadler, _Early Seeds of Axion Miniclusters_ , _JCAP_ 04 (2019) 012, [1809.09241].
  * [81] M. Buschmann, J. W. Foster and B. R. Safdi, _Early-Universe Simulations of the Cosmological Axion_ , _Phys. Rev. Lett._ 124 (2020) 161103, [1906.00967].
  * [82] C. A. J. O’Hare, G. Pierobon, J. Redondo and Y. Y. Y. Wong, _Simulations of axionlike particles in the postinflationary scenario_ , _Phys. Rev. D_ 105 (2022) 055025, [2112.05117].
  * [83] W. H. Press and P. Schechter, _Formation of Galaxies and Clusters of Galaxies by Self-Similar Gravitational Condensation_ , _ApJ_ 187 (Feb., 1974) 425–438.
  * [84] J. F. Navarro, C. S. Frenk and S. D. M. White, _The structure of cold dark matter halos_ , _The Astrophysical Journal_ 462 (May, 1996) 563.
  * [85] B. Eggemeier and J. C. Niemeyer, _Formation and Mass Growth of Axion Stars in Axion Miniclusters_ , _Phys. Rev. D_ 100 (2019) 063528, [1906.01348].
  * [86] J. Chen, X. Du, E. Lentz and D. J. E. Marsh, _Relaxation times for Bose-Einstein condensation by self-interaction and gravity_ , 2109.11474.
  * [87] K. Kirkpatrick, A. E. Mirasola and C. Prescod-Weinstein, _Relaxation Times for Bose-Einstein Condensation in Axion Miniclusters_ , _Phys. Rev. D_ 102 (2020) 103012, [2007.07438].
  * [88] J. Chen, X. Du, E. W. Lentz, D. J. E. Marsh and J. C. Niemeyer, _New Insights into the Formation and Growth of Boson Stars in Dark Matter Halos_ , _Phys. Rev. D_ 104 (2021) 083022, [2011.01333].
  * [89] J. H.-H. Chan, S. Sibiryakov and W. Xue, _Condensation and Evaporation of Boson Stars_ , 2207.04057.
  * [90] A. S. Dmitriev, D. G. Levkov, A. G. Panin and I. I. Tkachev, _Self-similar growth of Bose stars_ , 2305.01005.
  * [91] M. Gorghetto, E. Hardy and G. Villadoro, _More Axion Stars from Strings_ , 2405.19389.
  * [92] M. P. Hertzberg, M. Tegmark and F. Wilczek, _Axion Cosmology and the Energy Scale of Inflation_ , _Phys. Rev. D_ 78 (2008) 083507, [0807.1726].
  * [93] M. Gorghetto, E. Hardy and G. Villadoro, _Axions from Strings: the Attractive Solution_ , _JHEP_ 07 (2018) 151, [1806.04677].
  * [94] M. Kawasaki, T. Sekiguchi, M. Yamaguchi and J. Yokoyama, _Long-term dynamics of cosmological axion strings_ , _PTEP_ 2018 (2018) 091E01, [1806.05566].
  * [95] K. Saikawa, J. Redondo, A. Vaquero and M. Kaltschmidt, _Spectrum of global string networks and the axion dark matter mass_ , 2401.17253.
  * [96] H. Kim, J. Park and M. Son, _Axion Dark Matter from Cosmic String Network_ , 2402.00741.
  * [97] S. Borsanyi, Z. Fodor, K. H. Kampert, S. D. Katz, T. Kawanai, T. G. Kovacs et al., _Lattice qcd for cosmology_ , 2016.
  * [98] N. Blinov, M. J. Dolan, P. Draper and J. Kozaczuk, _Dark matter targets for axionlike particle searches_ , _Phys. Rev. D_ 100 (2019) 015049, [1905.06952].
  * [99] A. E. Nelson and H. Xiao, _Axion Cosmology with Early Matter Domination_ , _Phys. Rev. D_ 98 (2018) 063516, [1807.07176].
  * [100] P. W. Graham and H. Ramani, _Constraints on Dark Matter from Dynamical Heating of Stars in Ultrafaint Dwarfs. Part 1: MACHOs and Primordial Black Holes_ , 2311.07654.
  * [101] P. W. Graham and H. Ramani, _Constraints on Dark Matter from Dynamical Heating of Stars in Ultrafaint Dwarfs. Part 2: Substructure and the Primordial Power Spectrum_ , 2404.01378.
  * [102] V. Iršič, M. Viel, M. G. Haehnelt, J. S. Bolton and G. D. Becker, _First constraints on fuzzy dark matter from Lyman- $\alpha$ forest data and hydrodynamical simulations_, _Phys. Rev. Lett._ 119 (2017) 031302, [1703.04683].
  * [103] N. Dalal and A. Kravtsov, _Excluding fuzzy dark matter with sizes and stellar kinematics of ultrafaint dwarf galaxies_ , _Phys. Rev. D_ 106 (2022) 063517, [2203.05750].
  * [104] V. Iršič, H. Xiao and M. McQuinn, _Early structure formation constraints on the ultralight axion in the postinflation scenario_ , _Phys. Rev. D_ 101 (2020) 123518, [1911.11150].
  * [105] M. A. Amin and M. Mirbabayi, _A Lower Bound on Dark Matter Mass_ , _Phys. Rev. Lett._ 132 (2024) 221004, [2211.09775].
  * [106] D. Croon, D. McKeen, N. Raj and Z. Wang, _Subaru-HSC through a different lens: Microlensing by extended dark matter structures_ , _Phys. Rev. D_ 102 (2020) 083021, [2007.12697].
  * [107] Q. Zhu, E. Vasiliev, Y. Li and Y. Jing, _Primordial black holes as dark matter: constraints from compact ultra-faint dwarfs_ , _Monthly Notices of the Royal Astronomical Society_ 476 (01, 2018) 2–11, [https://academic.oup.com/mnras/article-pdf/476/1/2/24196942/sty079.pdf].
  * [108] S. L. Zoutendijk, J. Brinchmann, L. A. Boogaard, M. L. Gunawardhana, T.-O. Husser, S. Kamann et al., _The MUSE-Faint survey: I. Spectroscopic evidence for a star cluster in Eridanus 2 and constraints on MACHOs as a constituent of dark matter_ , _Astronomy and Astrophysics - A &A_ 635 (2020) A107.
  * [109] J. Stegmann, P. R. Capelo, E. Bortolas and L. Mayer, _Improved constraints from ultra-faint dwarf galaxies on primordial black holes as dark matter_ , _Mon. Not. Roy. Astron. Soc._ 492 (2020) 5247–5260, [1910.04793].
  * [110] T. D. Brandt, _Constraints on MACHO Dark Matter from Compact Stellar Systems in Ultra-Faint Dwarf Galaxies_ , _Astrophys. J. Lett._ 824 (2016) L31, [1605.03665].
  * [111] K. Schutz, _Subhalo mass function and ultralight bosonic dark matter_ , _Phys. Rev. D_ 101 (2020) 123026, [2001.05503].
  * [112] T. Zimmermann, J. Alvey, D. J. E. Marsh, M. Fairbairn and J. I. Read, _Dwarf galaxies imply dark matter is heavier than $\mathbf{2.2\times 10^{-21}}\,\mathbf{eV}$_, 2405.20374.