# ALP dark matter with non-periodic potentials: parametric resonance, halo
formation and gravitational signatures

Aleksandr Chatrchyan  Cem Eröncel  Matthias Koschnitzke  Géraldine Servant

###### Abstract

Axion-like particles (ALPs) are leading candidates to explain the dark matter
in the universe. Their production via the misalignment mechanism has been
extensively studied for cosine potentials characteristic of pseudo-Nambu-
Goldstone bosons. In this work we investigate ALPs with non-periodic
potentials, which allow for large misalignment of the field from the minimum.
As a result, the ALP can match the relic density of dark matter in a large
part of the parameter space. Such potentials give rise to self-interactions
which can trigger an exponential growth of fluctuations in the ALP field via
parametric resonance, leading to the fragmentation of the field. We study
these effects with both Floquet analysis and lattice simulations. Using the
Press-Schechter formalism, we predict the halo mass function and halo spectrum
arising from ALP dark matter. These halos can be dense enough to produce
observable gravitational effects such as astrometric lensing, diffraction of
gravitational wave signals from black hole mergers, photometric microlensing
of highly magnified stars, perturbations of stars in the galactic disk or
stellar streams. These effects would provide a probe of dark matter even if it
does not couple to the Standard Model. They would not be observable for halos
predicted for standard cold dark matter and for ALP dark matter in the
standard misalignment mechanism. We determine the relevant regions of
parameter space in the (ALP mass, decay constant)-plane and compare
predictions in different axion fragmentation models.

DESY-23-060

##  1 Introduction

One prominent idea to explain dark matter is to introduce a light scalar
particle called axion [1, 2], which arises as the pseudo-Nambu-Goldstone boson
[3, 4, 5] of a new spontaneously broken global $\mathrm{U}(1)$-symmetry. The
latter was first introduced as extension of the Standard Model by Peccei and
Quinn [6] not to solve the puzzle of dark matter, but the strong CP-problem in
QCD. In this work, we will focus on a more general class of particles, usually
referred to as axion-like particles (ALPs), which do not generally solve the
strong CP-problem, but are still highly motivated, as they arise naturally in
many Standard-Model extensions, in particular in string theory [7, 8, 9].

While interactions between axions and other particles are assumed to be too
weak for axions to get into thermal equilibrium with the rest of the early-
universe plasma, the most promising cosmological mechanism to produce axions
or ALPs is the vacuum misalignment (or vacuum realignment) mechanism in which
the axion or the ALP field, modelled as classical scalar field due to its
bosonic nature and high occupation numbers, has a non-zero initial field value
and non-zero potential energy in the early universe, leading to oscillations
of the field that let it act as (dark) matter component.

In most of the models, axions and ALPs have a periodic potential, which arises
due to non-perturbative instanton effects [10]. In this work we focus on an
alternative type of potentials, which are non-periodic and allow larger
displacements of the ALP field from the minimum. Such potentials can be
generated from interactions with strongly coupled Yang-Mills gauge fields [11,
12, 13, 14]. Non-periodic potentials are also motivated by axion-monodromy
[15, 16, 17]. Oscillations in such potentials can lead to a strong growth of
fluctuations of the field due to parametric resonance [18, 19]. This process
can lead to fragmentation of the field in which all of the energy density is
transferred to the fluctuations [20].

Recently, there has been a lot of interest in the study of parametric
resonance in ALP models. Refs. [21, 22, 23] show that parametric resonance can
be effective also for a periodic potential if the initial angle is very close
to the top of the potential. Furthermore, Refs. [20, 24] demonstrate that the
parametric resonance is very efficient if the ALP field has a large initial
kinetic energy as in the Kinetic Misalignment Mechanism [25, 26]. Similar
effects have also been observed in the axiverse models where two ALP fields
have similar masses [27, 28], when the discrete symmetry of the ALPs is broken
by a quadratic monomial [29, 30, 31], when the ALP has an $\alpha$-attractor-
type potential [32], and also in models where the ALP potential is
temperature-dependent, such as the QCD axion [33, 34]. Very recently, Ref.
[35] claims that oscillations of the Hubble parameter induced by the ALP
oscillations can also cause fragmentation of the low-wavelength ALP modes.

The growth of fluctuations has important consequences. Once gravitational
interactions become important, the overdense regions, corresponding to the
fluctuations, collapse to halos. This process can be described by the Press-
Schechter formalism [36] and its modifications. For the amplified scales, the
halos can differ significantly from the generic prediction for cold dark
matter (CDM) halos at the respective size and could lead to observable effects
purely caused by gravitational interactions, similar to what was found in ref.
[22] for an ALP field with a cosine potential, and in ref. [37] for an ALP
field with initial kinetic energy.

The outline of this work is the following: The next section comprises a
detailed investigation of the vacuum misalignment mechanism, where we focus on
the evolution of fluctuations of the ALP field during this process. Using the
semi-analytic Floquet analysis we explain why parametric resonance is
inefficient in the case of the cosine potential, which serves as a motivation
for our studies of non-periodic potentials. In section 3 we compute the energy
density power spectrum of the ALP field, resulting from parametric resonance
and fragmentation, using both linear evolution of the modes, as well as fully
non-linear lattice simulations. In section 4 we take the results of the
preceding chapter and examine how the fluctuations decouple from the Hubble
flow and evolve to virialised halos after gravitational interaction becomes
important. The main result, to be discussed in section 5, is the
identification of the parts of the parameter space of ALP mass $m_{a}$ and
symmetry breaking scale $f_{a}$ in which we predict halos that have observable
gravitational effects, assuming that the ALPs make up all of dark matter.
Finally, we compare this scenario with others that lead to fragmentation of
the ALP field in section 6. There, we estimate that all the scenarios we have
considered predict dense halos in a roughly similar region in the
$(m_{a},f_{a})$-parameter space which we show in figure 1 with a blue band
labelled as "Dense Halo Region". In this plot, we also show the experimental
constraints and projections assuming that the ALP field has a KSVZ-like
coupling [38, 39] to the photon and the neutron. For simplicity, we consider
ALPs with constant mass. The different predictions for axion fragmentation due
to a temperature-dependent ALP mass were discussed in ref. [24].

Throughout this work, unless otherwise noted, we work in units where
$c=\hbar=1$, use the mostly minus convention for $g_{\mu\nu}$, the reduced
Planck mass $M_{\mathrm{pl}}=\sqrt{1/(8\pi G)}$ and
$H_{0}=h100\,\mathrm{km}\,\mathrm{s}^{-1}\,\mathrm{Mpc}^{-1}$ with $h=0.68$.

Figure 1: A sketch of the region in the ALP parameter space where dense halos
are expected to form, together with all the experimental constraints and
projections on ALPs assuming a coupling to the electromagnetic field with eq.
(6.8). We have obtained this region by combining the regions where dense halos
are expected from the Kinetic Misalignment Mechanism and from ALPs with non-
periodic potentials considered in this work. For caveats about this plot, see
section 6. All the data for the constraints and projections are compiled from
ref. [40].

##  2 From misalignment mechanism to parametric resonance

In this section, we investigate the significance of parametric resonance for
ALP dark matter produced from the vacuum misalignment mechanism. Both
periodic, as well as non-periodic ALP potentials are considered.

After a brief review of the essentials of the mechanism, as well as
estimations of the relic density in section 2.1, we focus on the dynamics of
ALP fluctuations in the subsequent subsections. The corresponding equations of
motion are presented in section 2.2. We then discuss the growth of
fluctuations in a periodic, “cosine” potential, and in a non-periodic, “power-
law” potential in sections 2.3 and 2.4, respectively. We demonstrate why in
the first case Hubble friction prevents the fluctuation modes from
experiencing an amplification unless the misalignment field value is very
close to the top of the potential. Then, we show how this constraint is
avoided in the non-periodic case.

###  2.1 ALP dark matter from the misalignment mechanism

Ignoring interactions with all other particles and assuming a minimal coupling
to gravity, the action of an ALP field is given by

$S_{\phi}=\int d^{4}x\sqrt{-g}\left[\frac{1}{2}\partial_{\mu}\phi
g^{\mu\nu}\partial_{\nu}\phi-V(\phi)\right]\,,$ (2.1)

where $g$ is the determinant of the FLRW metric $g_{\mu\nu}$.

We assume that the ALPs were present during inflation. In this case, similar
to the conventional misalignment mechanism for the QCD axion, inflation leaves
the ALP field with an almost homogeneous value inside our Hubble patch, which
is essentially frozen before $H\sim m_{a}$. We denote this initial
misalignment value by $\phi_{i}$. Throughout this work, we assume
$m_{a}=\,$cst and $\dot{\phi}_{i}=0$. The field obeys the classical equation
of motion (EOM):

$\ddot{\phi}+3H\dot{\phi}-\frac{1}{a^{2}}\nabla^{2}\phi+V^{\prime}(\phi)=0\,,$
(2.2)

where we denoted $\partial V/\partial\phi$ as $V^{\prime}(\phi)$. In the
radiation-dominated era, once the Hubble friction term becomes subdominant,
the field rolls to the minimum of the potential and oscillates around it at
later times. Near the minimum of the potential, where
$V\approx\frac{1}{2}m_{a}^{2}\phi^{2}$, the equation of state of such an
oscillating ALP field averages to $w=0$ so that it behaves as a (dark) matter
component.

Assuming that oscillation starts in the radiation era, and $m_{a}=\,$cst, a
general expression for the relic energy density of ALPs today is given by

$\boxed{\Omega_{a,0}=\frac{1}{3}(\Omega_{r,0})^{3/4}\frac{g_{s}(T_{0})}{g_{s}(T_{\mathrm{osc}})}\left(\frac{g_{\rho}(T_{\mathrm{osc}})}{g_{\rho}(T_{0})}\right)^{3/4}\left(\frac{1}{H_{\mathrm{osc}}}\right)^{3/2}\left(\frac{1}{H_{0}}\right)^{1/2}\left(\frac{1}{M_{\mathrm{pl}}}\right)^{2}V(\phi_{i})\mathcal{Z},}$
(2.3)

where $H_{\mathrm{osc}}$ is the Hubble scale at the onset of oscillation,
$g_{s}$ and $g_{\rho}$ are the effective degrees of freedom in entropy and
energy, respectively, while $\Omega_{r,0}$ is today’s density parameter of
radiation. For this estimate one uses that the energy density scales
approximately as
$\rho_{a}=\rho_{a,\,\mathrm{osc}}\left({a_{\mathrm{osc}}}/{a}\right)^{3}$ for
$a>a_{\mathrm{osc}}$, where $\rho_{a,\,\mathrm{osc}}=V(\phi_{i})$ is the
energy density of the field before the onset of oscillation. $\mathcal{Z}$
incorporates corrections to this estimate. For a harmonic potential, or close
to the minimum of a general potential, $\mathcal{Z}$ can be found to be

$\mathcal{Z}=\frac{8}{\pi}\quantity(\Gamma(5/4))^{2}\simeq 2.1,$ (2.4)

when setting $H_{\mathrm{osc}}=m_{a}$. Eq. (2.3) can differ strongly from the
harmonic estimate due to anharmonic effects, as it will be discussed in the
next subsections.

####  2.1.1 Relic density for a periodic potential

In this section we consider the usual periodic axion potential

$V(\phi)=m_{a}^{2}f_{a}^{2}[1-\cos(\phi/f_{a})]\,,$ (2.5)

where $f_{a}$ is the symmetry breaking scale, also referred to as the decay
constant. For small initial field values, the evolution of the field in this
potential defines the well-studied ’standard misalignment mechanism’. For
small initial field values, $\phi_{i}/f_{a}\lesssim 1$, the potential looks
almost harmonic and the relic density can be found from eq. (2.3) and eq.
(2.4). Since the axion mass $m_{a}$ and the initial field value determine when
the axion field begins to oscillate and decay, in a scenario with small
initial field value, the axion can only match the relic density for large
decay constants $f_{a}$. To explain dark matter from the axion-like-particle
in a larger part of the parameter space, i.e. also in the region of small
$f_{a}$, where dark matter from an ALP with small initial field value would be
underproduced, the onset of oscillations of the homogeneous mode (or zero
mode) has to be delayed, such that more of the maximal energy budget of
$m_{a}^{2}f_{a}^{2}$ can be used.

In a model with the cosine potential, the onset of oscillations is delayed if
the axion field starts close to top of the potential. This mechanism is called
Large Misalignment Mechanism or extreme axion and was described in ref. [22]
and [23]. Examples on how such an initial condition could be realised can be
found in refs. [41, 42]. Defining

$\theta=\phi/f_{a},$

Ref. [22] found as empirical expression for the onset of oscillation for
values close to the top with $\left|\pi-\theta_{i}\right|<10^{-2}$,

$t_{\mathrm{osc}}=\frac{1}{m_{a}}\ln\left(\frac{1}{\pi-\left|\theta_{i}\right|}\frac{2^{1/4}\pi^{1/2}}{\Gamma(5/4)}\right),$
(2.6)

which we confirmed numerically for tunings up to
$\absolutevalue{\pi-\theta_{i}}=10^{-12}$. Using
$H_{\mathrm{osc}}=1/(2t_{\mathrm{osc}})$, we can plug this into eq. (2.3) and
find for the relic density:

$\boxed{\begin{split}\Omega_{a,0}=&~{}\frac{2^{3/2}}{3}(\Omega_{r,0})^{3/4}\frac{g_{s}(T_{0})}{g_{s}(T_{\mathrm{osc}})}\left(\frac{g_{\rho}(T_{\mathrm{osc}})}{g_{\rho}(T_{0})}\right)^{3/4}\left(\frac{m_{a}}{H_{0}}\right)^{1/2}\left(\frac{f_{a}}{M_{\mathrm{pl}}}\right)^{2}\\\
&\times\left[\ln\left(\frac{1}{\pi-\left|\theta_{i}\right|}\frac{2^{1/4}\pi^{1/2}}{\Gamma(5/4)}\right)\right]^{3/2}\left[1-\cos(\theta_{i})\right].\end{split}}$
(2.7)

We show the parameter region where an axion-like particle with a cosine
potential, initially at rest, can make up all of the DM,
$\Omega_{a,0}h^{2}=\Omega_{\mathrm{DM,0}}h^{2}=0.12$, in figure 2.

Figure 2: The required initial field value $\phi_{i}/f_{a}$ to match the dark
matter relic density at each point in the $(m_{a},f_{a})$-plane for ALPs with
a cosine potential. The relic density is given by eq. (2.3) with the cosine
potential and $\mathcal{Z}=2.1$ for small initial field values and eq. (2.7)
for initial field values close to the top.

####  2.1.2 Relic density for a non-periodic potential

Following refs. [19, 14] we take as a generic non-periodic potential for the
ALP field:

$\boxed{V(\phi)=\frac{m_{a}^{2}f_{a}^{2}}{2p}\left[\left(1+\frac{\phi^{2}}{f_{a}^{2}}\right)^{p}-1\right].}$
(2.8)

As in the previous case, $m_{a}$ is the mass near the minimum, while the decay
constant $f_{a}$ is supposed to make a link to the UV completion of the ALP.
The parameter $p$ can be adjusted to incorporate different models. For $p=1$
we recover the generic quadratic potential. The values $0<p<1$ lead to axion-
monodromy potentials [15, 16, 17], while the values $p<0$ represent potentials
with a plateau at large field values [32]. The latter can arise when the axion
is coupled to pure strongly coupled Yang-Mills gauge fields, which has been
shown by considering a $SU(N)$ gauge theory in the large $N$ limit [11, 12,
13, 14]. For more details on possible UV-completions see appendix A. A
comparison of the potentials for different $p$-values is shown in figure 3.

Figure 3: Comparison of the cosine potential (blue solid line) with the non-
periodic potential eq. (2.8) for different $p$ values. $p=1$ corresponds to a
harmonic potential (black dotted line), while the potential is flatter for
$p=1/2$ (orange dashed line) and for negative values for $p$, such as $p=-1/2$
(green dashed-dotted line), approaches a plateau for large values of
$\phi/f_{a}$.

A good estimate for the beginning of oscillations is given by [32]

$H_{\mathrm{osc}}\cong\sqrt{\frac{V^{\prime}(\phi_{i})}{\phi_{i}}}\,.$ (2.9)

The right hand side reduces to $m_{a}$ in the case of a quadratic potential,
but is smaller for $p<1$, meaning that the onset of oscillations is delayed in
this case.

Considering axion masses such that oscillation starts in radiation era, hence
starting from eq. (2.3) and using eq. (2.9), we find:

$\boxed{\Omega_{a,0}=\frac{1}{6p}(\Omega_{r,0})^{3/4}\frac{g_{s}(T_{0})}{g_{s}(T_{\mathrm{osc}})}\left(\frac{g_{\rho}(T_{\mathrm{osc}})}{g_{\rho}(T_{0})}\right)^{3/4}\left(\frac{m_{a}}{H_{0}}\right)^{1/2}\left(\frac{f_{a}}{M_{\mathrm{pl}}}\right)^{2}\frac{(1+\theta_{i}^{2})^{p}-1}{(1+\theta_{i}^{2})^{3(p-1)/4}}\mathcal{Z}_{p}(\theta).}$
(2.10)

Here $\mathcal{Z}_{p}(\theta)$ contains the corrections due to the
anharmonicity of the potential, which leads to a delayed relaxation of the
average equation of state to $w=0$. As a fit function to numerical results we
used
$\mathcal{Z}_{p}(\theta)=\mathcal{Z}_{1}[1+\beta_{p}((1+\theta^{2})^{\alpha_{p}}-1)]$
where we found $\alpha_{1/2}=0.25$, $\beta_{1/2}=1.33/\mathcal{Z}_{1}$,
$\alpha_{-1/2}=0.26$, $\beta_{-1/2}=4.18/\mathcal{Z}_{1}$, and
$\mathcal{Z}_{1}$ is given by eq. (2.4). The delayed onset of oscillations is
accounted for by using eq. (2.9).

For $p\neq 1$ and large $\theta_{i}$, this can significantly change the
$(m_{a},f_{a})$-parameter space for ALP DM with respect to the case of
periodic potentials. We find the initial field values to match the ALP relic
density to the one of dark matter in the $(m_{a},f_{a})$-plane and plot the
result for different values of $p$ in figure 4. In contrast to the cosine
potential, large initial field values make it possible to cover also the
parameter space of smaller $m_{a}$ and $f_{a}$.

Figure 4: The required initial field value $\phi_{i}/f_{a}$ to match the dark
matter relic density at each point in the $(m_{a},f_{a})$-plane according to
eq. (2.10) for ALPs with the non-periodic potential. We set $p=-1/2$ and
$p=1/2$ in the left and the right panels, respectively. Compare with figure 2,
which shows the same parameter space for ALPs with the cosine potential.

As it will be discussed in the next sections, when the field undergoes
fragmentation, a significant fraction of its energy is transferred from the
homogeneous mode into initially relativistic fluctuations. This effects leads
to an additional correction to the prefactor $\mathcal{Z}$ since the equation
of state parameter shifts for some time towards $w=1/3$, as one would expect
for a gas of relativistic particles. As we clarify in section 5.2, in the
interesting part of the parameter region this correction is of
$\mathcal{O}(1)$ [30, 24]. Therefore, for simplicity, we do not include this
correction in our calculations.

###  2.2 Fluctuations of the field

It is convenient to separate the field into its homogeneous mode and
fluctuations,

$\hat{\phi}(t,\mathbf{x})=\phi(t)+\hat{\delta\phi}(t,\mathbf{x})\,.$ (2.11)

The fluctuation field can be expressed in terms of time dependent mode
functions $u_{\vec{k}}(t)$ for the individual Fourier modes [20]:

$\hat{\delta\phi}(t,\vec{x})=\int\frac{d^{3}k}{(2\pi)^{3}}\left[\hat{a}_{\vec{k}}u_{\vec{k}}(t)e^{i\vec{k}\cdot\vec{x}}+\hat{a}_{\vec{k}}^{\dagger}u_{\vec{k}}^{*}(t)e^{-i\vec{k}\cdot\vec{x}}\right]\,,$
(2.12)

where the creation and annihilation operators $\hat{a}_{\vec{k}}$ and
$\hat{a}_{\vec{k}}^{\dagger}$ obey the usual commutation relations of scalar
fields [20]. Since absolute direction of the momenta is not important when
assuming homogeneity and isotropy of the ALP field, we are only interested in
absolute values of $\vec{k}$ and will from now on mostly use
$k\equiv\lvert{\vec{k}}\rvert$, especially for indices.

Due to the small mass of the considered ALPs, the fluctuations are usually
characterised by large occupation numbers. As a consequence, these
fluctuations can be treated as classical(-statistical) ones. Having this in
mind, in the following, we will drop the hats on $\delta\phi$.

To separate equations of motions for the homogeneous mode and the small
fluctuations, $\phi(t)\gg\delta\phi(t,\vec{x})$, we expand the potential with
respect to the fluctuations and write

$V(\phi+\delta\phi)=V(\phi)+V^{\prime}(\phi)\delta\phi+\frac{1}{2}V^{\prime\prime}(\phi)\delta\phi^{2}+\mathcal{O}(\delta\phi^{3})\,.$
(2.13)

Keeping only terms up to $\mathcal{O}(\delta\phi)$ in the potential and using
eq. (2.12), the result is an equation that decomposes into an integral part
including fluctuations and a non-integral part only containing the homogeneous
mode. This leads to the general EOMs for the homogeneous mode as well as the
mode functions,

$\displaystyle\ddot{\phi}(t)+3H\dot{\phi}(t)+V^{\prime}(\phi)=0\,,$ (2.14)
$\displaystyle\ddot{u}_{k}(t)+3H\dot{u}_{k}(t)+\left(\frac{\vec{k}^{2}}{a^{2}}+V^{\prime\prime}(\phi)\right)u_{k}(t)=0\,.$
(2.15)

Similarly, the total energy density separates into a part depending on the
homogeneous mode and one part that depends on the fluctuations,

$\langle{\rho}\rangle=\frac{1}{2}\dot{\phi}(t)^{2}+V(\phi)\\\
+\frac{1}{2}\int\frac{d^{3}k}{(2\pi)^{3}}\left[\lvert\dot{u}_{k}(t)\rvert^{2}+\left(\frac{\vec{k}^{2}}{a^{2}}+V^{\prime\prime}(\phi)\right)\lvert
u_{k}(t)\rvert^{2}\right]\,,$ (2.16)

The evolution of the mode functions is trivial in the case of a quadratic
potential, where $V^{\prime\prime}(\phi)$ is simply the mass squared. On the
other hand, the anharmonicities of the potential can introduce instabilities
for the fluctuations, which are discussed next.

###  2.3 Growth of fluctuations in a periodic potential

In this section we want to find out when the exponential growth of
fluctuations is possible for the case of the cosine potential from eq. (2.5).

To be able to work with the EOMs numerically and to make the equations
independent of the axion mass, we introduce a dimensionless time- and a
dimensionless momentum-variable (as in ref. [22, 24]):

$\tilde{t}=m_{a}t\,,\quad\quad\quad\tilde{k}^{2}=\frac{k^{2}}{2m_{a}a^{2}H}\,.$
(2.17)

From now on, the dot above a function will denote derivation with respect to
the time variable given in the argument. During radiation era we have $a\sim
t^{1/2}$ and $H=\frac{1}{2t}$, therefore $a^{2}H$ stays constant during this
epoch. From these definitions we can see that $\tilde{k}$ corresponds to the
ratio of the physical momentum to the axion mass at $\tilde{t}=1$, i.e. how
relativistic this mode is at this time, which is roughly when the homogeneous
mode begins to oscillate.

We can reformulate the equations of motion in terms of the dimensionless
quantities, valid during radiation era:

$\displaystyle\ddot{\theta}(\tilde{t})+\frac{3}{2\tilde{t}}\dot{\theta}(\tilde{t})+\sin(\theta(\tilde{t}))=0\,,$
(2.18)
$\displaystyle\ddot{u}_{\tilde{k}}(\tilde{t})+\frac{3}{2\tilde{t}}\dot{u}_{\tilde{k}}(\tilde{t})+\left(\frac{\tilde{k}^{2}}{\tilde{t}}+\cos(\theta(\tilde{t}))\right)u_{\tilde{k}}(\tilde{t})=0\,.$
(2.19)

Here we have used that $H=m_{a}/(2\tilde{t})$ and
$k^{2}/a^{2}=\tilde{k}^{2}m_{a}^{2}/\tilde{t}$ in radiation era.

####  2.3.1 Floquet analysis

Let us ignore expansion for a moment, i.e. set $H=0$ and $a=1$. This implies
that the momenta do not redshift and that the corresponding term in the EOM is
$k^{2}$ instead of $k^{2}/a^{2}$. One then arrives at the following equations
for the homogeneous mode and the fluctuations:

$\displaystyle\ddot{\theta}(t)+m^{2}\sin(\theta)=0\,,$ (2.20)
$\displaystyle\ddot{u}_{k}(t)+\left[k^{2}+m^{2}\cos(\theta)\right]u_{k}(t)=0\,.$
(2.21)

Since $\theta$ is periodic, this equation is a second-order differential
equation with periodic coefficients, so it has the form of the Hill’s equation
[43, 44]. According to the Floquet theorem [45], the solutions are of the form

$u_{k}(t)=u_{+}(k;t)e^{\mu_{k}t}+u_{-}(k;t)e^{-\mu_{k}t},$ (2.22)

where $u_{\pm}$’s are periodic functions in time, and $\mu_{k}$’s are in
general complex coefficients. A necessary, but not sufficient, condition for
parametric resonance is that the Floquet exponent has non-zero real part
$\Re{\mu_{k}}>0$. So determining whether parametric resonance occurs amounts
to calculating the Floquet exponents. For more details on Floquet theory see
appendix B.

Figure 5: The contour of the real part of the Floquet exponents for the cosine
potential calculated via eq. (2.23). Larger values of
$\mathfrak{R}[\mu_{k}]/m$ (warmer colours) indicate stronger parametric
resonance. The dashed line separates the stable and the unstable regions as
determined by eq. (2.25).

For the cosine potential, these exponents are calculated in ref. [21]:

$\mu_{k}=\frac{m\mathcal{I}}{\mathrm{K}(\sqrt{\epsilon})}\sqrt{2\kappa^{2}\quantity(\epsilon-\kappa^{2})\quantity(1-\epsilon+\kappa^{2})}$
(2.23)

where

$\kappa\equiv k/m.$

$\mathrm{K}$ is the complete elliptic integral of the first kind, and
$\mathcal{I}$ denotes the following integral:

$\mathcal{I}=\int_{0}^{\pi/2}\frac{\differential{\vartheta}}{1+\quantity(1-2\epsilon+2\kappa^{2})\sin^{2}\vartheta}\frac{2\sin^{2}\vartheta}{\sqrt{\quantity(1+\sin^{2}\vartheta)\quantity[1+(1-2\epsilon)\sin^{2}\vartheta]}}.$
(2.24)

In these expressions $\epsilon$ denotes the ratio of the homogeneous mode
energy density to the height of the potential barrier $2m_{a}^{2}f_{a}^{2}$,
so it satisfies $0<\epsilon<1$. From the Floquet exponent eq. (2.23) we can
immediately read off the condition for $\Re{\mu_{k}}>0$ as

$0<\quantity(\frac{k}{m})^{2}<\epsilon.$ (2.25)

The modes for which the above condition is true are referred to as the modes
inside the instability band. We see that the width of the instability band
depends on the energy in the homogeneous mode. If $\epsilon$ is large, the
homogeneous mode oscillates with a larger amplitude which allows the ALP field
to probe the non-quadratic parts of the potential. This enhances the
parametric resonance effect as we can see in figure 5, where we show a plot of
the Floquet exponents for the periodic potential without the Hubble expansion.
We see that a large $\epsilon$ not only increases the width of the instability
bands, but also renders the parametric resonance more efficient.

In the case of non-zero $H$, the evolution of the mode functions can be
approximated by [24]

$u_{k}(t)\sim u_{k}(t_{i})A_{k}(t)N_{k}(t),$ (2.26)

where $u_{k}(t_{i})$ is the initial condition before parametric resonance
becomes effective, $A_{k}(t)$ encompasses the decay of the mode amplitude due
to the redshift which can be estimated using the WKB approximation as

$A_{k}(t)\propto\omega_{k}^{1/2}(t)a^{-3/2}(t)\mbox{\quad
where\quad}\omega_{k}(t)=\sqrt{\frac{k^{2}}{a^{2}(t)}+m_{a}^{2}},$ (2.27)

and $N_{k}(t)$ is the total growth of the mode due to parametric resonance
which is obtained by integrating the Floquet exponent over time:

$N_{k}(t)=\exp\quantity(\int_{t_{i}}^{t}\differential
t^{\prime}\mu_{k}\quantity(\kappa(t^{\prime}),\epsilon(t^{\prime}))),\quad\kappa(t)=\frac{k}{m_{a}a(t)}.$
(2.28)

In terms of dimensionless variables

$\kappa_{\rm osc}\equiv\frac{k}{m_{a}a_{\rm osc}},\quad\tau\equiv 2H_{\rm
osc}t,\quad\tilde{u}_{k}\equiv\frac{\mu_{k}}{m_{a}},$ (2.29)

eq. (2.28) reads

$\boxed{N_{k}(\tau)=\exp\quantity(\frac{m_{a}}{2H_{\rm
osc}}\int_{\tau_{i}}^{\tau}d\tau^{\prime}\tilde{\mu}_{\kappa}\quantity(\kappa(\tau^{\prime}),\epsilon(\tau^{\prime})))\equiv\exp\quantity(\frac{m_{a}}{2H_{\rm
osc}}\mathcal{B}_{\kappa}(\tau)).}$ (2.30)

By calculating $\mathcal{B}_{\kappa}$ numerically, we find that it becomes at
most $\mathcal{O}(0.5)$. Therefore, the efficiency of the fragmentation is
determined primarily by the hierarchy between the ALP mass and the Hubble
scale at the onset of oscillations. This estimate shows that in the standard
misalignment mechanism with a small initial field value $\phi_{i}\lesssim
f_{a}$, the parametric resonance is not effective. However, in the Large
Misalignment Mechanism, the tiny potential gradient at the top of the cosine
potential delays the onset of oscillations, see section 2.1.1, so that a
hierarchy between $m_{a}$ and $H_{\rm osc}$ renders the parametric resonance
effective. Such a delay also occurs in the Kinetic Misalignment Mechanism,
where the ALP field has a large initial kinetic energy [25, 26] which yields
very efficient parametric resonance [24].

####  2.3.2 Tachyonic instability

Let us again briefly ignore the Hubble expansion, but also assume that the
homogeneous mode is still frozen at the initial field value $\theta_{i}$,
although this is not actually possible when $H=0$. If the initial field value
would be close to $-\pi$, the second derivative of the potential,
$\cos(\theta_{i})$, gets close to $-1$. A negative value of the second
derivative implies a tachyonic instability [46], which leads to growth of
fluctuations, i.e. particle production. This can be seen by assuming
$\cos(\theta_{i})=-1$, with which eq. (2.19) becomes

$\ddot{u}_{k}(\tilde{t})+\left(\tilde{k}^{2}-1\right)u_{k}(\tilde{t})=0\,.$
(2.31)

The solution to this is

$u_{k}(\tilde{t})=u_{k}(\tilde{t}_{0})e^{\pm
i\sqrt{\tilde{k}^{2}-1}(\tilde{t}-\tilde{t}_{0})}\,.$ (2.32)

For $\tilde{k}>1$, i.e. for relativistic modes, this results in oscillation of
the fluctuation. This stays true even if the second derivative of the
potential would be closer to zero, i.e. if the initial field value would be
smaller. For $\tilde{k}<1$, the term in the root is negative, which implies
that the mode grows exponentially.

It is obvious that this analysis does not give a realistic result for the
expanding universe. To find out if this exponential growth could have
cosmological implications, we have to compare this (dimensionless) imaginary
frequency $\sqrt{\tilde{k}^{2}-1}$ with the dimensionless Hubble-factor
$H/m_{a}$, and check if the latter is large compared to the former:

$\frac{\sqrt{\tilde{k}^{2}-1}\,m_{a}}{H(\tilde{t})}=\sqrt{\tilde{k}^{2}-1}\cdot{2\tilde{t}}=\mathcal{O}(1)\cdot\tilde{t}\stackrel{{\scriptstyle!}}{{>}}1\,.$
(2.33)

This shows that the tachyonic instability can only be cosmologically relevant
if $\tilde{t}>1$. We saw that $\tilde{t}=1$ is approximately the dimensionless
time at which the homogeneous mode begins to oscillate in the standard
misalignment mechanism. This means, $\cos(\theta)$ does not stay negative and
exponential growth of the fluctuations will only be possible under certain
conditions, which we will investigate in the next subsection.

Figure 6: Numerical solutions to eq. (2.18) (homogeneous mode) with
$\theta_{i}=\pi/2+0.5$, for which the second derivative of the potential is
negative, and to eq. (2.19) (fluctuations) for different values of
$\tilde{k}$. Naively we could expect exponential growth for $\tilde{k}=0.1$
and $0.5$ due to the tachyonic instability, but this would only apply if the
homogeneous mode would stay frozen. Since it starts to oscillate before
$\tilde{t}=1$, we see no growth of fluctuations. After the onset of
oscillation all modes decay due to Hubble friction.

To summarise, when the homogeneous mode oscillates harmonically near the
minimum and is damped due to the expansion, no exponential growth of
fluctuations due to the tachyonic instability is possible. Figure 6 shows a
numerical solution of the equations of motions, with no growth of
fluctuations, since the homogeneous mode does not stay frozen long enough.
This prediction changes if the oscillations start closer to the top of the
potential.

####  2.3.3 Growth of fluctuations in the Large Misalignment Mechanism

In the Large Misalignment Mechanism introduced in section 2.1.1, for an
initial value of $\theta_{i}$ very close to $\pi$, like $\theta_{i}=\pi-0.01$,
the damping term in eq. (2.18) stays dominant for longer and the homogeneous
mode can still be frozen for $\tilde{t}>1$. Since $\cos(\theta_{i})\approx-1$,
eq. (2.19) approximately becomes

$\ddot{u}_{k}(\tilde{t})+\frac{3}{2\tilde{t}}\dot{u}_{k}(\tilde{t})+\left(\frac{\tilde{k}^{2}}{\tilde{t}}-1\right)u_{k}(\tilde{t})=0\,.$
(2.34)

When the damping term gets unimportant, i.e. $\tilde{t}>1$, a mode begins to
oscillate with frequency $\sqrt{{\tilde{k}^{2}}/{\tilde{t}}-1}$ if the square
root is real. If this becomes imaginary, i.e. $\tilde{k}^{2}/\tilde{t}<1$,
modes show exponential growth. Because of the redshift of physical modes, also
modes with $\tilde{k}>1$ can grow if the oscillations are delayed long enough.
We show this in figure 7.

Due to the large amplitude of the oscillation of the homogeneous mode, which
makes the higher order terms in the expansion of the cosine potential become
important, at later times also parametric resonance leads to growth of
fluctuations, which was analysed in section 2.3.1.

Figure 7: Homogeneous mode and mode functions with similar normalisation to
figure 6. Oscillation of the homogeneous mode is delayed beyond $m_{a}t>1$ due
to the large misalignment $\theta_{i}=\pi-0.01$. The large misalignment makes
it possible for fluctuations to grow exponentially due to the tachyonic
instability already before they start to oscillate. This differs from figure
6, where no growth at early times could happen because oscillation of the
homogeneous mode started too early.

###  2.4 Growth of fluctuations in a non-periodic potential

Now we consider the potential from eq. (2.8) for which, due to the non-
periodicity, it is possible to have large initial field values for the field,
i.e. $\phi_{i}\gtrsim f_{a}$. We are particularly interested in this regime,
since it can involve an efficient parametric resonance and fragmentation of
the ALP field. The motivation for large field values depends on the details of
the actual UV-completion and on pre-inflationary physics. We discuss some of
them in appendix A. Generally, the ALP field can be driven to large values as
a result of a random-walk behaviour during inflation, originating from the
quantum fluctuations of the ALP field, stretched to superhorizon scales during
inflation. These have an amplitude of the order of the inflationary Hubble
scale $H_{I}$ when they leave the horizon.111The variance of the field is in
equilibrium when the average field value is approximately $H_{I}^{2}/m_{a}$.
This shows that for small masses initial field values of $\phi_{i}\gg f_{a}$
originating from this process should be possible even under our assumption
that $H_{I}\ll f_{a}$. [47] Note however that such displacements during
inflation can be noticeable only in the case of a very long duration of
inflation, typically more than $10^{20}$ $e$-folds.

####  2.4.1 Floquet analysis

Analogous to section 2.3.1, for the analytical analysis of parametric
resonance with the given non-periodic potential we will first ignore the
expansion and set $H=0$ and $a=1$. The resulting EOMs for the dimensionless
homogeneous mode and the mode functions are

$\begin{split}0&=\ddot{\theta}(\tilde{t})+\frac{V^{\prime}(\phi)}{m_{a}^{2}f_{a}}\\\
\iff
0&=\ddot{\theta}(\tilde{t})+\theta(\tilde{t})\left[1+\theta^{2}(\tilde{t})\right]^{p-1}\,\\\
\end{split}$ (2.35)
$\begin{split}0&=\ddot{u}_{k}(\tilde{t})+\left[\frac{k^{2}}{m_{a}^{2}}+\frac{V^{\prime\prime}(\phi)}{m_{a}^{2}}\right]u_{k}(\tilde{t})\\\
\iff
0&=\ddot{u}_{k}(\tilde{t})+\left[\frac{k^{2}}{m_{a}^{2}}+(1+\theta(\tilde{t})^{2})^{p-2}\left(1+\theta(\tilde{t})^{2}(2p-1)\right)\right]u_{k}(\tilde{t})\,.\end{split}$
(2.36)

Since the homogeneous mode is oscillating periodically (with period depending
on the initial amplitude), the second derivative of the potential in eq.
(2.36) is periodic in time. Therefore, we can again investigate the evolution
of the fluctuations with the help of the Floquet theorem. (see appendix B).

We show the real part of the Floquet exponents for different initial field
values and momenta in the left panel of figure 8 for $p=-1/2$. The instability
regions show for which values of $\phi_{i}/f_{a}$ and $k/m_{a}$ exponential
growth of the modes happens.

Figure 8: Parametric resonance in the ALP model with non-periodic potential
from eq. (2.8) with $p=-1/2$ and $\phi_{i}=15f_{a}$. The left plot shows the
real part of the Floquet exponents and how the two modes move through the
$(\kappa$,$\phi)$-plane as they redshift and the homogeneous mode gets damped.
The right plot shows numerical results for the absolute value of the
homogeneous mode as a dotted line, and an analytical fit of the envelope in
blue. We also show the absolute values of exponentially growing modes with
$\tilde{k}=4$ and $8$, respectively. Details of this plot are discussed in the
text of section 2.4.1.

It goes beyond the scope of this analysis to attempt a rigorous Floquet
analysis for $H\neq 0$, that involves a time dependent physical momentum in
the equations of motion. Nevertheless, we can use the analysis from the last
subsection to analytically infer which modes will show exponential growth.
Since the field behaves as matter component after the onset of oscillation,
the homogeneous mode roughly dilutes as $\phi\sim t^{-3/4}$, as soon as it
begins to oscillate. To find the beginning of oscillation, we again use eq.
(2.9). For the non-periodic potential this results in
$\tilde{t}_{\mathrm{osc}}=\left(2\sqrt{(1+\phi_{i}^{2})^{p-1}}\right)^{-1}$.
Hence, we conclude that the amplitude of the oscillation of the homogeneous
mode behaves as
$\phi_{\mathrm{ampl.}}({\tilde{t}})=\phi_{i}\left(\tilde{t}/(A\tilde{t}_{\mathrm{osc}})\right)^{-3/4}$
for $\tilde{t}>A\tilde{t}_{\mathrm{osc}}$, where $A=\mathcal{O}(1)$ is
introduced to compensate for the anharmonicity of the oscillation and fit the
numerical result for the homogeneous mode (compare to section 2.1.2). It is
obvious that in radiation era, where $a\sim t^{1/2}$, the dimensionless
physical momentum $\kappa\equiv k/(m_{a}a)$ of one given mode dilutes as
$\kappa\sim t^{-1/2}$. Together with $\phi\propto a^{-3}$ this results in a
curve in the $(\kappa$,$\phi)$-plane on which the different modes move towards
smaller values of the homogeneous mode and physical momentum during the
expansion. From figure 8, we can estimate which modes will grow and how
strongly as different modes probe different parts of the instability region.
We can compare this to numerical solutions for the homogeneous mode and
fluctuations eq. (2.14) and eq. (2.15), which is shown for $p=-1/2$ for two
exemplary modes in figure 8. The left plot shows the real part of the Floquet
exponent in the $(\kappa$,$\phi)$-plane and analytically derived curves that
illustrate how modes “move” through this plane, i.e. what the value of the
amplitude of the oscillation of the homogeneous mode is when a mode has
reached a certain physical momentum. These curves are outside the shown region
at $\tilde{t}_{\mathrm{osc}}$ (at $\phi/f_{a}=15$ and
$\kappa=\tilde{k}/\tilde{t}_{\mathrm{osc}}$) and move towards the origin. The
curves end when $\tilde{t}$ reaches the end of the numerical simulation. The
exponential growth in the modes that is seen on the right hand side
corresponds to the respective time these modes enter the instability region
that is depicted on the left hand side.

####  2.4.2 Tachyonic instabilities

Tachyonic instabilities lead to growth of fluctuations when the homogeneous
mode is still frozen and sits at a point in the potential at which the second
derivative is negative. In the Large Misalignment Mechanism this was an
important effect, since the top of the potential, where the field is fine-
tuned to, has the largest negative value of the second derivative of the
cosine potential. In the non-periodic potential, depending on the value of
$p$, there may be no region of tachyonic instabilities, for example for
$p=0.5$, while for $p=-0.5$ the second derivative is negative for most field
values, but asymptotes to zero from below for larger ones, rendering the
tachyonic instability increasingly unimportant for large initial field values.

###  2.5 Regions of linear and non-linear growth

Eq. (2.15) is based on linearization of the equations of motion, assuming
small fluctuations. This neglects the interaction between different modes and
the backreaction of the modes onto the homogeneous mode. Ignoring these
effects is justified only as long as the fluctuations remain small relative to
the homogeneous mode, which is the case for sufficiently small initial field
values $\phi_{i}/f_{a}$. In contrast, for larger values of $\phi_{i}/f_{a}$,
non-linear effects and the backreaction onto the homogeneous mode eventually
become significant. This usually causes the rapid, non-perturbative process of
fragmentation [20], during which the energy from the homogeneous mode is
completely transferred into fluctuations. For the cosine potential, numerical
simulations show that this only happens for initial field values tuned close
to the top or if the field initially carries enough kinetic energy to cross
several maxima of the potential [22, 24].

To determine the transition between the two regimes, we track the components
of the scalar field energy density in the homogeneous mode and in
fluctuations.222Since the ratio of these two quantities,
$\rho_{\mathrm{fluc}}/\rho_{0}$, is independent of $m_{a}$ and $f_{a}$, we
only have to evaluate them for different initial values
$\theta_{i}=\phi_{i}/f_{a}$ of the homogeneous mode.

$\displaystyle\rho_{0}=\frac{1}{2}\dot{\phi}(t)^{2}+V(\phi)\,,$ (2.37)
$\displaystyle\langle\rho_{\mathrm{fluc}}\rangle=\frac{1}{2}\int\frac{d^{3}k}{(2\pi)^{3}}\left[\lvert\dot{u}_{k}(t)\rvert^{2}+\left(\frac{\vec{k}^{2}}{a^{2}}+V^{\prime\prime}(\phi)\right)\lvert
u_{k}(t)\rvert^{2}\right]\,.$ (2.38)

Note that the mode functions grow in the unstable phase, whereas at
sufficiently late times, after leaving the instability regions and becoming
non-relativistic, they behave as matter and decay in the same way as the
homogeneous mode. As initial conditions for the mode functions at the onset of
oscillation we take

$u_{k}(\tilde{t}_{\mathrm{osc}})=\frac{\sqrt{2\pi^{2}\,2.1\cdot
10^{-9}}}{k^{3/2}}\phi_{i}\,,$ (2.39)

which we will justify in section 3.1. As low momentum cutoff, we take the mode
that enters the horizon at the end of the simulation and the large momentum
cutoff at $k/a_{\mathrm{osc}}=5m_{a}$, since larger momenta will not be
affected by parametric resonance. We stop the simulation when
$\rho_{\mathrm{fluc}}/\rho_{0}$ grows to unity since we can be sure that the
linear description fails in that case. This is labelled as non-linear regime.
If $\rho_{\mathrm{fluc}}/\rho_{0}<1$ at sufficiently late times of the
simulation, we are in the linear regime. The initial field value that
separates the linear regime from the non-linear one is referred to as the
critical field value, $\theta_{\mathrm{crit}}=\phi_{\mathrm{crit}}/f_{a}$.

If we numerically evaluate $\theta_{\mathrm{crit}}$ for several $p$, we find
that it decreases for flatter potentials, i.e. as $p$ gets more negative,
while it diverges if $p$ gets close to unity, which corresponds to a quadratic
potential where no parametric resonance is possible. The result is shown in
figure 9.

Figure 9: Critical field value $\phi_{\mathrm{crit}}$ at which
$\rho_{\mathrm{fluc}}/\rho_{0}$ becomes larger than 1 for all
$\phi_{i}>\phi_{\mathrm{crit}}$, with the initial conditions given in section
2.5, for different $p$ values in the non-periodic potential. The error-bars
reflect the resolution of the numerical simulation. As indicated by the arrow,
the rightmost point only gives a lower bound on the actual critical value.

When we repeat this analysis for the Large Misalignment Mechanism and
investigate for which parts of the parameter space in the
$(m_{a},f_{a})$-plane fragmentation is expected, we find that we need a fine-
tuning stronger than $\absolutevalue{\pi-\theta_{i}}\approx 10^{-9}$.

Our goal is to investigate the amplification of ALP fluctuations and determine
whether this can lead to observable effects. The crucial quantity that we will
focus on is the energy density power spectrum, which we will discuss in the
next section. It will quantify the inhomogeneities of the field and feed into
the formalism that captures gravitational collapse of the overdense regions,
which we discuss in section 4.

##  3 The power spectrum

We start this section by specifying the initial conditions for the
fluctuations of the field. In order to properly describe the horizon entry of
the curvature perturbations, we include the corresponding source terms in the
linearized evolution equations. This is explained in section 3.1, where both
adiabatic and isocurvature fluctuations are discussed. In section 3.2 we
follow the subsequent evolution of ALP fluctuations and compute the resulting
energy density power spectra. In the cases of efficient parametric resonance,
the linear approximation eventually breaks down and we “glue” the linear
evolution to a non-linear lattice simulation of the field dynamics, as it was
done in [30, 31]. We close with section 3.3, where the role of oscillons that
are formed in the non-linear regime is explained.

###  3.1 Initial conditions for the fluctuations

The ALP field can have fluctuations of both adiabatic and isocurvature type,
which are discussed in this subsection.

Adiabatic fluctuations of ALPs are sourced by curvature perturbations. This
process can be described by including scalar metric perturbations in the field
evolution equations. In the Newtonian gauge, the perturbed metric has the form

$ds^{2}=[1+2\Phi(t,\vec{x})]dt^{2}-a^{2}(t)[1-2\Phi(t,\vec{x})]d\vec{x}^{2},$
(3.1)

where $\Phi(t,\vec{x})$ denotes the scalar metric perturbation field. In the
radiation-dominated era, the Fourier modes $\Phi_{\vec{k}}$ of the scalar
perturbations are related, in the super-horizon limit, to those of curvature
perturbations $R_{\vec{k}}$ via $\Phi_{\vec{k}}=-(2/3)R_{\vec{k}}$. Inflation
predicts a power-law form for the power spectrum of the curvature
perturbations [48, 22],

$\langle
R_{\vec{k}}R^{\star}_{\vec{k^{\prime}}}\rangle=(2\pi)^{3}\delta^{(3)}(\vec{k}-\vec{k^{\prime}})\frac{2\pi^{2}}{k^{3}}\Delta^{2}_{R}(k),\>\>\>\>\>\>\>\>\>\Delta^{2}_{R}(k)=A_{s}\Bigl{(}\frac{k}{k_{\star}}\Bigr{)}^{n_{s}-1}\,.$
(3.2)

The amplitude $A_{s}=2.1\cdot 10^{-9}$ at $k_{\star}=0.05/\rm{Mpc}$ and the
spectral tilt $n_{s}=0.965$ are inferred from the PLANCK CMB-measurements
[49]. For simplicity, the slight spectral tilt can be ignored in our
calculations. The solution for the dynamics of curvature perturbations in the
radiation era, derived in the linear theory, has the form333Note that the
backreaction of the ALP field onto the curvature perturbations is negligible.

$\Phi_{\vec{k}}(t_{k})=3\Phi_{\vec{k},i}\left(-\frac{\cos(t_{k})}{t_{k}^{2}}+\frac{\sin(t_{k})}{t_{k}^{3}}\right)\,,$
(3.3)

where $\Phi_{\vec{k},i}$ denotes the super-horizon value of the mode function
and $t_{k}$ is a momentum dependent, dimensionless time variable defined as
$t_{k}\equiv k/(\sqrt{3}aH)$ [50]. As can be seen, the modes $\Phi_{\vec{k}}$
decrease as $a^{-2}$ once they are deep inside the horizon, $t_{k}\gg 1$.

Returning to the ALP field, due to the smallness of the fluctuations $\Phi$
and $\delta\varphi$, at least at early times, we keep only terms that are of
zeroth or of linear order in any of these perturbations in the field evolution
equations. In this way one arrives at [23, 22]

$\boxed{\ddot{u}_{k}(t)+3H\dot{u}_{k}(t)+\left(\frac{k^{2}}{a^{2}}+V^{\prime\prime}(\phi)\right)u_{k}(t)=4\dot{\phi}(t)\dot{\Phi}_{k}(t)-2V^{\prime}(\phi)\Phi_{k}(t).}$
(3.4)

The left hand side of the above equation is the same as the one in non-
perturbed spacetime in eq. (2.15), while the right hand side introduces source
terms from the curvature perturbations. These terms can lead to fluctuations
in the ALP field, even if the latter are initially zero. The sourcing mostly
occurs close to the time at which the corresponding curvature mode enters the
horizon, i.e. when $H\sim k/a$. The simulation should therefore be started
earlier than $\tilde{t}_{\mathrm{osc}}$, such that all modes of interest are
superhorizon.

Given that the ALP field was present during inflation, it unavoidably carries
isocurvature fluctuations of the size of the inflationary Hubble scale $H_{I}$
(see ref. [51]). More specifically, the scale-invariant power spectrum of the
fluctuations can be determined from the Bunch-Davies mode functions and is
given by

$\Delta^{2}_{\delta\phi}(k)=\left(\frac{H_{I}}{2\pi}\right)^{2},$ (3.5)

in the superhorizon limit.

In the following it will be assumed that the isocurvature fluctuations are
smaller compared to the adiabatic ones and, thus, they will be neglected. In
this case one can work with mode functions $u_{k}$ that are initially zero. It
was also verified that similar results are obtained if the source terms in eq.
(3.4) are set to zero and, instead, reasonable initial values for the ALP mode
functions are chosen, as it was done in section 2.4.

###  3.2 Growth of fluctuations and the density power spectrum

In this section we examine the amplification of ALP fluctuations. The
discussion is split into two parts. In section 3.2.1, we consider the case of
$\theta<\theta_{\rm crit}$, where a linear analysis of the dynamics is
performed. The non-linear regime for $\theta>\theta_{\rm crit}$ is
investigated in section 3.2.2.

To link the evolution of the fluctuations to observable quantities it is
convenient to compute the energy density (contrast) power spectrum. The
density contrast is defined as

$\delta(\vec{x})\equiv\frac{\rho(\vec{x})-\overline{\rho}}{\overline{\rho}}\,,$
(3.6)

where all quantities are time-dependent. $\overline{\rho}$ is the averaged
energy density and in the quantum picture we identify this as
$\expectationvalue{\rho}$. The Fourier transform of $\delta(\vec{x})$ is
$\delta_{\vec{k}}=\rho_{\vec{k}}/\overline{\rho}$, where $\rho_{\vec{k}}=\int
d^{3}x\,\rho(x)e^{i\vec{x}\cdot\vec{k}}$ is the Fourier transform of
$\rho(\vec{x})$ and we demand $k\neq 0$. Again we will assume isotropy in
momentum space, i.e. care only for the magnitude of $\vec{k}$. The
dimensionful density contrast power spectrum is defined as [52]

$P_{\delta}(k)\equiv\frac{\expectationvalue{\delta_{k}^{2}}}{\mathcal{V}}=\frac{\expectationvalue{\absolutevalue{\rho_{k}}^{2}}}{\mathcal{V}\expectationvalue{\rho}^{2}}\,,$
(3.7)

where $\mathcal{V}$ is a large volume over which we integrate to keep the
quantities finite. We will see that it drops out in the end. We also introduce
the dimensionless density contrast power spectrum which is defined as

$\Delta_{\delta}(k)=\frac{k^{3}}{2\pi^{2}}P_{\delta}(k)\,.$ (3.8)

####  3.2.1 Linear regime

Numerically calculating eq. (3.4) for a range of momenta allows one to track
the density power spectrum $\Delta_{\delta}(k)$. Note that the term
$\dot{\Phi}$ in eq. (3.4) can be re-expressed as
$m_{a}\frac{t_{k}}{\tilde{t}}\frac{d\Phi_{\vec{k}}}{dt_{k}}$, where the
derivative with respect to $t_{k}$ can be found from eq. (3.3). This makes it
possible to express eq. (3.4) fully in terms of dimensionless quantities and
arrive at an EOM that is independent of $m_{a}$ and $f_{a}$.444Note that
$f_{a}$ is again a common factor in all terms of eq. (3.4). The enhancement of
the fluctuations thus depends only on the values of $p$ and $\theta_{i}$.

In the linear regime the density contrast perturbation can be computed using
the following expression [22],

$\delta_{\vec{k},\mathrm{lin.,\,pert.}}=\frac{\dot{\phi}\dot{u}_{\vec{k}}+V^{\prime}(\phi(t))u_{\vec{k}}-\dot{\phi}^{2}\Phi_{\vec{k}}}{\frac{1}{2}\dot{\phi}^{2}+V(\phi(t))}\,.$
(3.9)

In figure 10 we illustrate the dimensionless power spectrum, computed using
the expression above, at several scale factors for $\theta_{i}=3$ and
$p=-1/2$. The spectra are shown with solid lines. On the horizontal axis is
the comoving momentum, in units of the ALP mass multiplied by the scale factor
at $H_{\mathrm{osc}}$. As expected, the spectrum is peaked around the momenta
corresponding to the most unstable modes.

Figure 10: Dimensionless density contrast power spectra $\Delta_{\delta}(t,k)$
for ALP dark matter with $\phi_{i}=3f_{a}$ and $p=-1/2$ (solid lines) and for
standard CDM (dashed lines), shown at different values of $a$, considering
radiation dominated era, $a\sim t^{(1/2)}$, and adiabatic initial conditions.
The horizontal axis shows the comoving momentum $k/(m_{a}a_{\rm osc})$. The
curves converge at low momenta, while at high momenta the ALP spectrum is
suppressed.

It is instructive to compare the evolution of the density power spectrum to
the one of standard (pressureless) CDM. The latter has the following
analytical form during the radiation-dominated era [53]

$\delta_{k,\rm CDM}=C_{2}+9\Phi_{k,i}\ln[t_{k}],$ (3.10)

where $C_{2}=3\Phi_{k,i}/2$ for adiabatic fluctuations. As can be seen, the
density contrast grows logarithmically with the scale factor.

The solution (eq. (3.10)) is shown in figure 10 with dashed lines at the same
scale factors as the solid lines of the corresponding colour. One can
distinguish three momentum ranges. At low momenta the ALP power spectrum
converges to the one for CDM. At intermediate momenta $k\sim
m_{a}a_{\mathrm{osc}}$, the ALP power spectrum is enhanced due to parametric
resonance. Finally, at high momenta one observes a strong suppression of the
ALP power spectrum. More details about the evolution of density perturbations
for wave dark matter can be found in ref. [53].

####  3.2.2 Nonlinear regime

If $\theta_{i}>\theta_{\rm crit}$, the linear approximation is expected to be
valid only at early times. Once the fluctuations become large, their
backreaction onto the homogeneous mode becomes important, which slows down
their growth and, in most cases, leads to the fragmentation of the field.

In order to extend the discussion to this non-linear regime, we “glue” the
linearized evolution to the lattice simulation that solves the full EOM eq.
(2.2). More specifically, the linear simulation is stopped when the
fluctuations are still small, $\langle\delta\theta^{2}\rangle=10^{-4}$. The
values of the field and its derivative, as well as the power spectra computed
using the mode functions, are then extracted and used to initialise a Gaussian
random field on the lattice. More details about the lattice setup can be found
in appendix C as well as in ref. [30]. The field is evolved on a spatial
lattice with $512^{3}$ grid points, periodic boundary conditions and a fixed
comoving volume, using a standard leap-frog algorithm. In principle, it is
possible to incorporate metric perturbations into the lattice simulation, in
the form of an additional scalar field. However this would require the Fourier
transformation of the field after each time step, to evolve the mode functions
in momentum space according to eq. (3.3). On the other hand, in the regime
when fragmentation takes place, the metric perturbations become unimportant as
the field is no longer correlated with the metric perturbations. To this end
we switch off the metric perturbations once the simulation starts.

The ALP dynamics after fragmentation usually involves the formation of long-
lived compact objects, called oscillons [54, 55]. These are coherent and
localised configurations of the scalar field, held together by the attractive
self-interactions. They can be very long-lived, even though there is no
conserved charge that would ensure their stability and prevent them from
decaying (see also ref. [56] for the approximate conservation of the adiabatic
invariant of oscillons). As a consequence, some fraction of the total energy,
denoted by $F_{\mathbb{OSC}}$, is contained at late times in oscillons, while
the rest is in the form of unbound fluctuations.555The subscript
$\mathbb{OSC}$ will refer to oscillons, and should not be confused with the
subscript $\mathrm{osc}$, referring to the onset of oscillations. The value of
$F_{\mathbb{OSC}}$ depends on $p$, as well as on the initial field value
$\theta_{i}$, and in some cases can reach around $70\%$ [57].

The separation into the contributions mentioned above can be observed in the
field power spectrum. At sufficiently late times, the spectrum splits into two
parts: one at high momenta, due to oscillons, and one at lower momenta due to
fluctuations. The oscillon peak "blue-shifts" with time to higher momenta, as
the oscillons maintain fixed physical size. In contrast, the low-momentum
contribution from the fluctuations eventually becomes time-independent,
maintaining constant comoving momenta. Such behaviour is expected due to the
Hubble expansion since the self-interactions become weaker as the ALP field is
driven to the vicinity of the quadratic minimum of the potential. The time for
this “freeze-out” depends on the initial field value, i.e. it is larger for
larger values of $\theta_{i}$.

As a result of this behaviour, the oscillon peak gets well separated from the
low-momentum fluctuations at late times. This is illustrated in figure 11,
where several snapshots of the dimensionless field power spectrum
$\Delta_{\theta}(k)$ are shown. The multiplication by the scale factor $a^{3}$
is done for convenience, to scale out the overall dilution due to the Hubble
expansion. The curves are shown at scale factors in the range from
$19a_{\mathrm{osc}}$ to $32a_{\mathrm{osc}}$. $\theta_{i}=4$ and $p=-1/2$ were
used for this plot, although the pattern is rather generic. The height of the
oscillon peak shows some oscillatory behaviour but on average does not change
much within the simulation time. The relative strength of the two
contributions depends on $F_{\mathbb{OSC}}$.

Figure 11: Snapshots of the dimensionless field power spectrum
$\Delta_{\theta}(t,k)$ for $\theta_{i}=4$ and $p=-1/2$, demonstrating the
oscillon peak at high momenta and the fluctuation peak at small momenta. The
colours from light to dark blue correspond to scale factors
$a/a_{\mathrm{osc}}$ ranging from $19$ to $32$. The inset compares the linear
(dashed lines) and the lattice (solid lines) evolution at scale factors
$a/a_{\mathrm{osc}}$ between $4$ and $24$, illustrating how non-linear effects
slow down and spread out the growth of low-momentum fluctuations.

As in the linear regime, we are interested in the energy density power
spectrum. Here it can be extracted from the lattice simulation by explicitly
computing the energy density of the scalar field and evaluating its
correlation functions. It is observed that the density power spectrum is
dominated by the contribution from oscillons, in the form of a strong peak in
momentum space. This is illustrated in figure 12, where the solid lines show
the dimensionless density contrast power spectrum at scale factors from
$14a_{\mathrm{osc}}$ to $a=28a_{\mathrm{osc}}$ for $\theta_{i}=5$ and
$p=-1/2$. A semi-analytical expression for this peak was found in ref. [57],
where it was parametrised in terms of an extended Poisson distribution,

$P_{\delta}(k;a)=\frac{F_{\mathbb{OSC}}^{2}}{n_{\mathbb{OSC}}}\frac{\langle
M_{\mathbb{OSC}}^{2}\exp\\{-r_{\mathbb{OSC}}^{2}k^{2}/(2a^{2})\\}\rangle}{\langle
M_{\mathbb{OSC}}\rangle^{2}}K(k/k_{c})\,$ (3.11)

with $K(x)=1-(2/x)^{2}\sin^{2}(x/2)$. Here $n_{\mathbb{OSC}}$,
$r_{\mathbb{OSC}}$ and $M_{\mathbb{OSC}}$ are the average number density (per
comoving volume), the (physical) radius and the mass of the oscillons. $k_{c}$
is a cut-off scale for the oscillon power spectrum, which is determined by
fitting the numerical data. It is related to the horizon scale at the oscillon
formation.

Figure 12: Snapshots of the energy density contrast power spectrum
$\Delta_{\delta}(t,k)$ for $\theta_{i}=5$ and $p=-1/2$ at the scale factors
when oscillons have already formed, ranging from $a=14a_{\mathrm{osc}}$ to
$a=28a_{\mathrm{osc}}$. The solid lines show the full (non-linear) power
spectrum, dominated by oscilllon contribution. The dashed lines correspond to
the Gaussian approximation for the power spectrum, according to eq. (3.13)
which, at small momenta, is dominated by the contribution from fluctuations.

###  3.3 The role of oscillon decay

The lifecycle of an oscillon can be separated into a long phase of slow
“evaporation”, via the emission of scalar radiation, and a very rapid “decay”
in the end, when the oscillon energy is converted into a burst of scalar waves
[58]. Power-law potentials with $p<1$, considered in this work, are known to
support very long-lived oscillons. In refs. [19, 59], the lifetime was
estimated to be $\tau_{\mathbb{OSC}}\sim 10^{10}m_{a}^{-1}$ for $p=-1/2$,
$\tau_{\mathbb{OSC}}\sim 10^{8}m_{a}^{-1}$ for $p=1/2$, and even longer for
$p\rightarrow 1$. In contrast, in the case of the QCD axion the corresponding
objects, often referred to as axitons, can survive only at most for $\sim
10^{3}m_{a}^{-1}$ [60]. The very long lifetimes of oscillons make it
impossible to extend the lattice simulations to their decay stage.
Nevertheless, they do eventually decay and, in a large region of the parameter
space given by

$m_{a}>10^{-19}\,\mathrm{eV}\left(\frac{\widetilde{\tau}_{\mathbb{OSC}}}{10^{8}}\right),$
(3.12)

with $\tilde{\tau}_{\mathbb{OSC}}=m\tau_{\mathbb{OSC}}$, this happens before
the epoch of matter and radiation equality. It is thus important to understand
how the decay is expected to impact the power spectra. This is discussed
below, with a more detailed analysis of the decay process to be included in a
future publication.

For simplicity we ignore the slow evaporation of oscillons, assuming that
their energy is radiated away in the very end. In this last stage the oscillon
peak in the field power spectrum will already be separated by several orders
of magnitude from the low-momentum fluctuation modes. The decay of oscillons
usually produces ALPs with wavelengths $(k/a)\gtrsim 0.1m_{a}$, which simply
free-stream away [58]. It is thus expected that this process does not affect
the field power spectrum at much smaller momenta.

We want to estimate the energy density power spectrum at late times, after the
decay. At this stage the self-interactions are unimportant and, assuming that
the field power spectrum is nearly Gaussian, it can be expressed as [60]

$P_{\delta}(q)=2\frac{\int\frac{d^{3}k}{(2\pi)^{3}}n_{k}n_{|k+q|}}{\left(\int\frac{d^{3}k}{(2\pi)^{3}}n_{k}\right)^{2}}\,,$
(3.13)

in the deep non-relativistic limit. Here

$n_{k}+\frac{1}{2}=\frac{|\dot{u}_{k}|^{2}+\omega_{k}^{2}|u_{k}|^{2}}{2\omega_{k}}$
(3.14)

denotes the occupation numbers. We discuss this in more detail in appendix D,
where we demonstrate that the numerically extracted power spectrum after the
decay of oscillons is indeed well approximated by (3.13).

In figure 12 we show several power spectra, computed using the above Gaussian
expression, in dashed lines. The solid lines in that figure correspond to the
full non-linear power spectra from the same simulation. As can be observed,
$P_{\delta}(q)$ at small momenta becomes constant after some time. From the
form of eq. (3.13) we can conclude that, by these times, it is dominated by
the contribution from fluctuations at small $k$’s, rather than by the time-
dependent oscillon contribution. This is expected to hold also after the decay
of oscillons and, therefore, $P_{\delta}(q)$ should not change much any
further. This observation allows us to estimate the power spectrum at small
momenta and late times by its time-independent value from the simulation.

To summarise, even though our simulations do not resolve the dynamics of the
oscillons until their late decay, we can make predictions for the low momentum
part of the energy power spectrum after their decay.

##  4 ALP halos

Having computed the density power spectra for ALPs in a non-periodic
potential, in this section we study how the overdense regions, corresponding
to the amplified fluctuations of the ALP field, decouple from the Hubble flow
and form bound objects. We have already briefly mentioned the formation of
long-lived oscillons, which are sustained by the self-interactions. Our focus
of this section is on the more important gravitational interaction, and the
formation of objects such as halos and solitons.

Most of the structure formation takes place in the matter-dominated era, when
the backreaction of the ALPs on the metric perturbations can no longer be
neglected. Our analysis is based on the Press-Schechter (PS) formalism [36].
We start by reviewing the spherical collapse model from ref. [61] and the
Meszaros equation [62] in section 4.1. Modifications to the formalism,
relevant at the length scales associated to fragmentation, are discussed as
well. In section 4.2 we apply the formalism to our model and, using the power
spectra of the density field from the previous section, compute the halo mass
functions. Virialization of the halos is summarised in section 4.3, where we
also find the halo spectrum.

###  4.1 Spherical collapse

We follow ref. [61] and consider a spherical overdense region of matter in the
homogeneous background of radiation. The equation of motion for the shell of
radius $r$ is given by

$\ddot{r}=-\frac{8\pi G}{3}\rho_{R}r-\frac{GM}{r^{2}}\,,$ (4.1)

where $\rho_{R}$ refers to the energy density of the radiation component and
$M$ is the conserved mass of the matter inside the shell.

As it was done in ref. [61], this equation is investigated by switching to
conformal time $d\eta=dt/a$ and introducing the following parameterisation of
$r$:

$r\equiv a(\eta)\xi(\eta)R\,,$ (4.2)

where $R$ is the comoving size of the spherical region and $\xi(\eta)$ is the
deceleration parameter. The latter is set to $\xi=1$ initially and describes
how an overdense region decouples from the Hubble flow. One also expresses the
mass in terms of the average matter density $\rho$ and the initial over-
density parameter $\delta_{\mathrm{in}}$,

$M=\frac{4\pi}{3}\rho_{\mathrm{EQ}}a_{\mathrm{EQ}}^{3}(1+\delta_{\mathrm{in}})R^{3}\,,$
(4.3)

Note that the term $(\rho a^{3})$ remains constant as soon as the ALPs become
non-relativistic and, in the above expression, it is evaluated at the time of
matter-radiation equality, as denoted by the index $\rm EQ$. For simplicity
the universe is taken to be flat, consisting only of the matter and radiation
components. The scale factor can then be written as
$a(\eta)=a_{\mathrm{EQ}}\left[2(\eta/\eta_{*})+(\eta/\eta_{*})^{2}\right],$
where $\eta_{*}^{-2}=2\pi G\rho_{\mathrm{EQ}}a_{\mathrm{EQ}}^{2}/3$ [61].
Inserting everything into eq. (4.1) and switching from the time variable to
the scale factor $x=a/a_{\rm EQ}$, one arrives at a dimensionless equation for
the deceleration parameter [61]:

$x(1+x)\frac{d^{2}\xi}{dx^{2}}+\left(1+\frac{3}{2}x\right)\frac{d\xi}{dx}+\frac{1}{2}\left(\frac{1+\delta_{\mathrm{in}}}{\xi^{2}}-\xi\right)=0\,.$
(4.4)

Solving this equation numerically, with $\xi=1$ and $d\xi/dx=0$ as initial
conditions, one can find the scale factors at the turnaround
$x_{\mathrm{ta}}$, after which the shell starts to shrink, and the collapse
scale factor $x_{\mathrm{coll}}$. The condition for the turnaround is given by
$\dot{r}=0$ which, using the definition (4.2), can be expressed as
$\dot{\xi}=-\xi H$. In terms of the dimensionless variable $x$ this reads as

$\left.\frac{d\xi}{dx}\right\lvert_{x_{\mathrm{ta}}}+\left.\frac{\xi}{x}\right\lvert_{x_{\mathrm{ta}}}=0\,.$
(4.5)

We find numerically that this is fulfilled at

$\delta_{\mathrm{in}}\cdot x_{\rm turn}\approx 0.71\equiv C_{x}\,,$ (4.6)

where we will use the last definition later. We also evaluated the product of
$\delta_{\mathrm{in}}$ with the numerically extracted collapse scale factor
$x_{\mathrm{coll}}$, defined by $r(x_{\mathrm{coll}})\rightarrow 0$ or
$\xi(x_{\mathrm{coll}})\rightarrow 0$. This gives

$\delta_{\mathrm{in}}\cdot x_{\mathrm{coll}}\approx 1.13\,.$ (4.7)

Linear growth: The actual over-density parameter of the spherical region
$\delta$ will grow during the collapse. It is important to investigate how
this parameter evolves in time, under the assumption of it being small,
$\delta\ll 1$. Since the mass inside the collapsing shell remains constant at
value given by eq. (4.3), one has

$(1+\delta_{\mathrm{in}})=\xi^{3}(x)(1+\delta(x))\,.$ (4.8)

If we now assume $\delta_{\mathrm{in}}\ll 1$, the left-hand side in the above
equation can be set to unity. Solving this for $\xi$ gives
$\xi=\left(1/(1+\delta)\right)^{1/3}$, which can be expanded to
$\xi=1-\delta/3+\mathcal{O}(\delta^{2})$. Inserting $\xi=1-\delta/3$ into eq.
(4.4) we arrive at the Meszaros equation [62] for $\delta$:

$x(1+x)\frac{d^{2}\delta}{dx^{2}}+\left(1+\frac{3}{2}x\right)\frac{d\delta}{dx}-\frac{3}{2}\delta=0\,.$
(4.9)

For a given positive initial over-density $\delta_{\mathrm{in}}$ and initial
condition $\dot{\delta}_{i}\approx 0$, the analytical solution of the Meszaros
equation is

$\delta(x)=\delta_{\mathrm{in}}\left(1+\frac{3}{.}{2}x\right)$ (4.10)

Plugging eq. (4.7) into this result, and assuming that $\delta_{\mathrm{in}}$
is small, such that $1.1/x_{\mathrm{coll}}\ll 1$, gives the well-known value
for the over-density in the linearized theory at which the collapse happens:

$\delta_{c}\equiv\delta(x_{\mathrm{coll}})=\frac{1.1}{x_{\mathrm{coll}}}\left(1+\frac{3}{2}x_{\mathrm{coll}}\right)\approx
1.69.$ (4.11)

Jeans length and scale-dependent growth: The linear growth of density
perturbations in eq. (4.9) is scale independent, which is true for CDM
perturbations. However, in the case of ALP DM, the growth is modified at high
momenta, above a certain Jeans length. To account for these effects due to the
gradient pressure, the Meszaros equations can be modified to [37]

$\boxed{x(1+x)\frac{d^{2}\delta_{k}}{dx^{2}}+\left(1+\frac{3}{2}x\right)\frac{d\delta_{k}}{dx}-\Bigl{[}\frac{3}{2}-\frac{\tilde{k}^{4}}{x}\Bigr{]}\delta_{k}=0.}$
(4.12)

where $\tilde{k}$ is defined in eq. (2.17). Modes with $\tilde{k}\ll 1$ behave
as in the case of CDM, while modes with $\tilde{k}\gg 1$ start growing only
after $x\sim\tilde{k}^{4}$, resulting in a suppression of the power spectrum
at high momenta.

###  4.2 The Press-Schechter formalism and halo mass function

The standard Press-Schechter formalism has been introduced in ref. [36] and
was modified since then (see ref. [63], ref. [64] and references therein). The
original Press-Schechter formalism was developed for dark matter that has been
pressureless for the whole history of the universe, with scale independent
power spectrum that has been imprinted by inflation.

In the PS formalism one introduces the smoothed over-density field
$\delta_{R}(\vec{x})$,

$\delta_{R}(\vec{x})=\int
d^{3}x^{\prime}W_{R}(\vec{x}-\vec{x^{\prime}})\delta(\vec{x^{\prime}})\,,$
(4.13)

where $W_{R}(\vec{x})\equiv W(\vec{x},R)$ is a window function that goes to
zero very fast if $|\vec{x}|>R$. An intuitive choice would be the so-called
spherical top-hat, $W_{\mathrm{sph.\,TH}}(\vec{x},R)=3/(4\pi
R^{3})\Theta(R-\vec{x}),$ where $\Theta(x)$ is the Heaviside step function.
$R$ is called the smoothing scale.

One can evaluate the variance of the smoothed density contrast,
$\sigma_{R}^{2}$. A quick calculation gives

$\sigma_{R}^{2}\equiv\frac{1}{\mathcal{V}}\int
d^{3}x\expectationvalue{\delta_{R}(\vec{x})^{2}}=\frac{1}{2\pi^{2}}\int
dk\,k^{2}\absolutevalue{\tilde{W}_{R}(\vec{k})}^{2}P(k)\,,$ (4.14)

where we used $\mathcal{V}^{-1}P(k)\equiv\expectationvalue{\delta(k)^{2}}$ and
introduced the Fourier transform of the window function
$\tilde{W}_{R}(\vec{k})\equiv\tilde{W}(\vec{k},R)=\int
d^{3}xW(\vec{x},R)e^{-i\vec{k}\cdot\vec{x}}$.

Under the assumption of initially small over-densities, one can relate the
smoothing scale $R$ directly to a halo mass $M$ by setting $M=\rho
V_{R}R^{3}$, where $V_{R}$ is a volume factor that depends on the window
function, where for the spherical top-hat it is given by $4\pi/3$. This
implies that $\sigma_{R}=\sigma(R)$ basically becomes a function of $M$, hence
we can write $\sigma(M)$.

Assuming that $\delta_{R}$ follows a Gaussian distribution (for any value of
$R$) and identifying $\sigma_{R}$ as its variance, the probability to find a
fluctuation in the range $[\delta_{R},\delta_{R}+d\delta_{R}]$ is [36]

$F_{\mathrm{sm}}(\delta_{R};R)=\frac{1}{\sqrt{2\pi\sigma_{R}^{2}}}\exp\left(-\frac{\delta_{R}^{2}}{2\sigma_{R}^{2}}\right)\,.$
(4.15)

The fraction of collapsed objects with mass $>M$ is then given by

$F(>M)=\frac{1}{\sqrt{2\pi}}\int_{\delta_{c}/\sigma_{R}}^{\infty}e^{-t^{2}/2}dt\equiv\mathrm{erfc}[\delta_{c}/\sigma_{R}]\,,$
(4.16)

where $\mathrm{erfc}$ denotes the error function and $\delta_{c}\approx 1.69$
is the critical density at which collapse happens. Due to the one-to-one
correspondence between $R$ and $M$, $\sigma_{R}$ can be understood as a
function of $M$. Thus one can calculate:

$\frac{dF}{dM}=-\frac{1}{\sqrt{2\pi}}\frac{\delta_{c}}{\sigma_{R}^{2}}\frac{d\sigma_{R}}{dM}\exp\left(\frac{-\delta_{c}^{2}}{2\sigma_{R}^{2}}\right)\,.$
(4.17)

This is the fraction of collapsed objects in the range $[M,M+dM]$. The
fractional energy density in these objects is then
$\rho_{0}\cdot\frac{dF}{dM}$, while finally the number density of objects
inside the given mass range is $\rho_{0}\cdot\frac{dF}{dM}/M$. We arrive at
[36]

$\frac{dn}{dM}=-\frac{\rho_{0}}{M}\sqrt{\frac{2}{\pi}}\frac{\delta_{c}}{\sigma_{R}^{2}}\frac{d\sigma_{R}}{dM}\exp\left(\frac{-\delta_{c}^{2}}{2\sigma_{R}^{2}}\right)\,.$
(4.18)

In this result, a seemingly ad hoc factor of 2 is introduced. This is done
because the Gaussian distribution of $\delta_{R}$ also includes negative over-
densities. This region is not accounted for when evaluating eq. (4.16), but
nevertheless these under-densities will be included in collapsed halos since
over-dense regions will be surrounded by under-dense regions which are dragged
into the halo by the higher gravitational potential in the collapsing regions.
The factor of 2 can be rigorously derived in the excursion set formalism [64].

Figure 13: Matter power spectrum today for CDM (black line) vs ALP DM of mass
$m_{a}=10^{-18}\mathrm{eV}$ with $p=-1/2$ (left) and $p=1/2$ (right).
Different colours correspond to different values of the initial misalignment
angle $\theta_{i}=\phi_{i}/f_{a}$. The power spectra agree with CDM on large
scales, are enhanced on scales at which fluctuations were amplified, and are
suppressed on small scales due to wave effects.

To compute the halo mass function, we start with the density contrast power
spectra that were discussed in section 3. These serve as initial conditions
for the Meszaros equation including the Jeans length eq. (4.12) and are
evolved until today. For the non-linear regime we use the density power
spectrum without the oscillon contribution, as it was described in the
previous section, together with $d\delta_{k}/dx=0$ for the initial conditions.
The resulting power spectra are shown in figure 13 for
$m_{a}=10^{-18}\mathrm{eV}$, both $p=-1/2$ and $p=1/2$, and several values of
the misalignment angle. The black line corresponds to the CDM power spectrum
linearly evolved until today [65],

$P_{\mathrm{CDM}}(k)=2\pi\delta_{H}^{2}\frac{k^{n_{s}}}{H_{0}^{n_{s}+3}}T^{2}(k)\,.$
(4.19)

Here $\delta_{H}$ is a normalising constant we set to $1.9\times 10^{-5}$,
$n_{s}$ is the spectral index of primordial perturbations, which we set to 1
and $T(k)$ is the BBKS transfer function, given by [66]

$T(k)=\frac{\ln(1+0.171\kappa)}{0.171\kappa}\left(1+0.284\kappa+(1.18\kappa)^{2}+(0.399\kappa)^{3}+(0.490\kappa)^{4}\right)^{-1/4}\,,$
(4.20)

where $\kappa=k/k_{\mathrm{EQ}}$ and
$k_{\mathrm{EQ}}=0.073\Omega_{m,0}h^{2}\,\mathrm{Mpc}^{-1}$.

We then compute the variance of the filtered field $\sigma_{R}(a)$, where in
this work we always use the spherical top-hat window function. The halo mass
function today is then calculated using eq. (4.18). It is shown in figure 14,
where we plot the function $X_{M}\equiv(M/\rho_{\mathrm{DM}})dn/d\ln M=dF/d\ln
M$ for $m_{a}=10^{-18}\mathrm{eV}$, both $p=-1/2$ and $p=1/2$, and the same
values of the misalignment angle as in figure 13. The black line again
corresponds to the CDM halo mass function today.

Figure 14: The halo mass function $X_{M}=(M/\rho_{\mathrm{DM}})dn/d\ln M$ for
ALP DM with several values of the misalignment angle
$\theta_{i}=\phi_{i}/f_{a}$. We set $m_{a}=10^{-18}\mathrm{eV}$ using both
$p=-1/2$ (left) and $p=1/2$ (right). The mass function is enhanced on
intermediate masses compared to the CDM mass function, shown as the black
line, due to the enhancement of fluctuations. At small masses the halo mass
function is correspondingly suppressed compared to CDM.

###  4.3 Halo spectrum

While $\xi=0$ is a useful condition to identify the collapse time from eq.
(4.4), the condition itself is not physical. Instead of collapsing into a
singularity, the overdensity is expected to virialize into a halo, in which
the potential and the kinetic energies are settled to a stable equilibrium
[36, 61, 64, 67]. In this section we discuss the properties of such halos.

The radius of the virialized halo is half of the turnaround radius, $r_{\rm
vir}=r_{\rm ta}/2$, which follows directly from the virial theorem and energy
conservation. Hence the final virial density $\rho_{\rm vir}$ is eight times
the density at turnaround [68]. $\rho_{\rm vir}$ can also be expressed in
terms of the background energy density of ALPs. Remembering that the total
mass in the halo stays constant,

$M=\frac{4\pi}{3}\rho_{\mathrm{EQ}}a_{\mathrm{EQ}}^{3}(1+\delta_{\mathrm{in}})R^{3}=\frac{4\pi}{3}\rho_{\rm
ta}r_{\rm ta}^{3}=\frac{4\pi}{3}\rho_{\rm ta}a_{\rm ta}^{3}\xi^{3}_{\rm
ta}R^{3},$ (4.21)

the virial density can be expressed as

$\rho_{\mathrm{vir}}=8\rho_{\rm
ta}=8C_{\rho}\rho_{\mathrm{EQ}}(1+\delta_{\mathrm{in}})\delta_{\mathrm{in}}^{3},$
(4.22)

where we defined $C_{\rho}=1/(\xi_{\rm ta}^{3}C_{x}^{3})$ and used
$C_{x}\equiv x_{\rm ta}\delta_{\mathrm{in}}$ from above. Numerically we find
that $C_{\rho}\approx 17$. Inserting this into the above equation we arrive at
[61, 67]

$\rho_{\rm vir}\approx
136\rho_{\mathrm{EQ}}(1+\delta_{\mathrm{in}})\delta_{\mathrm{in}}^{3}\,.$
(4.23)

Using $\delta_{\mathrm{in}}\approx 1.1/x_{\rm coll}$ we can estimate the
virial density of the halo right after collapse as

$\rho_{\rm vir}(x_{\rm coll})\approx 136\overline{\rho}_{m,\rm
EQ}\quantity(1+\frac{1.1}{x_{\rm coll}})\frac{1.1^{3}}{x_{\rm
coll}^{3}}=136\overline{\rho}_{\rm coll}1.1^{3}\quantity(1+\frac{1.1}{x_{\rm
coll}}).$ (4.24)

For small fluctuations collapse happens during matter domination, so the
$x_{\rm coll}^{-1}$ term can be neglected. Then

$\rho_{\rm vir}(x_{\rm coll})\approx 136\overline{\rho}_{\rm
coll}1.1^{3}\approx 181\overline{\rho}_{\rm coll},\quad\delta_{\mathrm{in}}\ll
1,$ (4.25)

For collapse during radiation domination we obtain

$\rho_{\rm vir}(x_{\rm coll})\approx 136\overline{\rho}_{m,\rm
EQ}\frac{1.1^{4}}{x_{\rm coll}^{4}}\approx 200\overline{\rho}_{R,\rm
EQ}\frac{a_{\rm EQ}^{4}}{a_{\rm coll}^{4}}\approx 200\overline{\rho}_{\rm
coll},\quad\delta_{\mathrm{in}}\gg 1.$ (4.26)

The conclusion is that the virialized density of the halo right after collapse
is approximately $200$ times the mean density of the universe at collapse,
regardless of the cosmological era in which the collapse happens.

It has been confirmed, by means of $N$-body simulations, that CDM forms halos
with a characteristic, red-shift dependent profile, dubbed Navarro-Frenk-
White-(NFW)-profile. It is given by [69]

$\rho(r,z)=\frac{4\rho_{s}}{(r/r_{s})(1+r/r_{s})^{2}}\,,$ (4.27)

where $\rho_{s}\equiv\rho(r_{s})$ is the scale density, and $r_{s}$ is the
scale radius defined by $\left.d\ln\rho(r)/d\ln r\right|_{r_{s}}=-2$ that
determines how strongly the density is concentrated within the halo. The
numerical simulations show that $\rho_{s}$ depends on the red-shift at which
the halos forms, but it stays constant after the halo is formed. For the
remainder of the section we review an approach to find $\rho_{s}$ for a given
collapsing region.

An equivalent definition of the NFW-profile is

$\rho(r,z)=\frac{\delta_{\mathrm{char}}\rho_{c}(z)}{(r/r_{s})(1+r/r_{s})^{2}}\,,$
(4.28)

where $\rho_{c}$ is the critical density of the universe and
$\delta_{\mathrm{char}}$ is the characteristic density.

The mass inside a certain radius $r$ can be found by integrating eq. (4.27) up
to $r$, which gives [53]

$M(r)=4\pi\int_{0}^{r}dr^{\prime}\,r^{\prime
2}\rho(r^{\prime})=16\pi\rho_{s}r_{s}^{3}\tilde{f}(r/r_{s})\,,$ (4.29)

with

$\tilde{f}(x)=\ln(1+x)-\frac{x}{x+1}\,.$ (4.30)

The mass integrated up to $r_{s}$ is called the scale mass $M_{s}$. We
discussed that the density of a virialized halo is roughly 200 times the mean
density of the universe at the time of formation. This motivates the
definition of $r_{200}$, which is the radius at which the average density of
the enclosed halo is 200 times the critical density:

$200\rho_{c}=\left(\frac{4\pi}{3}r_{200}^{3}\right)^{-1}4\pi\int_{0}^{r_{200}}dr\,r^{2}\rho(r),$
(4.31)

such that at the time of formation, this corresponds roughly to the virial
density. The mass inside this region is defined to be

$M_{200}=\frac{4\pi}{3}(200\rho_{\mathrm{crit}})r_{200}^{3}.$ (4.32)

Since we assume that the mass of the halo is conserved during the collapse,
the spherical mass $M$ from eq. (4.3) can be identified with $M_{200}$, which
is constant by definition.

One now defines the concentration parameter $c_{200}$ by

$c_{200}\equiv\frac{r_{200}}{r_{s}}.$ (4.33)

We can use eq. (4.31) to find the characteristic density as a function of the
concentration parameter:

$\delta_{\mathrm{char}}=\frac{200}{3}\frac{c_{200}^{3}}{\ln(1+c_{200})-c_{200}/(1+c_{200})}.$
(4.34)

From the definitions, this equation of $\delta_{\mathrm{char}}$ leads to a
time independent $\rho_{s}$. We now use the formalism from refs. [70, 67] and
use these definitions to relate the quantities of the collapsing ALP regions
to the empirical results of the scale density of halos formed at different
redshifts.

Let us assume we identify a halo at redshift $z_{\mathrm{id}}$, which most of
the time will be today, $z_{\mathrm{id}}=0$. We now assign the formation
redshift $z_{\mathrm{form}}(M_{200},\mathfrak{f})$ which is defined as the
redshift at which half the mass of the halo was first contained in progenitors
more massive than some fraction $\mathfrak{f}<1$ of the mass $M_{200}$.666We
note that $z_{\mathrm{form}}$ is distinct from $z_{\mathrm{coll}}$ discussed
in the spherical collapse model. We can find this with the Press-Schechter
formalism by

$\mathrm{erfc}\left[X(z_{\mathrm{form}})-X(z_{\mathrm{id}})\right]=\frac{1}{2},$
(4.35)

with

$X(z)=\frac{\delta_{c}(z)}{\sqrt{2\left[\sigma^{2}(\mathfrak{f}M_{200},z)-\sigma^{2}(M_{200},z)\right]}},$
(4.36)

where $\delta_{c}(z)$ is the red-shift dependent critical density for
collapse, and $\sigma^{2}(M_{200},z)$ is the variance at redshift $z$
evaluated with a spherical tophat as window function for the smoothing scale

$R=\left(\frac{3M_{200}}{4\pi\rho_{m,0}}\right)^{1/3},$ (4.37)

assuming here that the ALPs make up all of dark matter and ignoring the
baryonic contribution to $\rho_{m}$.

The connection to the numerical simulations is now made by a redshift
dependent proportionality factor $\mathfrak{C}$ which relates the critical
density of the universe at formation to the scale density:

$\rho_{s}=\mathfrak{C}(z_{\mathrm{form}})\rho_{c}(z_{\mathrm{form}}).$ (4.38)

To summarise this section: We use the variances found in section 3 and plug
them into eq. (4.35), choose $\mathfrak{f}=10^{-2}$, to find
$z_{\mathrm{form}}$ for given smoothing scale, which sets the halo mass. We
use the fit function from the appendix C of ref. [37] to find $\mathfrak{C}$
for given $z_{\mathrm{form}}$. The relation
$\rho_{s}=\delta_{\mathrm{char}}(z)\rho_{c}(z)$ immediately gives the
characteristic density today, leading to the concentration parameter
$c_{200}$. Since we know $r_{200}$ from the identification of $M_{200}$ with
the mass of the halo, we find $r_{s}$ and $M_{s}$ from $c_{200}$. Scale mass
$M_{s}$ and $\rho_{s}$ are the two parameters we are looking for in the halo
spectrum and which we are going to use to find observational prospects in the
next chapter.

Solitons: The NFW profile fails to describe the halos formed by ALP dark
matter at small scales due to the wavelike nature of the ALP. It has been
shown numerically that the central profiles of ALP halos do not obey a “cuspy”
$r^{-1}$ profile [71], but form a (pseudo-) soliton, being the eigenstate of
the time-independent Schrödinger-Poisson equation. A fit for the soliton
profile is given by [72]:

$\rho_{\mathrm{sol}}(r)=\rho_{\mathrm{sol},0}\left[1+(2^{1/8}-1)\left(\frac{r}{r_{1/2}}\right)^{2}\right]^{-8},$
(4.39)

where

$\rho_{\mathrm{sol},0}\approx
1.9\times\left(\frac{10^{-17}\mathrm{\,eV}}{m_{a}}\right)^{2}\left(\frac{\mathrm{pc}}{r_{1/2}}\right)^{4}\,M_{\odot}\mathrm{pc}^{-3}$
(4.40)

is the core density and $r_{1/2}$ is the radius at which the density is half
of the core density. From this form of the profile, one can calculate, how the
scale mass and scale density of the solitons are related:

$\rho^{\mathrm{sol}}_{s}\approx 1.2\times
10^{-4}\left(\frac{m_{a}}{10^{-17}\mathrm{\,eV}}\right)^{6}\left(\frac{M_{s}^{\mathrm{sol}}}{M_{\odot}}\right)^{4}M_{\odot}\mathrm{pc}^{-3}.$
(4.41)

Far away from the center of the halo, the profile is expected to still follow
the NFW-profile. We will neglect a possible interpolation between the soliton
core at small scales and the NFW-profile further out. Therefore, we are going
to make the assumption that for a given scale mass $M_{s}$ the scale density
of solitons $\rho_{\mathrm{sol}}(m_{a},M_{s})$ is an upper bound on the scale
density and that denser halos cannot form for that mass. Halo spectra for
$p=1/2$ and $p=-1/2$ at a benchmark ALP mass and for different initial field
values can be seen in figure 15.

Figure 15: Halo spectra, showing the scale density $\rho_{s}$ as function of
the scale mass $M_{s}$, for axion-like particles with mass
$m_{a}=10^{-18}\mathrm{\,eV}$ and potential eq. (2.8) with $p=1/2$ (left) and
$p=-1/2$ (right). We compare different initial field values $\theta_{i}$. The
variance $\sigma_{R}$ has been calculated with the spherical window function.
The soliton line from eq. (4.41) for this mass is shown as straight grey line.
This line acts as an upper bound for the scale density rendering densities
above the soliton line as unphysical. This means the actual halo spectrum will
follow the line that corresponds to the minimum of soliton bound and the
calculated scale densities. We show the parts of the calculation of the scale
density that exceed the soliton line as transparent. The stars indicate where
the halo mass function, shown in figure 14, peaks. Also shown is the CDM-line
in black.

##  5 Observational prospects from gravitational interactions

Now that we have described the ALP halo spectra for a certain mass and given
initial conditions, we want to proceed by investigating in which part of the
parameter space in the $(m_{a},f_{a})$-plane gravitational interactions could
actually lead to observational signatures of the ALP halos.

When the halos form, they are much denser than their surrounding. This remains
true when they become part of larger dark matter halos. However, if, at later
times, they enter galaxies, much denser objects like stars will cause tidal
forces on the halos which could potentially disrupt them. Since part of the
observational signatures rely on the fact that the halos are inside of the
Milky Way or other galaxies, we need to make sure that they actually survive
until today. For halos smaller than the typical star separation,
$r_{s}<1\,\mathrm{pc}$, passing a star could actually disrupt the halo
completely, depending on its density. For larger halos, interaction with
single stars cannot disrupt the halo completely, but here one has to take into
account the tidal forces induced by the gravitational potential of the Galaxy.
To understand all of the halo population, in a full analysis one would need to
find the sub-halo mass function within the Galaxy, as has been done e.g. in
refs. [73, 74]. However, we are only interested if the densest halos survive.
Hence, without going into details we take the estimate from ref. [22] that a
halo smaller than $1\,\mathrm{pc}$ must have $\rho_{s}>10^{3}\rho_{\odot}$,
with $\rho_{\odot}\equiv 1.1\times 10^{-2}M_{\odot}\mathrm{pc}^{-3}$ being the
local dark matter density. Halos much larger than $1\,\mathrm{pc}$ only need
to be slightly denser than the local average, so we take
$\rho_{s}>11\rho_{\odot}$. Both of these estimates are rather conservative.
For $r_{s}\approx 1\,\mathrm{pc}$, we would actually need to interpolate
between the two regimes, but this goes beyond the scope of our work. We show
the parameter space in which halos inside the Galaxy likely did not survive
until today as light grey area in figure 16.

###  5.1 Signatures of ALP minihalos

There are several halo phenomena that could lead to observable effects. Halos
could act as gravitational lenses or they could affect the motion of stars
when they pass them on their movement through the Galactic halo. In the
following we discuss the different effects these gravitational interactions
could cause.

Local Gravitational Perturbations: One distinguishes three different cases of
a similar phenomenon: When dark matter subhalos traverse the Galaxy, they can
change the 6D phase space distribution of stars in the Milky Way disk [75] or
of stars in the stellar halo [76]. The third effect that is heavily
investigated is the disruption of stellar streams by dark matter subhalos [77,
78, 79]. These streams are the result of tidal disruption of globular clusters
(spherical collections of stars) in the Milky Way halo, resulting in narrow
leading and trailing arms of stars, since the mass loss in the clusters is
deposited into orbits with slightly higher and lower orbital energies [80,
81]. These stellar streams have been experimentally confirmed for example in
refs. [82, 83].

Ref. [75] gives as an approximate formula for the velocity kick a star gets
from the passing of the halo:

$\Delta\vec{v}\approx-\hat{b}\,0.5\,\mathrm{km~{}s}^{-1}\left(\frac{M(b)}{10^{7}M_{\odot}}\right)\left(\frac{166\,\mathrm{km~{}s}^{-1}}{V}\right)\left(\frac{\mathrm{kpc}}{b}\right)\,,$
(5.1)

where $V$ is the relative velocity between the subhalo and the star we will
assume to be $166\,\mathrm{km\,s}^{-1}$, $b$ is the impact distance, and
$M(b)$ is the mass of the subhalo that is enclosed within that impact
distance. Such a velocity kick could be observed by the space telescope GAIA,
whose on-going mission is set to make the most detailed and largest map of
stars in the Milky Way to date [84]. We follow ref. [22] and calculate the
parameter region in the $(M_{s},\rho_{s})$-plane where the halo can give stars
in the Galactic disk velocity kicks that are larger than
$2\,\mathrm{km~{}s}^{-1}$, which we present as dotted light green line in
figure 16, while the region where the stars could be ejected from the disk, or
even the Galaxy, i.e. $\Delta v\gtrsim 100\,\mathrm{km\,s}^{-1}$, is shown as
grey area in the upper right corner. To find these lines, we take the minimum
velocity kick $\Delta\vec{v}$ we want to achieve, insert the given mass of the
halo into eq. (5.1) as $M(b)$ and then find the corresponding $b$ that still
allows for the velocity kick we want to achieve. This velocity kick will only
be possible if all the mass of the halo actually is contained inside the
impact distance $b$, so we impose that the scale radius must fulfil
$r_{s}<b/3$ which for given mass automatically determines the minimum density
$\rho_{s}$ such that all of the halo lies withing the impact distance. We also
present the parameter space for which the halos could generate velocity kicks
large enough to disrupt a very cold stream, resulting in change in speed of
stars in the stream of $\Delta v\gtrsim 0.5\,\mathrm{km\,s}^{-1}$ for impact
distances and scale radii with $\max\\{b,r_{s}\\}\gtrsim 10\,\mathrm{pc}$.
This area is shown as dashed dark-green line in figure 16.

Figure 16: The thick lines show three of the halo spectra with the densest
characteristic halos, for $p=-1/2$, $\theta_{i}=4$ and a mass of
$m_{a}=10^{-18}\,\mathrm{eV}$ in red, and for a mass of
$m_{a}=10^{-10}\,\mathrm{eV}$ in green, and for $p=1/2$, for a mass of
$m_{a}=10^{-14}\,\mathrm{eV}$ and $\theta_{i}=10$ in blue. The CDM line is
given in black. The light grey area depicts the region in which halos in the
Galaxy likely have been tidally disrupted until today. We also show the
following parameter regions in which halos could be observable by their
gravitational interactions: Local Gravitational Perturbations: The dashed
green (dotted lime-green) line depicts the parameter space above which halos
could significantly change the velocity of stars in stellar streams (the
Galaxy’s disk), while the grey area shows where dark matter halos could lead
to ejection of stars from the disk or even the Galactic halo. The cyan area
shows the best-fit parameter space that could explain features in the GD-1
stream. Astrometric Weak Gravitational Lensing: The dotted purple line shows
the boundary between transient effects (to the left) and enduring effects (to
the right), where the following phenomena can only test the enduring regime.
Above the dotted (dashed) blue line the test statistic $\mathcal{T}_{\mu}$
predicts that local effects in the angular velocity $\mu$ of stars due to
lensing effects are observable by GAIA (SKA), while above the dotted (dashed)
red lines, corresponding to the global test statistic $\mathcal{C}_{\alpha}$,
correlations of the angular acceleration $\alpha$ due to enduring
gravitational lensing should be observable by GAIA (THEIA). Photometric
Lensing: The orange region shows where photometric microlensing, leading to
irregularities in the light curve of highly magnified stars that transit a
caustic curve, may be observable. Diffraction of Gravitational Waves: The
region in which we expect distortions in the gravitational wave signals from
BH-BH mergers caused by dark matter halos that could be observable by aLIGO
(LISA) [ET] is shown by the dark-orange line (black line) [purple region with
dashed-dotted boundary]. On the top we show the corresponding ALP mass $m_{a}$
that produces a halo spectrum peaking at given $M_{s}$, see eq. (5.16).

One promising stellar stream, which shows characteristic features (including a
gap and an off-stream spur of stars) that could be explained by an encounter
with a dense dark matter subhalo is the stream GD-1. Ref. [85] claims that
such an encounter is the best explanation, since encounters with known
globular clusters, dwarf galaxies or molecular clouds in the Milky Way disk
are unlikely to explain the observations. The densities of the dark matter
halos needed to explain the features are in a mass range of about
$10^{6}\,M_{\odot}$ to $10^{8}\,M_{\odot}$ and would need to be more dense
than the CDM halo spectrum predicts. However, the halos predicted by ALPs with
non-period potential could fall into the region that could explain the
features in the GD-1 stream. We copied this region from ref. [85] and show it
as cyan area in figure 16.

We note that these observational prospects do not assume any halo mass
function or value of energy density in halos of certain mass. For these
observations, encounters with one halo would suffice, while of course a small
number density would diminish the possibility to actually observe these
velocity kicks or make the encounter of a halo with the GD-1 stream less
likely. The results from the halo mass function (see section 4.2) imply that
the number density of dense halos is large enough to make the possibility of
these encounters significant.

Astrometric Weak Gravitational Lensing: If the compact halos cross the line-
of-sight between Earth and stars or other luminous objects, they can induce
apparent motion of the latter through gravitational lensing, without the need
for multiple images or magnification [22]. This effect, called astrometric
weak lensing, was first considered in ref. [86, 87, 88, 89] for point-like
objects and cuspy minihalos. We use the results of ref. [90] who proposed a
program of searches for compact subhalos using the astrometric possibilities
of the space telescope GAIA [84], working since 2013, of the Square Kilometer
Array (SKA) [91], a radio interferometer whose construction has started in
2022, and the proposed space telescope Telescope for Habitable Exoplanets and
Interstellar/Intergalactic Astronomy (THEIA) [92, 93].

We have to distinguish two categories of astrometric lensing signatures: If
the minimum impact parameter between the lensed object and the line-of-sight
towards the lens is smaller than the change in impact parameter within a
typical multi-year astrometric survey, we expect transient effects, while
enduring effects correspond to the case where the change of the impact
parameter is small during the observation time. This means in the transient
regime we would expect to actually see a change in a gravitational lensing
event, while in the enduring regime, the change of the lensing happens so slow
that the instantaneous angular deflection is in practice unobservable since
the survey does not record the true position of the sources. In ref. [90], the
local velocity dispersion of dark matter $v=166\,\mathrm{km\,s}^{-1}$ is
picked as a typical rate of change in the impact parameter. Together with an
observation time $\tau$ of roughly 5 years, this results in
$v\tau\sim\mathcal{O}(10^{3}\mathrm{pc})$ as the size of the impact parameter
that marks the boundary between the transient and enduring regime. Since
usually the deflection angle is maximised for impact parameters near the scale
radius, i.e. $b\sim r_{s}$, a subhalo can produce a gravitational lensing
transient only if [22]

$\rho_{s}\approx 2\times
10^{10}\rho^{\odot}_{\mathrm{DM}}\left(\frac{M_{s}}{M_{\odot}}\right)\left(\frac{10^{-3}\,\mathrm{pc}}{v\tau}\right)^{3}\,,$
(5.2)

where again $\rho^{\odot}_{\mathrm{DM}}\approx 1.1\cdot
10^{-2}\,M_{\odot}/\mathrm{pc}^{3}$ is the assumed local dark matter density
in the Solar System. The line that marks the boundary between these two
regions, with transient (enduring) effects possible for smaller (larger)
masses at fixed density is shown in purple in figure 16, where we took
$v\tau=2\times 10^{-3}$.

Contrary to what one would maybe expect at first sight, it seems more
promising to observe enduring effect than transient effects [22]. Since the
density of a halo, for given mass, has to be very large to produce a transient
event, while we do not expect extraordinarily large densities for smaller
halos, the expected angular deflection would be too small to be detectable for
current astrometric observatories. As already mentioned, in the enduring
regime the angular deflection is in principle unobservable. Nevertheless, ref.
[90] studied, how the apparent angular velocity $\mu$ and angular acceleration
$\alpha$ could hint to dark matter halos. We present two of their core results
regarding dark matter halos: Stellar location and velocity data can be tested
for consistence with the existence of a lens. This is done by assuming the
velocity of a star in the absence of a lens $\dot{\vec{\theta}}$ is a Gaussian
variable with known variance, while the local, multi-source test statistic
$\mathcal{T}_{\mu}$, which takes stellar position and velocity data as input,
is taken to be the sum over all stars in the dataset of the scalar product of
$\dot{\vec{\theta}}$ with $\vec{\mu}$, where $\vec{\mu}$ is a template for the
predicted angular velocity in presence of a dark matter halo acting as lens.
The predicted sensibility, i.e. the maximal signal to noise ratio, in ref.
[90] is

$\left<\max\mathrm{SNR}_{\mathcal{T}_{\mu}}\right>\approx
0.4\,f^{1/3}\left[\frac{M_{s}}{10^{7}M_{\odot}}\right]^{2/3}\left[\frac{10\,\mathrm{pc}}{r_{s}}\right]\left[\frac{N_{0}}{10^{7}}\right]^{1/2}\left[\frac{0.01}{\Delta\Omega}\right]^{1/6}\left[\frac{200\mu\mathrm{as\,y}^{-1}}{\sigma_{\mu,\mathrm{eff}}}\right],$
(5.3)

where, here as well as for the following observational prospects, it is
assumed that the halo mass function is a delta peak around $M_{s}$. This means
we can write

$X_{M}(M)=fM\delta(M-M_{s}),$ (5.4)

where $f$ is the fraction of all dark matter that is assumed to be in these
halos. This is called the monochromatic approximation.

$\Delta\Omega$ is the solid angle covered by the survey, $N_{0}$ is the
product of the angular number density of stars and $\Delta\Omega$, and
$\sigma_{\mu,\mathrm{eff}}$ is the average angular velocity noise, i.e. the
effective standard error for the velocities [90]. For the GAIA observations of
the Large and Small Magellanic Clouds we followed ref. [90] and took
$N_{0}=10^{7}$, $\Delta\Omega=0.01$ and
$\sigma_{\mu,\mathrm{eff}}=200\mu\mathrm{as\,y}^{-1}$, for quasar observations
with the SKA we took $N_{0}=10^{8}$, $\Delta\Omega=4\pi$ and
$\sigma_{\mu,\mathrm{eff}}=1\mu\mathrm{as\,y}^{-1}$. The projected curves for
halo parameters where the resulting signal-to-noise ratio is larger than
unity, i.e. $\left<\max\mathrm{SNR}\tau_{\mu}\right>>1$, for all densities
above the curve, where we took $f=0.3$, are shown in figure 16 in blue, where
the dotted line shows the result for the GAIA observations while the dashed
line shows the projected region for the SKA (or comparable radio
observatories). Eq. (5.3) tells us that the smallest observable density scales
as $\rho_{s}\propto f^{-1}$.

The global, multi-sourced test statistic $\mathcal{C}_{\alpha}$ measures
small-angle correlations in angular accelerations that gravitational lensing
induces on a field of stellar motions [90]. The results for this are shown in
figure 16 in red, where the dotted line is a projection for GAIA observations
of stars in the Galaxy, while the dashed line is a projection for THEIA
observations of stars. Since the calculation of the curves is rather involved,
we graphically copied the results given in ref. [90] and note that also here
the result scales with $\rho_{s}\propto f^{-1}$ [90, 22].

Diffraction of Gravitational Waves: If the gravitational wave signal of a
merger of two black holes [94] passes through or near a dark matter halo on
its way to Earth, the signature of the lensing may allow for another probe of
the halo mass distribution. This would be possible even if the lens is not
strong enough to produce multiple images of the same merger, which would then
be detected as temporally distinct events. Instead, the lens would imprint
characteristic distortions in the waveform as well as in the amplitude of the
gravitational wave signal [95].

The dimensionless parameter $\omega$ characterises the strength of the
distortions, and it is given as a function of frequency and mass enclosed in
the lens by [22]

$\omega\cong
1.3(1+z_{L})\left(\frac{f_{\mathrm{GW}}}{10^{2}\,\mathrm{Hz}}\right)\left(\frac{M_{\mathrm{enc}}}{100\,M_{\odot}}\right)\,,$
(5.5)

where $z_{L}$ is the redshift of the lens, $f_{\mathrm{GW}}$ is the frequency
of the gravitational wave, and $M_{\mathrm{enc}}$ is the enclosed mass within
the impact parameter of the lens. The equation is valid if the effective lens
distance is of $\mathcal{O}(\mathrm{Gpc})$ and if the effective velocity
dispersion is of $\mathcal{O}(\mathrm{km\,s}^{-1})$. Ref. [95] explains how
advanced LIGO (aLIGO) [96] will be able to probe black hole mergers with
$\omega\sim\mathcal{O}(1)$ out to distances of a Gigaparsec. Since aLIGO
operates in frequencies of $(10-10^{3})\,\mathrm{Hz}$, this makes it possible
to detect halos of masses $\sim\mathcal{O}(10-10^{3})\,M_{\odot}.$ The
probability that a BH-BH merger is passing behind a halo with impact parameter
$b$ is roughly given as [95, 22]:

$P(b)\sim
0.045f\left[\frac{1+z_{L}}{2}\right]^{3}\left[\frac{D_{\mathrm{BH}}}{5\mathrm{\,Gpc}}\right]\left[\frac{10^{5}M_{\mathrm{\odot}}}{M_{s}}\right]\left[\frac{b}{1\mathrm{\,pc}}\right]^{2},$
(5.6)

where $D_{\mathrm{BH}}$ is the distance to the merger. To find observable
regions, we followed ref. [22] and set $D_{\mathrm{BH}}=5\mathrm{\,Gpc}$ and
took the redshift to the lens to be $z_{L}=0.3$. We then found for every scale
mass the corresponding $b_{\mathrm{min}}$ at which $P(b)$ is $1\%$. If the
resulting $b_{\mathrm{min}}$ was smaller as the so-called NFW smoothing scale
$2\pi/(2m_{a}v_{s})$ with $v_{s}=\sqrt{16\pi\ln(2)G\rho_{s}r_{s}^{2}}$, then
this smoothing scale was taken as $b_{\mathrm{min}}$ instead. It was then
checked if there is a $b>b_{\mathrm{min}}$ such that the mass enclosed in a
cylinder within the impact parameter $b$ leads to a distortion in the range
$0.5<\omega<5$ for some frequency in the range of aLIGO, $(10-10^{3})\,$Hz,
the planned space interferometer LISA [97], $(10^{-4}-10^{-1})\,$Hz or the
proposed Einstein Telescope (ET) [98], $(1-10^{4})\,$Hz. The last condition
was that the $b$ we found was not larger as ten times the so-called Einstein
radius of the lens, which was taken to be
$r_{E}=0.1\mathrm{\,pc}\left[M_{\mathrm{enc}}/(100M_{\odot})\right]$ [95]. We
show the regions where this calculation predicts observable effects for aLIGO,
LISA and ET in figure 16 as dark-orange, black and purple dashed-dotted line,
respectively. We note that while this prediction depends on $f$, again using
the monochromatic approximation eq. (5.4), a decrease in $f$ makes denser
halos unobservable first, because while they induce strong lensing, the
probability that they cross our line of sight to a binary gets smaller, i.e.
the smaller $f$ in eq. (5.6) can be compensated by a larger $b$. The low
density cutoff for the observable region is instead given by the fact that for
halos with small density, $b$ will be larger than the Einstein radius.

Photometric Microlensing: Photometric lensing, gravitational lensing events
that change the apparent brightness of a celestial object, has been proposed
in ref. [99]. How such a phenomenon could be observable with the help of stars
in the Milky Way, the Magellanic Clouds or Andromeda in the case where the
lensing objects are MACHO’s, e.g. planetary size objects or primordial black
holes, has been investigated in refs. [100, 101, 102, 103, 104]. At small
sizes the halos we discuss here are not dense enough to have the same effects
as the dense objects discussed in the references above. Ref. [105] discusses,
how such halos could nevertheless be observed via photometric microlensing.
The general idea is the following: At specific distances behind a
gravitational lens there are locations of highest possible magnification. The
connection of all these locations is called caustic curve, while the
projection of such a caustic to the sky is called the critical curve.
Recently, observations of stars highly magnified to around a factor of
$\mathcal{O}(10^{3})$ by this mechanism found near a caustic curve have been
reported, e.g. in ref. [106]. In these cases the lenses are clusters of
galaxies and the highly magnified star is expected to be close to the critical
curve of the lensing cluster. Ref. [105] proposes to look for small variations
to the magnifying event to find the footprint of dark matter minihalos within
the lensing cluster. Single stars inside the cluster lens can disrupt the
cluster’s critical curve into a network of micro-critical curves (projections
of micro-caustics) on angular scales of $\sim 10\,\mu\mathrm{as}$, while dark
matter minihalos with their sub-planetary masses and solar system sizes are
expected to produce even finer magnification variations scales of $\sim
10^{2}-10^{3}\,\mathrm{nas}$ [105]. In such a microlensing event, when a
highly magnified star crosses a micro-caustic, the observed light curve would
be altered by surface density fluctuations caused by dark matter minihalos
inside the cluster. Ref. [105] predicts that such an event could be monitored
through observations with the most powerful optical and infrared telescopes
based in space or on Earth.

We follow ref. [105] and ref. [53] to find the parameters of halos that could
be observable by this effect. The density power spectrum of halos on non-
linear scales can be found from

$P_{\rho}=\overline{\rho}\int\frac{dM}{M^{2}}\frac{dF}{d\ln
M}\left|\tilde{\rho}^{h}(q;M)\right|^{2},$ (5.7)

where $\overline{\rho}$ is the averaged density in some region of the cluster
and $\tilde{\rho}^{h}(q;M)$ is the Fourier transform of the halo profile. One
now assumes that slices of the lensing cluster all have the same fractional
mass distribution and only vary by their averaged density, with which the
power spectrum of the surface density becomes an integral along the line of
sight:

$P_{\Sigma}(q_{\perp})=\overline{\Sigma}\int\frac{dM}{M^{2}}\frac{dF}{d\ln
M}\left|\tilde{\rho}^{h}(q_{\perp};M)\right|^{2},$ (5.8)

where $q_{\perp}$ is the surface Fourier mode and $\overline{\Sigma}=\int
dL\overline{\rho}(L)$ is the mean surface density. It is assumed that inside
the cluster it is dominated by the cluster’s surface density, i.e.
$\overline{\Sigma}=\overline{\Sigma}_{\mathrm{cl}}$.

It is convenient to define the lensing convergence

$\kappa=\frac{\Sigma}{\Sigma_{\mathrm{crit}}},$ (5.9)

with $\Sigma_{\mathrm{crit}}=c^{2}/(4\pi GD_{\mathrm{eff}})$ being the
critical surface density for an isolated lens, at which it produces multiple
images, where we find $D_{\mathrm{eff}}=D_{L}D_{LS}/D_{S}$ with $D_{L}$,
$D_{LS}$ and $D_{S}$, which are the distances to the lens, to the source and
from source to lens, respectively.

The dimensionless convergence power spectrum is now given by

$\Delta_{\kappa}^{2}=\frac{q^{2}P_{\Sigma}}{2\pi\Sigma_{\mathrm{crit}}^{2}}.$
(5.10)

In the monochromatic assumption eq. (5.4), and using a NFW profile, ref. [53]
gives an analytic solution for eq. (5.10):

$\Delta_{\kappa}=\frac{1}{\ln{2/\sqrt{e}}}\frac{\sqrt{\Sigma_{\mathrm{cl}}fM_{s}}qr_{s}g(qr_{s})}{\Sigma_{\mathrm{crit}}r_{s}},$
(5.11)

where $g(qr_{s})$ is a function appearing in the Fourier transform of the NFW
profile. $qr_{s}g(qr_{s})$ reaches a maximum of 0.35 at $qr_{s}=0.77$ [53], so
we will assume that for all $r_{s}$ this maximum can be reached at a
particular $q_{\mathrm{max}}$ by scanning through all $q$. We will now assume
that $D_{\mathrm{eff}}=1\,\mathrm{Gpc}$, set
$\Sigma_{\mathrm{cl}}=0.8\Sigma_{\mathrm{crit}}$ and introduce the distance
the star seems to traverse during the microlensing event $d$, which is given
by the product of relative velocity between cluster and star, the event
duration and the magnification. As a benchmark value, $d$ is set to
$10\,\mathrm{AU}$. One assumes that the following four conditions must be
fulfilled for possible observations of the halos [53]: Firstly, the maximum
value of the convergence power spectrum must fulfil
$\mathrm{max}\Delta_{\kappa}(q)>10^{-3}$, so we could observe a brightness
fluctuation of $\mathcal{O}(1)$, if the star of interest is magnified by a
factor of $10^{3}$. Secondly, to ensure statistical treatment, the distance
the star seems to traverse should sweep over many halos, i.e. $f\pi
d^{2}\Sigma_{\mathrm{cl}}/M_{s}>10$, which leads to a sharp cutoff at a
certain mass scale. The third criterion is to assume that the clumps should be
smaller than $d$, so the density fluctuations induced during the lensing event
are actually significant, i.e. $r_{s}<d$, and as fourth criterion we use that
$r_{s}>10\,\mathrm{AU}$, so that the halos are not so small such that
fluctuations in the light curve they induce are washed out by the finite size
of the lensed star.
We show the region at which this leads to detectable effects of such
photometric lensing events for $f=0.3$ in figure 16 as orange area.

Observing ALP halos with Pulsar Timing Arrays: Pulsar Timing Arrays (PTAs),
like NANOGrav [107], PPTA [108] or EPTA [109] recently became an exciting
possibility to test gravitational physics. One observes a set of pulsars,
which generally have an extremely stable rotation, making them very accurate
astrophysical clocks. Correlation between time residuals, the difference
between expected and measured arrival time of the signal, can be analysed to
find hints of astrophysical or cosmological gravitational waves, with the
possibility that a stochastic GW background has already been observed [110].

The passing of dark matter halos could be inferred from analysing the pulsar
timing data, with a halo/halos passing through the line of sight causing a
Shapiro delay, or halos accelerating Earth or the pulsar directly, causing a
Doppler effect. The analysis of PTA data for hints of DM halos relies on
searching for the expected signal with a filter, where the latter is only
easily found for events involving a single halo and a single pulsar (called
deterministic), while the stochastic analysis involving the whole halo
population and a large number of pulsars is less well understood. This leads
to the fact that studies like refs. [111, 112, 113, 73] give different results
regarding which halo populations could be observable. The observational
prospects depend e.g. on the number of observed pulsars, the absolute
observation time, how often every single pulsar is observed, the halo mass
function of the halos as well their concentration parameters, while the
predicted observable regions disagree additionally due to different filters
and statistical methods. A simplified analysis with a monochromatic approach
for the halo mass function does not suffice for the stochastic scenario, while
observation with current and near-future PTAs of single halos is unrealistic
for the model we investigate. Therefore we do not show any observational
prospects for PTAs and postpone a dedicated analysis of the observational
prospects for ALP halos to future work.

###  5.2 Parameter regions of observable gravitational effects

Figure 17: The peak values of the scale mass as a function of the initial
angle $\theta_{i}$ for $m_{a}=10^{-18}\,$\mathrm{eV}$$ obtained by calculating
the halo spectrum as described in section 4.3.

Our aim is to show in which region of the $(m_{a},f_{a})$-parameter space the
corresponding halos can cause observable effects.. Ideally, this requires
calculating the halo spectrum for each point in the $(m_{a},f_{a})$-parameter
space, which is not numerically feasible, especially since we need to run
lattice simulations for a sizeable region. Instead, we decided to extrapolate
our results for $m_{a}=10^{-18}\,$\mathrm{eV}$$ using a technique which we
will now explain.

We need a simple way to obtain $(M_{s}^{\rm peak},\rho_{s}^{\rm peak})$ as a
function of model parameters $(m_{a},f_{a},p)$. In figure 17 we show the
values of the scale mass at which the halo spectrum gets peaked. We only show
the benchmarks at which $\rho_{s}^{\rm peak}>10^{-1}M_{\odot}\textrm{pc}^{-3}$
since they are relevant for observations. We see that the halo spectra for
both $p=1/2$ and $p=-1/2$ are peaked at the scale mass $M_{s}\sim 2\times
10^{3}M_{\odot}$. To find the scaling of the scale mass with respect to the
ALP mass $m_{a}$ we can use the following relation for $M_{s}$

$M_{s}=\frac{f(1)}{f(c_{200})}\frac{4\pi}{3}\rho_{m,0}a_{0}^{3}R^{3},$ (5.12)

which can be obtained by combining eq. (4.29) and eq. (4.32). The dark matter
density today $\rho_{m,0}a_{0}^{3}$ is the same for all benchmark points. We
can also assume that $f(c_{200})$ at the peak of the halo spectrum is not
sensitive to the model parameters. Then we can write

$M_{s}^{\rm peak}\propto R_{\rm peak}^{3}\propto k_{\rm peak}^{-3},$ (5.13)

where $k_{\rm peak}$ is the comoving momentum around which structure is
enhanced. We expect that the enhancement occurs around the same dimensionless
momenta

$\tilde{k}_{\rm peak}=\frac{k_{\rm{peak}}/a}{\sqrt{2m_{a}H}}$ (5.14)

in all of the parameter space. Since $\tilde{k}$ is conserved during radiation
era, we have $\tilde{k}_{\rm peak}\propto k_{\rm peak}/\sqrt{m_{a}}$ which
results in

$M_{s}^{\rm peak}\propto k_{\rm peak}^{-3}\propto m_{a}^{-3/2}.$ (5.15)

So we can write the following relation for the scale mass:

$M_{s}^{\rm peak}\approx 2\times
10^{3}M_{\odot}\quantity(\frac{m_{a}}{10^{-18}\,$\mathrm{eV}$})^{-3/2}.$
(5.16)

A similar relation has also been found in ref. [22].

|  $p=1/2$ |  $p=-1/2$
---|---|---
$\alpha_{<}$ |  $-6.22$ |  $-7.18$
$\alpha_{>}$ |  $34.6$ |  $10.2$
$\beta_{<}$ |  $12.7$ |  $15.6$
$\beta_{>}$ |  $-26.8$ |  $-9.20$
$\gamma$ |  $8.02$ |  $4.67$
Table 1: Best fit parameters of the fit function eq. (5.17) for the peak scale
density. We fixed the ALP mass to be $m_{a}=10^{-18}\,$\mathrm{eV}$$. Figure
18: Comparison of the fit function for the peak scale density eq. (5.17)
(solid lines) with the values obtained by the full numerical calculation
(circles). The parameters of the fit function are tabulated in Table 1.

For the peak scale density, we found that the following fit function provides
a reasonable approximation:

$\ln\rho_{s}^{\rm
peak}(\theta_{i})=\quantity(\alpha_{<}+\beta_{<}\,\tilde{\theta}_{i}^{2})h\quantity(\gamma(1-\tilde{\theta}_{i}))+\quantity(\alpha_{>}+\beta_{>}\,\sqrt{\tilde{\theta}_{i}})h\quantity(\gamma(\tilde{\theta}_{i}-1)),$
(5.17)

where

$h(x)=\frac{1}{2}+\frac{1}{2}\tanh{x}\mbox{\quad
and\quad}\tilde{\theta}_{i}\equiv\frac{\theta_{i}}{\theta_{i}^{\rm peak}}.$
(5.18)

Here $\theta_{i}^{\rm peak}$ is the value of the initial angle for which the
halo spectrum has the strongest peak. We found $\theta_{i}^{\rm peak}\approx
10$ and $\theta_{i}^{\rm peak}\approx 4$ for $p=1/2$ and $p=-1/2$,
respectively. The best-fit values for the remaining parameters are given in
table 1. A comparison between the numerically obtained peak values and this
fit function is given in figure 18.

Figure 19: Regions in the $(m_{a},f_{a})$-plane in which the ALP dark matter
can lead to observable gravitational effects. The corresponding initial field
values $\theta_{i}$ are shown by the grey contour lines. These projections
assume that 30% of all of dark matter is inside of characteristic halos with
masses within one decade around $M_{s}$. The upper plot shows $p=-1/2$, while
the lower plot shows $p=1/2$. The labelled regions show the constraint on the
minimum mass from the Lyman-$\alpha$ forest (for $p=1$), the minimum mass at
which oscillons decay before matter-radiation equality, the constraint from
heating of stellar orbits and the constraint from (astrophysical) black hole
superradiance. Figure 20: Same as figure 19, but using that 1% of all dark
matter is inside characteristic halos with masses within one decade around
$M_{s}$.

We can now combine the last two subsections by starting with a given tuple of
$m_{a}$ and $f_{a}$. When we assume that
$\Omega_{a,0}=\Omega_{\mathrm{DM},0}$, the initial tuple ($m_{a},f_{a}$) fixes
the initial misalignment $\theta_{i}=\phi_{i}/f_{a}$ to achieve the correct
relic abundance by eq. (2.10) and with this pair of $m_{a}$ and $\theta_{i}$,
we can find the approximate peak mass and peak density of the characteristic
halo via eq. (5.16) and eq. (5.17). The parameter space in which we expect
observable gravitational effects, as they have been discussed in the previous
subsections and presented in figure 16, from these halos can be seen in figure
19 for $X_{M}=0.3M\delta(M-M_{s})$ and in figure 20 for
$X_{M}=0.01M\delta(M-M_{s})$. Since the axion-like-particle we investigate
could be only gravitationally coupled to the SM, we consider the few
constrains that depend only on their self-interactions and coupling to
gravity. We added a line to show the currently strongest constraint from the
Lyman-$\alpha$ forest at $2\times 10^{-20}\,\mathrm{eV}$ [114]777See refs.
[115, 116, 117, 118] for earlier works on constraining the fuzzy dark matter
mass via Lyman-$\alpha$ observations., for which we point out that it has been
derived for ALPs with harmonic potential and could be modified for other
potentials [119], but still prospects for smaller masses should be taken with
a grain of salt. We also show the constraints from superradiance around
astrophysical black holes. In this process, dense clouds of fundamental
scalars accumulate around black holes, extracting the latter’s rotational
energies. For this process to be efficient and the clouds to be long-lived,
the scalar’s Compton wavelength has to be roughly the size of the black hole
horizon. Therefore, observing fast spinning black holes at certain masses
excludes the corresponding mass range of the axion-like-particle. Because
strong self-interaction quench the cloud, superradiance can only constrain
large $f_{a}$. We show the excluded region found in ref. [120]. We also show
the mass at which oscillons decay before matter-radiation equality, eq.
(3.12), meaning that for lower axion masses the presence of oscillons could
influence the evolution of the halos, which we did not include in our
analysis. In a recent work, ref. [121] claims that a fuzzy dark matter mass
smaller than $3\times 10^{-19}\,\mathrm{eV}$ is in conflict with the observed
stellar velocity dispersion inside the ultra-faint dwarf galaxies Segue 1 and
Segue 2. The main idea is that the wave interference effects produce large
fluctuations in the local density and gravitational potential inside the fuzzy
dark matter halos. These fluctuations cause gravitational heating of stellar
systems [122] and increase the velocity dispersion of the stars. We indicate
this bound in our plots with a green band with the label "stellar heating".

As can be seen, the ALP halos typically lead to observable effects in a region
of the parameter space where $\theta_{i}$ is in the range $1$-$100$. For these
values of the field, also the change in the equation of state due to the
transfer of energy into fluctuations is of $\mathcal{O}(1)$, as it was studied
in ref. [30] for a similar potential. This justifies neglecting this
$\mathcal{O}(1)$-correction in this work.

##  6 Comparison between axion fragmentation models

In this section, we compare our results with the halo predictions from the
Large Misalignment Mechanism (LMM) [22] and the Kinetic Misalignment Mechanism
(KMM) [26, 25] analyzed in refs. [24, 37].

We already discussed the LMM in section 2.1.1. Now we shall give a brief
review of ALP dark matter in KMM. See ref. [24] for more details. In KMM, the
ALP field is assumed to have an initial kinetic energy parameterised by the
"yield" parameter defined by

$Y=\frac{f_{a}\dot{\phi}(T)}{s(T)},$ (6.1)

where $\dot{\phi}(T)$ and $s(T)$ are the ALP velocity and the entropy of the
universe at temperature $T$, respectively. At early times, the ALP field is
dominated by its kinetic energy so its energy density scales as
$\rho_{\phi}\propto a^{-6}$. This implies that $\dot{\phi}\propto a^{-3}$, so
the yield parameter is constant if the entropy is conserved. The large initial
kinetic energy delays the onset of oscillations from its conventional value
$m_{a}\sim H_{\rm osc}$. Instead, the oscillations start when the kinetic
energy drops below the potential energy, i.e.

$\frac{1}{2}\dot{\phi}^{2}(T_{\ast})=2f_{a}^{2}m_{a}^{2},$ (6.2)

where $T_{\ast}$ is defined by this equation and is called the _trapping
temperature_. For a given $m_{a}$ and $f_{a}$, the parameters $Y$ and
$T_{\ast}$ can be fixed such that the ALP field gives the desired amount of
dark matter. In particular, the yield is fixed by

$h^{2}\Omega_{\phi,0}\approx h^{2}\Omega_{\rm
DM,0}\,\quantity(\frac{m_{a}}{5\times
10^{-3}\,$\mathrm{eV}$})\quantity(\frac{Y}{40}),$ (6.3)

where the trapping temperature is given by the relation

$\frac{T_{\ast}}{\Lambda_{b}}=(2\times
10^{8})^{1/3}\quantity(\frac{g_{s}(T_{\ast})}{72})^{-1/3}\quantity(\frac{\Lambda_{b}}{$\mathrm{GeV}$})^{1/3}\quantity(\frac{h^{2}\Omega_{\phi,0}}{h^{2}\Omega_{\rm
DM}})^{-1/3},$ (6.4)

where $\Lambda_{b}=\sqrt{m_{a}f_{a}}$. The delay of the onset of oscillations
is quantified by $m/H_{\ast}$ where $H_{\ast}=H(T_{\ast})$. The Kinetic
Misalignment only occurs when this quantity is larger than the conventional
value, i.e. when $m/H_{\ast}\gtrsim 3$. The field is completely fragmented if
$m_{\ast}/H_{\ast}\gtrsim 40$.

The $a^{-6}$ scaling of the ALP energy density in Kinetic Misalignment can be
problematic if the trapping happens very late. In this case, the ALP field
will have too much energy during Big Bang Nucleosynthesis (BBN) so that the
constraints from observations of primordial helium-4 and deuterium abundances
are violated. This constraint is given by

$T_{\ast}\gtrsim 20\,$\mathrm{keV}$\mbox{\quad or\quad}\Lambda_{b}\gtrsim
9\times 10^{-7}\times\quantity(\frac{h^{2}\Omega_{\phi,0}}{h^{2}\Omega_{\rm
DM}})^{1/2}.$ (6.5)

Thus, in contrast to an ALP with non-periodic potential or the Large
Misalignment Mechanism, the Kinetic Misalignment is incompatible with small
ALP masses such as $m_{a}\lesssim 10^{-16}\,$\mathrm{eV}$$.

Model | Line style |  $\mathbf{\it f_{a}}$ [GeV] | Model specific
---|---|---|---
Non-periodic Potential $p=+1/2$ (blue) | solid | $9\times 10^{13}$ | $\theta_{i}=10$
dashed | $1\times 10^{14}$ | $\theta_{i}=8$
dotted | $7\times 10^{13}$ | $\theta_{i}=12$
Non-Periodic Potential $p=-1/2$ (purple) | solid | $9\times 10^{13}$ | $\theta_{i}=5$
dashed | $1\times 10^{14}$ | $\theta_{i}=4$
dotted | $7\times 10^{13}$ | $\theta_{i}=6$
Kinetic Misalignment (red) | solid | $9\times 10^{13}$ | $m/H_{\ast}\approx 20$
dashed | $1\times 10^{14}$ | $m/H_{\ast}\approx 12$
dotted | $7\times 10^{13}$ | $m/H_{\ast}\approx 27$
Large Misalignment (green) | solid | $9\times 10^{13}$ | $\absolutevalue{\pi-\theta_{i}}\approx 1\times 10^{-4}$
dashed | $1\times 10^{14}$ | $\absolutevalue{\pi-\theta_{i}}\approx 6\times 10^{-3}$
dotted | $7\times 10^{13}$ | $\absolutevalue{\pi-\theta_{i}}\approx 2\times 10^{-6}$
Post-inflationary
(yellow) | solid | $2\times 10^{14}$ |
Table 2: The tabulated list of parameters for the benchmark models in figure
21. Here $\theta_{i}=\phi_{i}/f_{a}$ denotes the initial value of the ALP
field, and $H_{\ast}$ is the Hubble at the time when the ALP field gets
trapped in the Kinetic Misalignment Mechanism.

For the comparison, shown in figure 21, we fix the ALP mass to
$m_{a}=10^{-14}\,$\mathrm{eV}$$ and choose several different values for the
decay constant $f_{a}$. Assuming that the ALP makes up all of dark matter,
this fixes the initial field value $\phi_{i}$ in the misalignment model with
the non-periodic potential and the LMM, and the yield and the trapping
temperature in the KMM. All the model parameters are also tabulated in table
2.

In the post-inflationary scenario, the initial angle $\theta_{i}$ is not a
free parameter anymore. Therefore, for a given ALP mass $m_{a}$, there is
single value for the ALP decay constant $f_{a}$ such that ALPs make all of
dark matter. The relic abundance for this case has recently been calculated in
ref. [123] as888This result includes contribution from the decay of
topological defects such as strings and domain walls. However, due to the
computational limitations they were able to simulate up to $\kappa\approx 8$
where $\kappa\equiv\ln(m_{s}/H)$ is the string tension, and $m_{s}$ is the
mass of the radial mode. In realistic scenarios, $kappa$ is much larger, for
example $\kappa\sim 70$ for the QCD axion. How the extrapolation should be
done is still an open problem. See refs. [124, 125, 126] for more details.

$h^{2}\Omega_{a}\simeq
0.019\quantity(\frac{g_{\rho}(T_{1})}{70})^{3/4}\quantity(\frac{g_{s}(T_{1})}{70})^{-1}\quantity(\frac{m_{a}}{$\mathrm{\SIUnitSymbolMicro
eV}$})^{1/2}\quantity(\frac{f_{a}}{10^{12}\,$\mathrm{GeV}$})^{2},$ (6.6)

where $T_{1}$ is the temperature such that $m_{a}=(8/5)H(T_{1})$. From this