aainstitutetext: Deutsches Elektronen-Synchrotron DESY, Notkestr. 85, 22607
Hamburg, Germanybbinstitutetext: II. Institute of Theoretical Physics,
Universität Hamburg D-22761, Germany

# ALP Dark Matter Mini-Clusters from Kinetic Fragmentation

Cem Eröncel a,b Géraldine Servant<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

We show that very compact axion mini-clusters can form in models where axion-
like-particle (ALP) dark matter is produced via the kinetic misalignment
mechanism, which is well-motivated in pre-inflationary $U(1)$ symmetry
breaking scenarios. This is due to ALP fragmentation. We predict denser halos
than what has been obtained so far in the literature from standard
misalignment in post-inflationary $U(1)$ breaking scenarios or from large
misalignment. The main reason is that adiabatic fluctuations are significant
at early times, therefore, even if amplification from parametric resonance
effects is moderate, the final size of ALP fluctuations is larger in kinetic
misalignment. We compare halo mass functions and halo spectra obtained in
kinetic misalignment, large misalignment and standard misalignment
respectively. Our analysis does not depend on the specific model realization
of the kinetic misalignment mechanism. We present our results generally as a
function of the ALP mass and the ALP decay constant only. We show that a
sizable region of this ALP parameter space can be tested by future experiments
that probe small-scale structures.

DESY-22-115

##  1 Introduction

Axion-like-particles (ALPs) are prime candidates for dark matter, and they
have been subject to extensive theoretical research in recent years (see
Marsh:2015xka ; Irastorza:2018dyq ; Grin:2019mub ; Choi:2020rgn for recent
reviews). Typically, they are pseudo-Nambu-Goldstone bosons of a global $U(1)$
symmetry that is spontaneously broken at some high scale $f$ which is referred
as the ALP decay constant. In this picture, ALPs are massless after symmetry
breaking due to the continuous shift symmetry. They obtain a mass at a lower
scale $\Lambda_{b}\ll f$ via their interactions with a strong gauge sector
that breaks the continuous shift symmetry explicitly. This explicit breaking
creates a periodic potential
$V(\theta,T)=\Lambda_{b}^{4}(T)\quantity[1-\cos(\theta)]$ for the ALP field
$\theta$ which leads to an ALP mass $m(T)\propto\Lambda_{b}^{2}(T)/f$. By far
the best-known ALP model is the QCD axion, see DiLuzio:2020wdo for a recent
review, that arises as a by-product of the Peccei-Quinn solution to the Strong
CP problem Peccei:1977hh ; Weinberg:1977ma ; Wilczek:1977pj . In this case,
the strong sector is QCD so that $\Lambda_{b}\sim\Lambda_{\rm QCD}$, thus the
QCD axion mass $m_{0}\equiv m(T=0)$ and $f$ are related to each other. In a
generic ALP model, both $m_{0}$ and $f$ are free parameters.

The couplings of ALPs to Standard Model particles depend on the UV completion
of a particular ALP model, however they typically scale as $f^{-1}$. The weak
coupling and the small mass make ALPs cosmologically stable, hence they are
well-motivated dark matter candidates Preskill:1982cy ; Dine:1982ah ;
Abbott:1982af . The most widely studied production mechanism goes by the name
_misalignment mechanism_. After the spontaneous breaking of $U(1)$, the ALP
field selects a random initial angle $\theta_{i}=[-\pi,\pi)$ on each Hubble
patch. The later evolution is governed by
$\ddot{\theta}+3H(t)\dot{\theta}+m^{2}(t)\sin\theta=0$ which is derived from
the Klein-Gordon equation with the Friedmann-Lemaitre-Robertson-Walker (FLRW)
metric by neglecting the spatial gradient terms. Here $H(t)$ is the Hubble
scale, and $m(t)$ is the axion mass which is in general temperature-dependent.
The initial velocity is usually assumed to be vanishing. In this case, the ALP
field is initially frozen due to the Hubble friction. Later it starts to
oscillate around the minimum $\theta=0$ once the Hubble scale drops below the
ALP mass. After this point, the ALP field behaves effectively as pressureless
cold matter, and contributes to the dark matter relic density.

Recently, some works modified the usual assumption that the ALP field is
static in the early universe, and studied the consequences of a large initial
kinetic energy for the ALP field Co:2019jts ; Chang:2019tvx . The kinetic
energy causes a delay in the onset of oscillations so that the ALP dark matter
parameter space can be expanded to lower values of the ALP decay constant,
offering very promising prospects for upcoming experiments, see in particular
the recent overview in Eroncel:2022vjg . This variation is referred as the
_kinetic misalignment mechanism (KMM)_ , whereas the original case is called
the _standard misalignment mechanism (SMM)_. Such a kinetic energy can arise
if the $U(1)$ symmetry is explicitly broken at a very high scale by the high-
dimensional operators which induces a kick for the angular (axion) mode once
the radial mode (saxion) starts oscillating Co:2019jts ; Co:2019wyp ;
Co:2020dya ; Co:2020jtv ; Co:2021lkc ; Gouttenoire:2021jhk . Another
possibility, named _trapped misalignment_ , is that the ALP potential has a
non-trivial temperature dependence such that the ALP field first starts
oscillating around the high-temperature minimum $\theta=\pm\pi$, then move to
the low-temperature minimum $\theta=0$ after some critical temperature $T_{c}$
DiLuzio:2021gos ; DiLuzio:2021pxd . A large kinetic energy can also be
realized in models with more than one ALP fields Daido:2015bva ; Daido:2015cba
.

In Eroncel:2022vjg , we pointed out that ALP fluctuations play a crucial role
in the KMM setup. When the ALP field is rolling over the potential barriers
and when it is oscillating around the minimum, the fluctuations can grow
exponentially via parametric resonance. This effect is called axion
fragmentation Fonseca:2019ypl . It was demonstrated in Eroncel:2022vjg that in
large part of the KMM parameter space all the energy in the homogeneous mode
is converted into ALP fluctuations. This means that ALP dark matter consists
of particles that are mildly relativistic at the end of fragmentation, and
cool down later via redshift. This is in contrast with the SMM case where the
dark matter is an extremely cold coherently oscillating scalar field. In this
work, we continue our study of ALP fluctuations in the KMM setup, and discuss
the observational consequences of fragmentation for the ALP miniclusters. We
calculate the halo spectrum, typical densities of the ALP halos for a given
halo mass, and show that a sizable region of the $(m,f)$-ALP parameter space
can be tested by future experiments that probe small-scale structures.

Most of the literature on ALP miniclusters assumes the scenario where the
$U(1)$ symmetry breaking, Peccei-Quinn (PQ) in the case of the QCD axion,
takes place after inflation, referred as the _post-inflationary scenario_. In
this case, the ALP field takes uncorrelated initial values in causally
disconnected Hubble patches. Since our universe contains many of these Hubble
patches, the ALP field does have $\mathcal{O}(1)$ density fluctuations in the
post-inflationary scenario. These fluctuations gravitationally collapse very
early, and form dense objects called ALP miniclusters Hogan:1988mp ;
Kolb:1993zz ; Kolb:1993hw ; Kolb:1994fi . In recent years, this process has
been studied in detail both semi-analytically Hardy:2016mns ; Enander:2017ogx
; Fairbairn:2017dmf ; Fairbairn:2017sil , and numerically Bertschinger:1985pd
; Zurek:2006sy ; Ricotti:2009bs ; PhysRevD.96.123519 ; Delos:2017thv ;
Vaquero:2018tib ; Buschmann:2019icd ; Eggemeier:2019khm ; OHare:2021zrq . For
a recent review see Niemeyer:2019aqm . A recent discussion of how the
properties of the QCD axion miniclusters in the post-inflationary scenario are
modified in the presence of kinetic misalignment can be found in
Barman:2021rdr .

In the case when $U(1)$ breaking happens before or during inflation, all
Hubble patches are inflated to huge sizes so that the initial angle is the
same in all the observable universe. This is referred as the _pre-inflationary
scenario_. Even in this scenario the ALP field has some fluctuations.
Unavoidable ones are the adiabatic fluctuations that arise due to the
temperature fluctuations of the radiation bath. If ALPs are present during
inflation, and they are massless or light compared to the inflation scale
$H_{I}$, they also pick up quantum fluctuations given by $\delta\theta\simeq
H_{I}/(2\pi f_{I})$ where $f_{I}$ is the effective decay constant during
inflation which might be different from its value today Linde:1991km . The
ratio of the isocurvature fluctuations to the adiabatic ones are highly
constrained by the Planck 2018 data Aghanim:2018eyx ; Planck:2018jri .
Throughout this work we will consider the pre-inflationary scenario, and
assume that the inflationary isocurvature fluctuations are negligible.

Even though the initial size of the fluctuations is small in the pre-
inflationary scenario, in some cases fluctuations can be enhanced
significantly due to various instabilities Johnson:2008se . A well-known
example related to the periodic potential is when the initial angle
$\theta_{i}$ is very close to the top of the potential111Such an initial
condition can be realized with dynamical mechanisms without tuning Co:2018mho
; Takahashi:2019pqf ; Huang:2020etx ., i.e. $\absolutevalue{\pi-\theta_{i}}\ll
1$, referred as the _large misalignment mechanism (LMM)_ Arvanitaki:2019rax or
_extreme axion_ Zhang:2017dpp . In this case, the onset of oscillations is
delayed from the conventional condition $m\sim 3H$ due to the tiny potential
gradient at the top. As we will see in detail in Section 3.3, this delay
causes the fluctuations to grow exponentially via tachyonic instability and
parametric resonance PhysRevD.50.7690 ; Greene:1998pb ; Zhang:2017flu ;
Zhang:2017dpp ; Cedeno:2017sou ; Arvanitaki:2019rax . This feature also arises
in ALP models with non-periodic potentials for sufficiently large initial
angles Jaeckel:2016qjp ; Alonso-Alvarez:2018tus ; Olle:2019kbo ;
Berges:2019dgr ; Chatrchyan:2020pzh ; Chatrchyan:2022askdja . The common
conditions in all these cases are the delay of the onset of oscillations and a
large initial field value. The latter causes the ALP field to probe the non-
quadratic parts of its potential which causes instabilities, while the former
ensures that the field amplitude decays slowly so that the fragmentation
becomes efficient. We expect that the effects studied in this work are
relevant for any ALP model which predict a delay in the onset of oscillations,
such as the recently proposed _frictional misalignment_ Papageorgiou:2022prc .

The enhancement of fluctuations does also occur in KMM Eroncel:2022vjg ;
Fonseca:2019ypl . The main reason is the delayed onset of oscillations due to
the initial kinetic energy. As a result, the physics responsible for
fragmentation in KMM and LMM models is quite similar. However, there are
subtle differences which we study thoughtfully in Section 3.3. We demonstrate
that in the LMM scenario, some modes exhibit significant growth even before
the homogeneous mode starts oscillating since they become tachyonic. Due to
this effect, the relative amplification is stronger in LMM compared to KMM.
However, the initial size of the fluctuations is much smaller in LMM due to
the fact that the homogeneous mode behaves as dark energy $w\approx 0$ in the
early time limit, and dark energy has no adiabatic perturbations Marsh:2015xka
. On the other hand, the equation of state of the homogeneous mode in KMM is
$w\approx 1$, thus adiabatic perturbations are not zero at early times. We
shall see that because of this difference in the initial conditions, the final
size of the fluctuations are larger in KMM compared to LMM, even though the
amount of amplification is stronger in the latter.

The fragmentation process modifies the matter power spectrum significantly at
the scales $k\sim m_{\rm osc}a_{\rm osc}$ that get amplified most strongly.
Here $m_{\rm osc}$ and $a_{\rm osc}$ are the ALP mass and the scale factor at
the onset of oscillations. Therefore, fragmentation hypothesis can be tested
via experiments that probe the matter power spectrum Brito:2022lmd . An
important consequence of fragmentation is that the peak of the matter power
spectrum can reach $\mathcal{O}(1)$ values right after the fragmentation.
These fluctuations gravitationally collapse very early in the matter era, and
can create very dense and compact ALP miniclusters like in the post-
inflationary scenario Arvanitaki:2019rax . Our primary goal in this work is to
calculate the properties of these miniclusters semi-analytically, and comment
on their discovery prospects. We compare the predictions of all the production
mechanisms that we mentioned so far: SMM, LMM, KMM, post-inflationary
scenario. Such a comparison can be the first step towards an ambitious goal of
inferring the ALP production mechanism from the small-scale power spectrum
observations. Our analysis is model-independent. The link to specific model
implementations of KMM will be presented in Eroncel:2022a .

In Eroncel:2022vjg we demonstrated that fragmentation occurs both in
temperature-dependent and independent potentials. In this work, we concentrate
on ALP models with temperature-independent potentials. The constant ALP mass
is denoted by $m$. A discussion of the temperature-dependent potential
including a precise analysis for the QCD axion will be presented in a future
publication.

ptBenchmark $(m,f)$$T_{\ast}$Initial conditionsEq. (3.17) for KMMEq. (3.40)
for SMM/LMMSolve for the mode functionsEqs. (3.10) and (3.11)Density
contrastEq. (3.4)Power spectrumEq. (3.22)Use the ansatz (4.13)for the mode
functions Eroncel:2022vjg Density contrastEq. (4.7)Power spectrumEqs.
(4.9)(4.11)(4.12)(4.24)Late time evolution (well after fragmentation) Sec.
3.2, Fig. 4, Eq. (3.27)Gravitational collapseSec. 5.1, Eq. (5.13)Press-
Schechter formalism Sec. 5.2Halo mass function Sec. 5.3, Fig. 18Halo
spectrumSec. 5.4, Fig. 20Eq. (2.5)ptIncompletefragmentation$m/H_{\ast}\lesssim
40$ptCompletefragmentation$m/H_{\ast}\gtrsim 40$ptEvolution during incomplete
fragmentationptEvolution during complete fragmentation Figure 1: Summary of
the procedure followed in this paper.

The main programme of the paper is sketched in Figure 1. The ultimate goal is
to predict the observable ALP DM small-scale structures that can form at late
times from gravitational collapse as signatures of the kinetic misalignment
production mechanism. In Section 2, we briefly review the key results of
Eroncel:2022vjg which are relevant for this work. The most important quantity
for this paper is the trapping temperature $T_{\ast}$ that can be calculated
for any benchmark point $(m,f)$. From this, one can derive the initial
conditions for the ALP fluctuations. This is exposed in detail in Appendix A.
The main quantity that is needed to determine the late-time evolution of the
fluctuations is their power spectrum. There are two approaches for this
calculation one can follow depending on the relative values of the ALP mass
$m$ and the Hubble rate at trapping $H_{*}$ that determine whether
fragmentation is complete or incomplete. In Section 3, we calculate the power
spectrum corresponding to the ALP density contrast in KMM with fragmentation
using the cosmological perturbation theory. We also do the same exercise for
LMM, and compare the two mechanisms. In Section 4, we describe a semi-
analytical estimate of the power spectrum in the case of complete
fragmentation where the cosmological perturbation theory breaks down. In this
section, we also make a comparison with the post-inflationary scenario. Next,
In Section 5, the formation of the dark matter halos is studied analytically
via the Press-Schechter (PS) formalism. We use the results obtained for the
power spectrum to derive the halo mass function (HMF), and the halo spectrum.
We briefly describe the observational prospects in Section 6, and finally
conclude in Section 7. More technical details are presented in appendices.
Important equations are inside frames. Those equations which are new are in
addition in blue background.

##### Notation:

We use the metric convention $\operatorname{diag}\quantity(-,+,+,+)$, and
$M_{Pl}\approx 2.435\times 10^{18}$ GeV denotes the reduced Planck mass.

##  2 Brief review of ALP fragmentation in kinetic misalignment

In this section, we summarize the main findings of Eroncel:2022vjg . We
consider an ALP field $\theta$ with the Lagrangian

$\mathcal{L}=-\frac{f^{2}}{2}g^{\mu\nu}\partial_{\mu}\theta\partial_{\nu}\theta-V(\theta)=-\frac{f^{2}}{2}g^{\mu\nu}\partial_{\mu}\theta\partial_{\nu}\theta-m^{2}f^{2}\quantity(1-\cos\theta).$
(2.1)

In KMM, the evolution of the ALP homogeneous mode $\Theta$ can be divided into
two regimes. When the ALP kinetic energy $f^{2}\dot{\Theta}^{2}/2$ dominates
over the size of the potential barrier $2m^{2}f^{2}$, then the ALP field rolls
with the velocity $\dot{\Theta}\propto a^{-3}$ so that its energy density
redshifts as $\rho_{\Theta}\propto a^{-6}$. Shortly after its kinetic energy
becomes sub-dominant compared to barrier height, it starts to oscillate around
one of the minima, and behaves as Cold dark matter (CDM).

For the early evolution, one can introduce the quantity called the _yield_
defined by

$Y=f^{2}\dot{\Theta}(T)/s(T)$ (2.2)

where $s(T)$ is the entropy density of the universe. Assuming entropy
conservation, this quantity is conserved at early times when
$\rho_{\Theta}\propto a^{-6}$. The ALP relic density today can be expressed as

$h^{2}\Omega_{\Theta,0}\approx h^{2}\Omega_{\rm DM}\quantity(\frac{m}{5\times
10^{-3}\,$\mathrm{eV}$})\quantity(\frac{Y}{40}),$ (2.3)

where we took $\Omega_{\rm DM}=0.12$ Aghanim:2018eyx . This result is also
valid for a temperature-dependent potential with $m$ replaced by $m_{0}$.

The ALP field gets _trapped_ once its kinetic energy becomes comparable to the
barrier height. This defines the _trapping temperature_ $T_{\ast}$:

$\frac{1}{2}f^{2}\dot{\Theta}^{2}(T_{\ast})=2f^{2}m^{2}(T_{\ast}).$ (2.4)

This depends on the barrier height at zero temperature, and on the high-
temperature scaling of the potential. For a temperature-independent potential
it simplifies to Eroncel:2022vjg

$\boxed{\frac{T_{\ast}}{\Lambda_{b}}\approx\quantity(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_{\Theta,0}}{h^{2}\Omega_{\rm
DM}})^{-1/3}}\ ,$ (2.5)

where $g_{s}$ denotes the number of effective degrees of freedom in the
entropy.

The different fragmentation regimes in KMM depend strongly on the hierarchy
between the ALP mass and Hubble scale at trapping, i.e. $m/H_{\ast}$, which
can be obtained after solving for $T_{\ast}$ in (2.5). These regimes can be
summarized as follows:

  1. 1.

Standard misalignment ($m/H_{\ast}<3$): The onset of oscillations is not
delayed from its conventional value $m(T_{\rm osc})=3H(T_{\rm osc})$, and
there is no kinetic misalignment.

  2. 2.

Kinetic misalignment with incomplete fragmentation ($3<m/H_{\ast}\lesssim
4\times 10^{1}$):  The onset of oscillations is delayed due to the non-zero
initial velocity, but the fragmentation is weak such that the energy density
in fluctuations is always sub-dominant.

  3. 3.

Complete fragmentation after trapping ($4\times 10^{1}\lesssim
m/H_{\ast}\lesssim 5\times 10^{2}$): The ALP field is completely fragmented,
and the fragmentation ends after it would have been trapped by the potential
in the absence of fragmentation, i.e. $T_{\rm end}<T_{\ast}$.

  4. 4.

Complete fragmentation before trapping ($5\times 10^{2}\lesssim
m/H_{\ast}\lesssim 5\times 10^{3}$): The ALP field is completely fragmented,
and the fragmentation ends before it would have been trapped by the potential
in the absence of fragmentation, i.e. $T_{\rm end}>T_{\ast}$.

  5. 5.

Non-perturbative ($5\times 10^{3}\lesssim m/H_{\ast}$): The ALP field is non-
perturbative already at the beginning of fragmentation due to the large
adiabatic perturbations, so our analytical solution breaks down.

In this work, we will denote the second region as the _linear regime_ ,
corresponding to incomplete fragmentation, while the third and fourth regions
will be called the _non-linear regime_ , corresponding to complete
fragmentation.

Our primary goal in this paper is to study the observational consequences of
fragmentation. For this, we need to calculate the power spectrum of the ALP
density perturbations at late times. Our pipeline to calculate the power
spectrum depends whether the corresponding benchmark point is in the linear or
non-linear regime.

  * •

For the benchmark points in the linear regime, we will use the cosmological
perturbation theory to solve the mode functions numerically without relying on
the analytical theory that we have developed in Eroncel:2022vjg . We do this
beacuse in this regime $m/H_{\ast}$ is not very large so the analytical method
of Eroncel:2022vjg is less reliable. Also, the exact evolution can be obtained
without too much difficulty by solving a couple of differential equations
numerically. The procedure is described in detail in Section 3.

  * •

For the benchmark points in the non-linear regime, the use of the cosmological
perturbation theory is not possible. Even though precise results will require
a non-perturbative analysis such as a lattice simulation, we will rely on the
analytical theory in Eroncel:2022vjg to calculate the late-time behavior of
the mode functions. We then calculate the power spectrum directly from these
mode functions using the method in Enander:2017ogx . This is the topic of
Section 4.

##  3 ALP density perturbations in the incomplete fragmentation regime

In this section, we calculate the power spectrum corresponding to the ALP
energy density at late times. This is the necessary ingredient in order to
study the halo spectrum and the gravitational signatures that we will work out
in the next sections. In Section 3.1 we will calculate the evolution of the
mode functions during the fragmentation from which we will obtain the
evolution of the density contrast. In Section 3.2 we describe an effective
description of the density contrast at late times. Finally, in Section 3.3 we
will compare the power spectra of KMM and LMM models at late times and discuss
the similarities and differences.

###  3.1 Early evolution of fluctuations: before and during fragmentation

Our starting point is the Friedmann-Lemaitre-Robertson-Walker metric including
curvature perturbations $\Phi,\Psi$ in the Newtonian gauge222We adopt the
conventions used in maggiore2008gravitational .:

$\differential{s}^{2}=a^{2}(\eta)\quantity{-\quantity[1+2\Psi(\eta,\mathbf{x})]\differential{\eta}^{2}+\quantity[1+2\Phi(\eta,\mathbf{x})]\delta_{ij}\differential{x}^{i}\differential{x}^{j}},$
(3.1)

The anisotropic stress vanishes for scalar fields at the linear order in
perturbation theory which sets $\Psi=-\Phi$. At early times, we assume that
the dominant component of the energy density is radiation, so $\Phi$ is an
external field for the ALP.

We expand the ALP field as

$\theta(\eta,\mathbf{x})=\Theta(\eta)+\delta\theta(\eta,\mathbf{x}),$ (3.2)

and keep only the terms that are linear in fluctuations $\delta\theta$ and
$\Phi$. Since the fluctuations are small, the average energy density is given
by the energy density of the homogeneous mode:

$\expectationvalue{\rho}=\rho_{\Theta}=\frac{f^{2}}{2}\quantity(\frac{\Theta^{\prime}}{a})^{2}+f^{2}m^{2}\quantity(1-\cos\Theta),$
(3.3)

where $\prime$ denotes derivative with respect to conformal time. The density
contrast is

$\delta\equiv\frac{\delta\rho_{\theta}}{\rho_{\Theta}}=\frac{f^{2}}{\rho_{\Theta}}\quantity[\frac{\Theta^{\prime}\delta\theta^{\prime}}{a^{2}}+\Phi\frac{\Theta^{\prime
2}}{a^{2}}+m^{2}\sin\Theta\delta\theta].$ (3.4)

The homogeneous mode satisfies the following equation of motion in physical
time $t$:

$\ddot{\Theta}+3H\dot{\Theta}+m^{2}\sin\Theta=0,$ (3.5)

while the equation of motion for the mode functions $\theta_{k}$ reads

$\ddot{\theta}_{k}+3H\dot{\theta}_{k}+\quantity(\frac{k^{2}}{a^{2}}+m^{2}\cos\Theta)\theta_{k}=2\Phi_{k}m^{2}\sin\Theta-4\dot{\Phi}_{k}\dot{\Theta},$
(3.6)

where $\Phi_{k}$ are the Fourier modes of the curvature perturbations.
Following Arvanitaki:2019rax , we introduce the following dimensionless
quantities:

$t_{m}\equiv\frac{m}{2H}\simeq mt,\quad
t_{k}\equiv\frac{k/a}{\sqrt{3}H},\quad\widetilde{k}^{2}\equiv\frac{k^{2}/a^{2}}{2mH},$
(3.7)

where the approximate expression for $t_{m}$ is valid during radiation
domination. The momentum variable $\widetilde{k}$ is constant during radiation
era, therefore it can be used to uniquely identify the momentum mode. The
variable $t_{k}$ is convenient to express the exact solution for the curvature
perturbations $\Phi_{k}$ in radiation era:

$\Phi_{k}(t_{k})=3\Phi_{k}(0)\quantity(\frac{\sin t_{k}-t_{k}\cos
t_{k}}{t_{k}^{3}}),$ (3.8)

where $\Phi_{k}(0)$ is a stochastic variable that depends on the amplitude of
the scalar primordial fluctuations $A_{s}$ by

$\expectationvalue{\absolutevalue{\Phi_{k}(0)}^{2}}=\quantity(\frac{2}{3})^{2}\expectationvalue{\absolutevalue{\mathcal{R}_{k}(0)}^{2}}=\quantity(\frac{2}{3})^{2}\quantity(\frac{2\pi^{2}}{k^{3}})A_{s}\quantity(\frac{k}{k_{\star}})^{n_{s}-1}.$
(3.9)

Here $\mathcal{R}_{k}$ is the comoving curvature perturbation,
$k_{\star}=0.05\,\textrm{Mpc}^{-1}$ is the pivot scale, and $n_{s}$ is the
spectral tilt. In this work we take $n_{s}=1$ for simplicity, and set
$A_{s}=2.1\times 10^{-9}$ consistent with the Planck 2018 measurements
(TT,TE,EE+lowE+lensing 68%) Aghanim:2018eyx .

In terms of the variables defined in (3.7), the equations of motion for the
homogeneous mode (3.5) and the mode functions (3.6) during the radiation era
take the form Arvanitaki:2019rax ; Zhang:2017dpp ; Zhang:2017flu 333Our
convention for the curvature perturbations differs by a sign compared to
Arvanitaki:2019rax .:

$\displaystyle\derivative[2]{\Theta}{t_{m}}+\frac{3}{2t_{m}}\derivative{\Theta}{t_{m}}+\sin\Theta$
$\displaystyle=0,$ (3.10)
$\displaystyle\derivative[2]{\theta_{k}}{t_{m}}+\frac{3}{2t_{m}}\derivative{\theta_{k}}{t_{m}}+\quantity(\frac{\widetilde{k}^{2}}{t_{m}}+\cos\Theta)\theta_{k}$
$\displaystyle=2\quantity(\Phi_{k}\sin\Theta-\frac{t_{k}}{t_{m}}\derivative{\Phi_{k}}{t_{k}}\derivative{\Theta}{t_{m}}),$
(3.11)

where we have used $H=(2t)^{-1}$. After choosing the correct initial
conditions for the homogeneous mode $\Theta$ and the mode functions
$\theta_{k}$, the evolution of the ALP field and its fluctuations can be
determined by numerically solving (3.10) and (3.11).

####  3.1.1 Initial conditions

To fix the initial conditions for the homogeneous mode, we used the $\varrho$
parameter that we have introduced in Eroncel:2022vjg :

$\varrho(t_{m})\equiv\frac{\rho_{\Theta}(t_{m})}{2m^{2}f^{2}}=\frac{1}{4}\quantity(\derivative{\Theta}{t_{m}})^{2}+\sin^{2}\frac{\Theta}{2}.$
(3.12)

The homogeneous mode is rolling when $\varrho>1$, and oscillating when
$\varrho<1$. In the former case, $\varrho$ satisfies the following relation
Eroncel:2022vjg :

$1=\tau^{3/2}\sqrt{\varrho(\tau)}\,\textrm{E}\quantity(\frac{1}{\sqrt{\varrho(\tau)}}),\quad\varrho>1$
(3.13)

where E is the complete elliptic integral of the second kind444We use the
following definition for the elliptic integral:
$\textrm{E}(k)=\int_{0}^{\pi/2}\differential{\varphi}\sqrt{1-k^{2}\sin^{2}\varphi}$
Note that most software packages such as Mathematica and scipy 2020SciPy-NMeth
use $m=k^{2}$ instead of $k$ as their argument when defining elliptic
integrals. , and

$\tau\equiv 2H_{\ast}t=(2H_{\ast}/m)t_{m}$ (3.14)

is an another dimensionless time variable defined such that $\tau=1$ coincides
with the time of trapping. At any $\tau$, this equation can be solved
numerically to get $\varrho(\tau)$. At early times when $\varrho\gg 1$, we can
approximate

$\varrho(\tau\ll
1)\approx\frac{4}{\pi^{2}}\frac{1}{\tau^{3}}=\frac{4}{\pi^{2}}\quantity(\frac{m}{2H_{\ast}})^{3}\frac{1}{t_{m}^{3}}.$
(3.15)

By plugging this result into (3.12) we get

$\evaluated{\derivative{\Theta}{t_{m}}}_{t_{m,i}}=2\sqrt{\frac{4}{\pi^{2}}\quantity(\frac{m}{2H_{\ast}})^{3}\frac{1}{t_{m}^{3}}-\sin^{2}\frac{\Theta_{i}}{2}}.$
(3.16)

Here $\Theta_{i}$ is the initial angle at the start of the simulation which is
irrelevant in the case of kinetic misalignment mechanism555For some choices of
$t_{m,i}$ and $\Theta_{i}$ the homogeneous mode can accidently stop on the top
of the potential. In this case, all the fragmentation effects we discuss in
this section will be amplified. We will neglect such accidental configurations
in this work.. We have chosen it as $\Theta_{i}=\pi/2$.

In order to fix the initial conditions for the mode functions, we solve (3.10)
and (3.11) analytically at early times when the ALP mass can be neglected by
assuming that all the modes are sourced from adiabatic perturbations and were
super-horizon when ALP homogeneous mode started to scale as kineation, i.e.
$\rho_{\Theta}\propto a^{-6}$. We describe the calculation in Appendix A.1,
and in much more detail in Eroncel:2022vjg . For $t_{k}<1$, we can express the
solution as a power series in $t_{k}$ as (see (A.7))
$\displaystyle\text{KMM:}\quad\theta_{k}(t_{m})$ $\displaystyle\approx
2t_{m}\derivative{\Theta}{t_{m}}\Phi_{k}(0)\sum_{i=0}^{4}c_{i}t_{k}^{2i},$
(3.17) where the coefficients $\quantity{c_{i}}$ are

$c_{0}=-\frac{1}{2},\quad c_{1}=\frac{23}{20},\quad
c_{2}=-\frac{491}{1680},\quad c_{3}=\frac{4427}{\numprint{151200}},\quad
c_{4}=-\frac{\numprint{146131}}{\numprint{93139200}}.$ (3.18)

The conversion between $t_{k}$ and $t_{m}$ can be done via

$t_{k}=\sqrt{\frac{4t_{m}}{3}}\widetilde{k}.$ (3.19)

The approximation (3.17) is quite accurate up to $t_{k}=1$. Therefore we have
chosen our initial time early enough so that for all simulated modes $t_{k}<1$
initially.

Figure 2: Evolution of the density contrast, obtained as the numerical
solution of (3.10) and (3.11) together with the initial conditions (3.17). The
amplitude is normalized with respect to the value of the primordial curvature
perturbations $\absolutevalue{\Phi_{k}(0)}=(2/3)\sqrt{A_{s}}\approx 3\times
10^{-5}$. The corresponding benchmark point, $m=10^{-15}\,$\mathrm{eV}$$ and
$f=1.3\times 10^{14}\,$\mathrm{GeV}$$, lies in the region where fragmentation
is not complete. The dashed lines show the CDM predictions which are
calculated using (3.20).

####  3.1.2 Density contrast evolution

The evolution of the density contrast is shown for a couple of modes in Figure
2, for the benchmark point $m=10^{-15}\,$\mathrm{eV}$$ and $f=1.3\times
10^{14}\,$\mathrm{GeV}$$, which lies in the second region (incomplete
fragmentation) described in Section 2, but is close to the boundary with the
third region. The dashed black vertical line shows the time at which the
homogeneous mode gets trapped $t_{m,\ast}=m/2H_{\ast}$. For comparison we also
show the CDM predictions with the colored dashed lines. These can be
calculated analytically in the radiation era, and they read Zhang:2017flu

$\delta_{\rm{CDM},k}=9\Phi_{k}(0)\quantity[\frac{\sin t_{k}}{t_{k}}+\frac{\cos
t_{k}}{t_{k}^{2}}-\frac{\sin t_{k}}{t_{k}^{3}}-\operatorname{Ci}(t_{k})+\ln
t_{k}+\gamma_{\rm E}-\frac{1}{2}],$ (3.20)

where $\operatorname{Ci}$ is the cosine integral, and $\gamma_{\rm E}$ is the
Euler-Mascheroni constant. Note that all modes approach to $3\Phi_{k}(0)$ in
the super-horizon limit. This is becuase in this limit the adiabatic
perturbations satisfy

$\frac{\delta_{i}}{1+w_{i}}=\frac{\delta_{j}}{1+w_{j}},$ (3.21)

for all species $i$ and $j$, where $\delta$ and $w$ are the density contrast
and the equation of state respectively. During the radiation era the radiation
perturbations are related to the curvature perturbations by $\delta_{r}\approx
2\Phi_{k}(0)$ in the super-horizon limit. Plugging this relation into (3.21)
and using the fact that the ALP equation of state is $w_{\theta}\approx 1$
before trapping gives $\delta\approx 3\Phi_{k}(0)$ at early times in the
super-horizon limit.

Later behavior can be interpreted as follows:

  * •

$\widetilde{k}=0.1$: This mode enters the horizon at
$t_{m}=3/4\widetilde{k}^{2}=75$ long after the homogeneous mode starts
oscillating and redshifts as matter. Therefore, this mode locks on the CDM
evolution once the homogeneous mode gets trapped.

  * •

$\widetilde{k}=5.0$: This mode experiences some growth slightly before the
homogeneous mode gets trapped, but after that it oscillates with a constant
amplitude which is smaller compared to the CDM prediction.

  * •

$\widetilde{k}=1.7$: This mode experiences significant growth due to the
parametric resonance effect explained in Eroncel:2022vjg , and its amplitude
is enhanced significantly compared to CDM.

Figure 3: The dimensionless power spectrum for different benchmarks points
(colored lines), and the CDM power spectrum (black line) calculated via (3.20)
at $t_{m}=m/2H=2\times 10^{3}$. For all benchmarks points, the ALP mass is
$m=10^{-15}\,$\mathrm{eV}$$. We observe that all power spectra converge to the
CDM prediction for low momenta. By decreasing the ALP decay constant $f$, we
increase the hierarchy between the ALP mass and the Hubble scale at trapping,
$m/H_{*}$, which implies more efficient fragmentation.

####  3.1.3 Power spectrum

By solving for the mode functions, we obtain the dimensionless power spectrum

$\mathcal{P}_{\delta}(k)=\frac{k^{3}}{2\pi^{2}}P_{\delta}(k)=A_{s}\quantity(\frac{2}{3})^{2}\widetilde{\delta}_{k}^{2},$
(3.22)

where $\delta_{k}=\widetilde{\delta}_{k}\Phi_{k}(0)$. We show it in Figure 3
for different $f$ values, and compare to the CDM prediction calculated via
(3.20). We set the ALP mass to $m=10^{-15}\,$\mathrm{eV}$$. Smaller ALP decay
constants correspond to larger hierarchies between the ALP mass and the Hubble
scale at trapping which implies more efficient fragmentation. This feature can
easily be observed in the plot. We can also see that all power spectra
converge to the CDM line at low momenta since these modes enter the horizon
long after the homogeneous mode gets trapped and redshifts as matter.

###  3.2 Evolution at late times: Well after fragmentation until today

In order to study the observational signatures such as the halo spectrum, we
need to evolve the power spectrum until today. Performing this by numerically
solving the mode functions is not feasible due to the rapid oscillations and
very long time scales. However, an effective description can be obtained for
sub-horizon modes by using the fact that long after the onset of oscillations
$a\gg a_{\ast}$, i.e. $t_{m}\gg 1$, the ALP energy density redshifts like
matter with an average equation of state $w\approx 0$. Therefore, one can use
the WKB approximation for both the homogeneous mode $\Theta$, and the mode
functions $\theta_{k}(t)$. One makes the ansatz Park_2012

$\displaystyle\Theta(t)$
$\displaystyle=a^{-3/2}\quantity[\Theta_{+}\cos(mt)+\Theta_{-}\sin(mt)],$
(3.23) $\displaystyle\theta_{k}(t)$
$\displaystyle=\theta_{+}(k,t)\cos(mt)+\theta_{-}(k,t)\sin(mt),$ (3.24)

where $\Theta_{\pm}$ are constants. With this approximation, one can directly
derive the evolution of the density constrast as Park_2012

$\ddot{\delta}_{k}+2H\dot{\delta}_{k}+\quantity[c_{s,\rm
eff}^{2}\quantity(\frac{k^{2}}{a^{2}}+16\pi G\rho_{r})-4\pi
G\rho_{\Theta}]\delta_{k}=8\pi G\rho_{r}\delta_{r,k},$ (3.25)

where $\rho_{r}$ is the energy density in the radiation, $\delta_{r}$ is the
density constrast of the radiation, and $c_{s,\rm eff}^{2}$ is the effective
sound speed of the ALP:

$c_{s,\rm
eff}^{2}=\frac{1}{4}\frac{k^{2}}{a^{2}m^{2}}\quantity(1+\frac{1}{4}\frac{k^{2}}{a^{2}m^{2}})^{-1}\approx\frac{1}{4}\frac{k^{2}}{a^{2}m^{2}},$
(3.26)

where the approximation holds for non-relativistic modes. Even though this
expression is derived in the axion comoving gauge, additional terms that arise
when converting this result into any other gauge decay on sub-horizon scales
Hlozek:2014lca . The source term on the RHS of (3.25) is proportional to the
density constrast of radiation which oscillates rapidly for sub-horizon modes.
Therefore, by averaging over these oscillations we can neglect the source term
for the sub-horizon modes. The $16\pi G\rho_{r}$ term on the LHS is negligible
during the matter era since $c_{s,\rm eff}^{2}\rho_{r}\ll\rho_{\Theta}$666Here
we assume that ALPs make all of dark matter so $\rho_{\Theta}=\rho_{\rm DM}$.
. Also, in radiation era we have $16\pi G\rho_{r}\approx 6H^{2}\ll
k^{2}/a^{2}$ for sub-horizon modes. In the end, we find that in both the
radiation and matter era, the evolution of the density constrast on sub-
horizon scales at late times can be described by

$\boxed{\ddot{\delta}_{k}+2H\dot{\delta}_{k}+\quantity(c_{s,\rm
eff}^{2}\frac{k^{2}}{a^{2}}-4\pi G\rho_{\Theta})\delta_{k}=0}$ (3.27)

The terms $4\pi G\rho_{\Theta}$ is responsible for the graviational collapse
and the structure formation during the matter era, while $c_{s,\rm
eff}^{2}k^{2}/a^{2}$ is the “pressure” term. The CDM evolution can be obtained
by setting the sound speed to zero. On large scales $c_{s,\rm
eff}^{2}k^{2}/a^{2}\ll 4\pi G\rho_{\Theta}$, density wins over the pressure so
the density constrast grows like in the CDM case; logarithmically during the
radiation era, and linearly during the matter era with respect to the scale
factor. However, at small scales the pressure wins, and the density contrast
oscillates with a constant amplitude. The wavenumber at which both become
equal is known as the _axion Jeans scale_ given by Marsh:2015xka

$k_{J}=\quantity(16\pi Ga\rho_{\Theta,0})^{1/4}\sqrt{m}=66.5\times
a^{1/4}\quantity(\frac{h^{2}\Omega_{\Theta}}{h^{2}\Omega_{\rm
DM}})^{1/4}\sqrt{\frac{m}{10^{-22}\,$\mathrm{eV}$}}\,\rm{Mpc}^{-1}.$ (3.28)

The take-away lesson is for ALPs there is scale-dependent growth, and ALP dark
matter differs from CDM on scales below the axion Jeans scale.

Figure 4: Evolution of the density contrast (solid lines) for three different
momentum modes from the time they enter the horizon until today. Up to the
vertical dotted line on the left, the evolution is obtained by solving the
full equations of motion (3.10) and (3.11), while the WKB approximation (3.27)
is used for the remaining evolution. The colored dashed lines show the CDM
evolution which we have obtained by matching the exact solution (3.20) in the
ratiation era with the WKB approximation with $c_{s}^{2}=0$. The lowest
momentum mode shown in blue locks on to the CDM evolution when the homogeneous
mode gets trapped, and follows it until today. The intermediate mode shown in
orange gets enhanced significantly at trapping, then oscillates with a
frequency given by the sound speed until it drops below the Jeans scale in the
matter era. After that it starts to grow linearly. The highest momentum mode
shown in green does not experience exponential growth at trapping, and
oscillates rapidly in its remaining evolution. The grey shaded region on the
top right denotes where the overdensities collapse gravitationally. Its
expression is derived in Section 5.1.

Following Arvanitaki:2019rax , we define

$y\equiv a/a_{\rm eq}.$ (3.29)

In terms of this variable (3.25) becomes777In this equation $\widetilde{k}$ is
constant and defined in radiation era via (3.7).

$y(1+y)\derivative[2]{\delta_{k}}{y}+\quantity(1+\frac{3}{2}y)\derivative{\delta_{k}}{y}+\quantity(\frac{\widetilde{k}^{4}}{y}-\frac{3}{2})\delta_{k}=0.$
(3.30)

Therefore, the momentum mode $\widetilde{k}$ drops below the Jeans scale at
$y=2\widetilde{k}^{4}/3$. We match the evolution of the density constrast
obtained by solving the mode function equation of motion (3.11) with the
differential equation given above at some $t_{m}^{\rm match}$ which is much
later than the onset of oscillations $t_{m}^{\rm match}\gg 1$, but long before
the matter-radiation equality $y_{\rm match}\ll 1$. By using the fact that
$y=2^{1/4}\sqrt{H_{\rm eq}t_{m}/m}$ in radiation era, we find the matching
conditions as

$\evaluated{\delta}_{t_{m}^{\rm match}}=\evaluated{\delta}_{y_{\rm
match}}\mbox{\quad
and\quad}\evaluated{2t_{m}\derivative{\delta}{t_{m}}}_{t_{m}^{\rm
match}}=\evaluated{y\derivative{\delta}{y}}_{y_{\rm match}}.$ (3.31)

The full evolution of the density contrast of three momentum modes
$\widetilde{k}=0.1$, $\widetilde{k}=1.73$, $\widetilde{k}=5.0$ is plotted in
Figure 4 for a benchmark point $m=10^{-15}\,$\mathrm{eV}$$ and $f=1.3\times
10^{14}\,$\mathrm{GeV}$$. The full solution of the equations of motion (3.10)
and (3.11) is matched with the solution of the WKB approximation (3.27) by
using (3.31) at $t_{m}=2\times 10^{3}$ or $y\sim 10^{-8}$. The momentum mode
$\widetilde{k}=0.1$ shown in blue continues to track the CDM evolution shown
via colored dashed lines. The mode $\widetilde{k}=1.73$ shown in orange
oscillate with a constant amplitude while it is above the Jeans scale. Once it
drops below the Jeans scale during the matter era, it starts to grow linearly.
The evolution of the mode $\widetilde{k}=5.0$ is similar, except it does not
experience exponential growth during trapping, and drops below the Jeans scale
much later.

###  3.3 Comparison with the Large Misalignment Mechanism

In this section, we compute the power spectrum of the ALP fluctuations in the
Large Misalignment Mechanism (LMM), and make the comparison with the Kinetic
Misalignment. This mechanism Zhang:2017dpp ; Arvanitaki:2019rax is a special
case of the Standard Misalignment Mechanism where the initial angle
$\Theta_{i}$ is very close to the top of the cosine potential, i.e.
$\absolutevalue{\pi-\Theta_{i}}\ll 1$.

The calculation is the same as in the KMM case with the exception of initial
conditions for the homogeneous mode, and the mode functions. We start by
discussing the homogeneous mode evolution. This will also give us the
expression for the relic density today, which we will use to determine the
model parameters. Then we state the initial conditions for the mode functions,
and finally show the comparison.

#### Evolution of the homogeneous mode

We again use the dimensionless time $t_{m}=m/2H\simeq mt$ as the dynamical
variable. The initial conditions are given by

$\lim_{t_{m}\rightarrow 0}\Theta(t_{m})=\Theta_{i},\quad\lim_{t_{m}\rightarrow
0}\Theta^{\prime}(t_{m})=0.$ (3.32)

If the initial angle is very close to the bottom of the potential
$\absolutevalue{\Theta_{i}}\ll 1$, then the cosine potential can be
approximated by a quadratic, which implies $\sin\Theta\approx\Theta$. In this
case the homogeneous mode equation of motion (3.10) with initial conditions
(3.32) has an exact solution:

$\Theta(t_{m})=2^{1/4}\Gamma\quantity(\frac{5}{4})\Theta_{i}\frac{J_{1/4}(t_{m})}{t_{m}^{1/4}},\quad\Theta_{i}\ll
1.$ (3.33)

In this case, the energy density is given by

$\rho_{\Theta}(t_{m})=\frac{m^{2}f^{2}}{2}\quantity[\Theta^{\prime}(t_{m})^{2}+\Theta(t_{m})^{2}]=m^{2}f^{2}\quantity(\Gamma\quantity(\frac{5}{4})\Theta_{i})^{2}\quantity[\frac{J_{1/4}^{2}(t_{m})+J_{5/4}^{2}(t_{m})}{\sqrt{2t_{m}}}].$
(3.34)

In the late time limit $t_{m}\gg 1$, in other words much after the onset of
oscillations, one can use asymptotic forms of the Bessel functions to show
that

$\rho_{\Theta}(t_{m}\gg 1)\approx
m^{2}f^{2}\frac{\sqrt{2}}{\pi}\quantity(\Gamma\quantity(\frac{5}{4})\Theta_{i})^{2}t_{m}^{-3/2}\propto
a^{-3},$ (3.35)

so the ALP field scales as cold dark matter at late times. Redshifting this
energy density until today we obtain

$\Omega_{\Theta,0}\approx\quantity(5.2\times
10^{-4})\quantity(\frac{g_{\rho}^{3/4}(T_{m})}{g_{s}(T_{m})})\sqrt{\frac{m}{H_{0}}}\quantity(\frac{f}{M_{\rm
pl}})^{2}\Theta_{i}^{2},$ (3.36)

where $T_{m}$ is the temperature at which $t_{m}=1$, i.e. $2H(T_{m})=m$, and
$g_{\rho}(T_{m}),g_{s}(T_{m})$ are the number of effective degrees of freedom
in the energy density and entropy respectively.

If the initial angle $\Theta_{i}$ is not close to the minimum, the expression
(3.36) receives corrections from the anharmonacity of the potential. For a
general initial angle, the relic density today is given by

$\boxed{\Omega_{\Theta,0}=\quantity(5.2\times
10^{-4})\quantity(\frac{g_{\rho}^{3/4}(T_{m})}{g_{s}(T_{m})})\sqrt{\frac{m}{H_{0}}}\quantity(\frac{f}{M_{\rm
pl}})^{2}\Theta_{i}^{2}\mathcal{F}_{\rm anh}(\Theta_{i}),}$ (3.37)

where $\mathcal{F}_{\rm anh}$ is the anharmonicity correction which should
behave as $\mathcal{F}_{\rm anh}(\Theta_{i})\rightarrow 1$ as
$\Theta_{i}\rightarrow 0$. By solving the equation of motion (3.5) numerically
with different initial angles, we found that the fit function

$\mathcal{F}_{\rm anh}^{\rm
fit}(\Theta_{i})=\quantity[1-\ln\quantity(1-\quantity(\frac{\absolutevalue{\Theta_{i}}}{\pi})^{3.829})]^{1.395},$
(3.38)

approximates the numerical solution to an accuracy better than $10\%$ for most
of the initial values. This function together with the numerically calculated
values and the relative errors are shown on the right plot of Figure 5.

Figure 5: Left figure: The evolution of the homogeneous mode in SMM/LMM with
different initial conditions obtained by numerically solving (3.10). When the
initial angle is very close to the top of the cosine potential, the ALP field
starts oscillating much later due to to the tiny potential gradient at the
top. The onset of oscillations are shown by vertical dashed lines, and
calculated via (3.39). Right figure: Upper plot shows the numerically obtained
anharmonicity factor values, and the plot of the fit function (3.38), while
the bottom plot show the relative errors of the fit function. In most of the
parameter space, the relative error is less than $10\%$.

As the initial angle becomes closer to the top of the potential, the
anharmonicity factor increases, thereby the relic density also increases. This
behavior is due to the fact that the tiny potential gradient at the top delays
the onset of oscillations as can be seen from the left plot of Figure 5. For
small displacements from the top $\absolutevalue{\pi-\Theta_{i}}\ll 1$, the
dimensionless time at which the oscillations start can be approximated by
Arvanitaki:2019rax

$t_{m}^{\rm
osc}=\ln\quantity[\frac{1}{\absolutevalue{\pi-\Theta_{i}}}\frac{2^{1/4}\pi^{1/2}}{\Gamma(5/4)}],$
(3.39)

which we show via dashed vertical lines in the left plot of Figure 5. As a
result, the duration of the matter-like scaling decreases which enhances the
relic density. Soon we will show that this delay of oscillations also
enchances the ALP fluctuations just like in the Kinetic Misalignment case.

#### Initial conditions for the mode functions

In the case of Standard/Large Misalignment Mechanism, the background ALP field
is frozen at early times, and thus behaves as dark energy. At the zeroth
order, the adiabatic initial conditions for the ALP perturbations do vanish
Marsh:2015xka . Higher order corrections can be derived by solving the
equations of motion for the homogeneous mode (3.10), and the mode functions
(3.11) at early times $t_{m}\ll 1$. The calculation is somewhat technical and
is presented in Appendix A.2, leading to equation (A.17), the early-time
super-horizon behavior of the mode functions:
$\text{SMM/LMM:}\quad\theta_{k}(t_{m}\ll
1)=\frac{2}{5}\Phi_{k}(0)\sin\Theta_{i}\quantity[t_{m}^{2}-\frac{22}{105}\widetilde{k}^{2}t_{m}^{3}+\mathcal{O}(t_{m}^{4})].$
(3.40)

#### The results and comparison with KMM

Figure 6: Power spectra of ALP fluctuations in the Standard/Large Misalignment
with different initial conditions. These results are independent of the ALP
mass and the decay constant, provided that the oscillations start deep in the
radiation era. We observe that if the initial angle is significantly close to
the top of the potential $\absolutevalue{\pi-\Theta_{i}}\lesssim 10^{-5}$,
there is an enhancement in the power spectrum around the modes
$\widetilde{k}\sim\mathcal{O}(1)$.

We first show the dimensionless power spectra at $t_{m}=m/2H=2\times 10^{3}$
with different initial conditions in Figure 6. These results depend only on
the initial conditions, and they are independent of the ALP mass and the ALP
decay constant provided that the oscillations start deep in the radiation era.
From the plot, we can observe that for a large range of initial values, the
power spectra have similar features. The low momentum modes $\widetilde{k}\ll
1$ behave as CDM, while larger modes $\widetilde{k}\gg 1$ are suppressed with
respect to CDM. Only when the initial angle is significantly close to the top
of the potential $\absolutevalue{\pi-\Theta_{i}}\lesssim 10^{-5}$, there is an
enhancement in the modes $\widetilde{k}\sim\mathcal{O}(1)$. The power spectrum
reaches to $\mathcal{O}(1)$ values when
$\absolutevalue{\pi-\Theta_{i}}\lesssim 10^{-10}$, after which the linear
perturbation theory becomes unreliable.

Figure 7: Comparison of the dimensionless power spectra of ALP fluctuations at
$t_{m}=m/2H=2\times 10^{3}$ in the Kinetic Misalignment mechanism (solid
lines), with Large Misalignment Mechanism (dashed lines). For all curves we
assume a constant ALP mass $m=10^{-15}\,$\mathrm{eV}$$. For LMM benchmarks,
the initial angles are determined via (3.37) such that the relic density is
the same as the dark matter density. For $f=6.0\times 10^{14}\,$\mathrm{GeV}$$
(blue lines) both power spectra are similar except there is a suppression of
the high-momentum modes in LMM which is absent in KMM. For $f=1.3\times
10^{14}\,$\mathrm{GeV}$$ (green lines), both models predict an enhancement for
the modes around $\widetilde{k}\sim\mathcal{O}(1)$, however the growth is much
larger in KMM. The reasons behind both of these features are explained in the
main text.

We now compare the power spectra in the Large and Kinetic misalignment
mechanisms. The results are shown in Figure 7 which assumes a constant ALP
mass equal to $m=10^{-15}\,$\mathrm{eV}$$. We have chosen two benchmark values
for the ALP decay constant; $f=6.0\times 10^{14}\,$\mathrm{GeV}$$, and
$f=1.3\times 10^{14}\,$\mathrm{GeV}$$. The former is close to the boundary of
the SMM-KMM transition in the language of Section 2888In Section 2, we defined
the SMM region as the region where $m_{\ast}/H_{\ast}<3$. This is slightly
different what we mean as SMM in this section. In the former, the ALP field
has an initial kinetic energy, but it is insufficient to delay the onset of
oscillations, while in the latter there is no initial kinetic energy. Despite
of this difference, there is no difference from a phenomenological
standpoint., while the latter is close to the boundary of incomplete-complete
fragmentation. To compare the power spectra in these two benchmark points with
their LMM counterparts, we kept the ALP decay constant the same, and
determined the initial angle $\Theta_{i}$ via (3.37) and (3.38) such that all
of dark matter is made of ALPs. For $f=6.0\times 10^{14}\,$\mathrm{GeV}$$
(blue lines), both models predict a similar power spectra except higher
momentum modes are suppressed in the case of LMM, while they are not in KMM.
On the other hand, for $f=1.3\times 10^{14}\,$\mathrm{GeV}$$ (green lines)
both models predict an enhancement of the modes with
$\widetilde{k}\sim\mathcal{O}(1)$, however the enhancement in the case of KMM
is much larger compared to LMM.

Figure 8: Comparison of the evolution of the gauge-invariant density contrast
of the mode $\widetilde{k}=1.7$ in the KMM (green) and LMM (light green)
models. The thin lines denote the numerical solution while the thick lines
show the analytically obtained super-horizon behavior given in (3.45) for KMM
and (3.46) for LMM. The brown dotted line in the middle shows the time at
which the mode enters the horizon. The green line on the right shows the time
of trapping in KMM, while the nearby light green dashed line is the onset of
oscillations in LMM. The dotted light green line shows the point when
$\widetilde{k}^{2}/t_{m}=1$ indicating the start of the tachyonic instability.
We see that even though the enhancement is much larger in LMM when the initial
angle is very close to the top, the final value is smaller than its KMM
counterpart due to the heavy suppression of the initial fluctuations from
small $\absolutevalue{\pi-\Theta_{i}}$, see (3.48).

The larger enhancement in KMM compared to LMM deserves an explanation. At
first, one might think that this is due to the fact that the onset of
oscillations is not delayed by the same amount. However this is not true.
Recall that in KMM the oscillations start at $\tau=2H_{\ast}t$ by definition.
In terms of $t_{m}$ this is

$t_{m}^{\rm osc}=\frac{m}{2H_{\ast}},\quad\text{KMM}.$ (3.41)

On the other hand, the onset of oscillations in LMM is given by (3.39). Then
for $f=1.3\times 10^{14}\,$\mathrm{GeV}$$ one gets $t_{m}^{\rm osc}\approx
13.4(12.7)$ for KMM(LMM), so in both cases the onset of oscillations are
delayed roughly by the same amount. The real reason behind the difference is
imprinted on the evolution of the density contrast before the onset of
oscillations, and can be studied analytically. We will make the comparison
with the gauge-invariant density contrast

$\delta^{\rm GI}=\delta-3\mathcal{H}(1+w)v,$ (3.42)

where $v$ is the _velocity potential_ which for the ALP field takes the form
Eroncel:2022vjg

$v=-\frac{f^{2}\dot{\Theta}\delta\theta}{a\rho_{\Theta}(1+w)}.$ (3.43)

Thus, in Fourier space the gauge-invariant density contrast is

$\delta_{k}^{\rm
GI}=\delta_{k}-3H\frac{f^{2}\dot{\Theta}\theta_{k}}{\rho_{\Theta}}.$ (3.44)

By using the early-time behavior of the KMM mode functions (3.17) and the
expression for the curvature perturbations (3.8), we can show that at leading
order

$\delta^{\rm GI}_{\rm
KMM}\approx\Phi_{k}(0)\times\begin{cases}12\,\widetilde{k}^{2}t_{m},&t_{k}\ll
1,\\\ 3t_{k}\sin(\sqrt{3}t_{k})&t_{k}\gg 1\end{cases}.$ (3.45)

In the SMM/LMM case, we can use (A.17) and (A.25) to obtain at leading order

$\delta^{\rm GI}_{\rm
SMM/LMM}\approx\Phi_{k}(0)\times\begin{cases}\dfrac{16}{175}\cos^{2}\quantity(\dfrac{\Theta_{i}}{2})\widetilde{k}^{2}t_{m}^{3},&t_{k}\ll
1\\\\[10.00002pt]
\dfrac{81}{100}\quantity(\dfrac{\sin^{2}\Theta_{i}}{\cos\Theta_{i}-1})\widetilde{k}^{-4}t_{k}^{2}\cos(t_{k}),&t_{k}\gg
1\end{cases}.$ (3.46)

By roughly comparing the density contrast for a mode $\widetilde{k}\sim 1$ at
$t_{m}\sim 1$ we find  $\frac{\delta_{\rm SMM/LMM}^{\rm GI}}{\delta_{\rm
KMM}^{\rm GI}}\sim 8\times 10^{-3}\cos^{2}\quantity(\frac{\Theta_{i}}{2}).$
(3.47) The enhancement in LMM happens when $\absolutevalue{\pi-\Theta_{i}}\ll
1$. In this case, on top of the two orders of magnitude numerical suppression,
there is also a drastic suppression from the cosine term:  $\frac{\delta_{\rm
LMM}^{\rm GI}}{\delta_{\rm KMM}^{\rm GI}}\sim 8\times
10^{-3}\times\frac{\absolutevalue{\pi-\Theta_{i}}^{2}}{4},\quad\absolutevalue{\pi-\Theta_{i}}\ll
1.$ (3.48) This behavior can clearly be seen in Figure 8. In this plot we
compare the evolution of the gauge-invariant density contrast of the mode
$\widetilde{k}=1.7$ in KMM and LMM models. The thin lines denote the numerical
solution, while the thick lines use super-horizon estimates (3.45) and (3.46).
We see that these estimates predict the evolution quite accurately until the
mode enters the horizon which is shown by the brown dotted line in the middle.
Around $t_{m}=1$, the LMM density constrast is around $12$-orders of magnitude
smaller compared to KMM, consistent with the estimate (3.48). We also see that
despite the fact that the final value of the density contrast is smaller in
the LMM case, the enhancement is significantly larger. The explanation for
this is the tachyonic instability. Recall that the mode function equation of
motion has an effective frequency term given by
$\omega_{k}^{2}=\widetilde{k}^{2}/t_{m}-\cos\Theta$. During the period when
the homogeneous mode is stuck at the top we have
$\cos\Theta\approx\cos\Theta_{i}\approx-1$. So the effective frequency of the
mode becomes tachyonic when $\widetilde{k}^{2}/t_{m}>1$. This moment is shown
by the vertical dotted light green line in Figure 8 labeled as “tachyonic”. We
see that the exponential growth starts shortly after this point.

Before closing this sub-section we also explain briefly why there is a
suppression in SMM/LMM power spectra at higher momentum values. These modes
are well sub-horizon when the oscillations start so their values slightly
before the oscillations can be read from (3.46) as $\delta^{\rm
GI}\propto\widetilde{k}^{-4}t_{k}^{2}\propto\widetilde{k}^{-2}t_{m}$.
Therefore, higher momentum modes are suppressed with $\widetilde{k}^{-2}$
before the onset of oscillations. After the oscillations start their
amplitudes are approximately constant which explains the suppression at high
momentum modes in SMM/LMM.

##  4 ALP density perturbations in the complete fragmentation regime

If the fragmentation is complete, the system dynamics is non-linear, and the
equation of motion for the mode functions become coupled to each other. In
this case, one cannot use the cosmological perturbation theory as we did in
the previous section. To study this regime precisely, one needs to use non-
perturbative methods such as lattice simulations which we leave for future
work. However, it is still possible to get some semi-analytical estimates by
approximating the axion potential at late times by

$V(\theta)\approx\frac{1}{2}f^{2}m^{2}\theta^{2};$ (4.1)

Even though this approximation clearly breaks down during the fragmentation,
it will eventually be a good one once the self-interactions become
negligible999One feature which cannot be studied with this approximation is
the formation of _oscillons_ which are localized field configurations of a
real scalar field sustained only by self-interactions Bogolyubsky:1976yu ;
Gleiser:1993pt ; Copeland:1995fq . For more information about the oscillons
and their consequences for dark matter see Olle:2019kbo and the references
therein.. Since what we need is the power spectrum at late times, this
approximation is adequate for our estimates.

In the case when the fluctuations are not necessarily small, the density
contrast is defined by

$\delta(t,\mathbf{x})\equiv\frac{\rho(t,\mathbf{x})-\expectationvalue{\rho}}{\expectationvalue{\rho}},$
(4.2)

where $\expectationvalue{\rho}$ is the average energy density. For a quadratic
potential (4.1) it takes the form

$\begin{split}\delta=\frac{f^{2}}{\expectationvalue{\rho}}\Bigg{\\{}&\frac{1}{2}\quantity[\dot{\theta}^{2}\quantity(1+2\Phi)-\expectationvalue{\dot{\theta}^{2}}]+\frac{1}{2a^{2}}\quantity[\gradient{\theta}\cdot\gradient{\theta}\quantity(1-2\Phi)-\expectationvalue{\gradient{\theta}\cdot\gradient{\theta}}]\\\
&+m^{2}\quantity(\theta^{2}-\expectationvalue{\theta^{2}})\Bigg{\\}}.\end{split}$
(4.3)

In deriving this expression we have neglected the $\mathcal{O}(\Phi^{2})$
terms, but the ALP field $\theta(t,\mathbf{x})$ is kept non-perturbative. We
now expand the ALP field as

$\theta(t,\mathbf{x})=\theta_{B}(t)+\theta_{F}(t,\mathbf{x})=\theta_{B}(t)+\int_{k>k_{\ast}}\frac{\differential[3]{k}}{(2\pi)^{3}}\hat{\theta}_{\mathbf{k}}\theta_{k}(t)e^{-i\mathbf{k}\cdot\mathbf{x}}.$
(4.4)

Here $\hat{\theta}_{\mathbf{k}}$’s are stochastic variables that carry the
statistical properties, while $\theta_{k}$’s are $c$-number functions which
depend only on the magnitude of the momentum $k\equiv\mathbf{k}$. We assume
that $\hat{\theta}_{\mathbf{k}}$’s satisfy

$\expectationvalue{\hat{\theta}_{\mathbf{k}}\hat{\theta}_{\mathbf{k^{\prime}}}^{\ast}}=(2\pi)^{3}\delta^{(3)}(\mathbf{k}-\mathbf{k}^{\prime}).$
(4.5)

The function $\theta_{B}$ describes the evolution of the long wavelength modes
$k<k_{\ast}$ for which $\gradient{\theta}\approx 0$. We took
$k_{\ast}=a_{\ast}H_{\ast}$, but this choice does not affect our results as
long as $k_{\ast}$ is much smaller than the relevant modes for the
fragmentation for which $k\sim ma_{\ast}$ Eroncel:2022vjg . The reason why we
have made such a split even in the case of complete fragmentation is the
following: When calculating the power spectrum in the case of complete
fragmentation, we will use the semi-analytical results for the mode functions
that we have derived in Eroncel:2022vjg . This method does not work for the
modes that are super-horizon during fragmentation, i.e.
$k_{\ast}<a_{\ast}H_{\ast}$.

By plugging the expansion (4.4) into (4.3) we find

$\begin{split}\delta=&\frac{f^{2}}{\expectationvalue{\rho}}\quantity[\frac{\theta_{B}^{\prime}\theta_{F}^{\prime}}{a^{2}}\quantity(1+2\Phi)+m^{2}\theta_{B}\theta_{F}+\Phi\frac{\theta_{B}^{2}}{a^{2}}]\\\
+&\frac{f^{2}}{\expectationvalue{\rho}}\Bigg{\\{}\frac{1}{2a^{2}}\quantity[\theta_{F}^{\prime
2}\quantity(1+2\Phi)-\expectationvalue{\theta_{F}^{\prime
2}}]+\frac{1}{2a^{2}}\quantity[\gradient{\theta_{F}}\cdot\gradient{\theta_{F}}\quantity(1-2\Phi)-\expectationvalue{\gradient{\theta_{F}}\cdot\gradient{\theta_{F}}}]\\\
&\qquad+\frac{m^{2}}{2}\quantity(\theta_{F}^{2}-\expectationvalue{\theta_{F}^{2}})\Bigg{\\}}\end{split}$
(4.6)

By approximating $1\pm 2\Phi\approx 1$, the density contrast reduces to

$\begin{split}\delta=&\frac{f^{2}}{\expectationvalue{\rho}}\quantity[\frac{\theta_{B}^{\prime}\theta_{F}^{\prime}}{a^{2}}+m^{2}\theta_{B}\theta_{F}+\Phi\frac{\theta_{B}^{2}}{a^{2}}]\\\
+&\frac{f^{2}}{\expectationvalue{\rho}}\Bigg{\\{}\frac{1}{2a^{2}}\quantity[\theta_{F}^{\prime
2}-\expectationvalue{\theta_{F}^{\prime
2}}]+\frac{1}{2a^{2}}\quantity[\gradient{\theta_{F}}\cdot\gradient{\theta_{F}}-\expectationvalue{\gradient{\theta_{F}}\cdot\gradient{\theta_{F}}}]+\frac{m^{2}}{2}\quantity(\theta_{F}^{2}-\expectationvalue{\theta_{F}^{2}})\Bigg{\\}}\end{split}$
(4.7)

We see that the density contrast has two contributions. The first term
contains terms which are linear in fluctuations, and it reduces to the linear
density constrast (3.4) for a quadratic potential if the fragmentation is
incomplete since in this case $\theta_{B}\rightarrow\Theta$, and
$\expectationvalue{\rho}\approx\rho_{\Theta}$. Therefore, we refer this
contribution to the density contrast as $\delta^{\rm lin}$. Since the second
term is quadratic in fluctuations, we refer it as $\delta^{\rm quad}$, and
write $\delta=\delta^{\rm lin}+\delta^{\rm quad}$. On the scales
$k\sim\mathcal{O}(1)\times ma_{\ast}$ where the fragmentation is efficient,
the quadratic contribution $\delta^{\rm quad}$ dominates over the linear
contribution $\delta^{\rm lin}$. However, at scales $k\ll ma_{\ast}$ the
quadratic contribution is small as we shall show soon, so the linear
contribution dominates.

After this splitting, the power spectrum takes the form

$P_{\delta}(k)\equiv\frac{1}{\mathcal{V}}\expectationvalue{\absolutevalue{\delta_{k}}^{2}}=\frac{1}{\mathcal{V}}\quantity[\expectationvalue{\absolutevalue{\delta_{k}^{\rm
lin}}^{2}}+\expectationvalue{\absolutevalue{\delta_{k}^{\rm
quad}}^{2}}+2\expectationvalue{\Re{\quantity(\delta_{k}^{\rm
lin})^{\ast}\delta_{k}^{\rm quad}}}],$ (4.8)

where $\mathcal{V}$ is a volume factor which is much larger than all scales in
the problem. The last term in this expression does contain terms of the form
$\expectationvalue{\hat{\Phi}\hat{\theta}\hat{\theta}}$ where $\hat{\Phi}$ and
$\hat{\theta}$ are the operators which encode the statistical properties of
the curvature perturbations and ALP fluctuations respectively. If the latter
are sourced by adiabatic perturbations, then $\hat{\theta}\sim\hat{\Phi}$.
Therefore, by assuming that the three-point function for the curvature
perturbations vanishes $\expectationvalue{\hat{\Phi}\hat{\Phi}\hat{\Phi}}=0$,
the last term in (4.8) becomes zero. Then the power spectrum becomes

$P_{\delta}(k)=P_{\delta}^{\rm lin}(k)+P_{\delta}^{\rm quad}(k),$ (4.9)

where $P_{\delta}^{\rm
lin}(k)=\mathcal{V}^{-1}\expectationvalue{\absolutevalue{\delta_{k}^{\rm
lin}}^{2}}$ and $P_{\delta}^{\rm
quad}(k)=\mathcal{V}^{-1}\expectationvalue{\absolutevalue{\delta_{k}^{\rm
quad}}^{2}}$. The linear contribution is given by

$P_{\delta}^{\rm
lin}(k)=\quantity(\frac{f^{2}}{\expectationvalue{\rho}})^{2}\quantity(\dot{\theta}_{B}\dot{\theta}_{k}+m^{2}\Theta_{B}\theta_{k}+\Phi_{k}\dot{\theta}_{B}^{2})^{2},$
(4.10)

and it reduces to the power spectrum in the linear case when
$\theta_{b}\rightarrow\Theta$, and
$\expectationvalue{\rho}\rightarrow\rho_{\Theta}$. The quadratic contribution
is derived in Enander:2017ogx , and it reads

$\boxed{P_{\delta}^{\rm
quad}(k)=2(2\pi)^{3}\dfrac{\int\differential[3]{q}\absolutevalue{F(q,q-k)}^{2}}{\quantity(\int\differential[3]{q}F(q,q))^{2}},}$
(4.11)

where

$\boxed{F(k,k^{\prime})=\dot{\theta}_{k}\dot{\theta}_{k^{\prime}}+\quantity(\frac{\mathbf{k}\cdot\mathbf{k}^{\prime}}{a^{2}}+m^{2})\theta_{k}\theta_{k^{\prime}},}$
(4.12)

The relevant contribution for the late time observables is the quadratic one.
So we discuss it first, and then comment on the linear contribution.

###  4.1 Quadratic contribution

We are interested in the late time behavior of the power spectrum, when all
modes of interest are non-relativistic. In the non-linear regime (complete
fragmentation) we will use the following ansatz to express the mode functions
Eroncel:2022vjg :

$\theta_{k}(t)=\sqrt{P_{\theta}(k)}A_{k}(t)N_{k}(t)\cos(mt+\varphi_{k}),$
(4.13)

where $\varphi_{k}$ is a momentum-dependent phase factor, $P_{\theta}(k)$ is
the initial power spectrum before the fragmentation

$P_{\theta}(k)=\frac{2\pi^{2}}{k^{3}}\quantity(\frac{1}{3})^{2}A_{s}\quantity(\frac{\dot{\Theta}_{i}^{2}}{H_{i}^{2}}),$
(4.14)

and the redshift factor $A_{k}(t)$ is given by

$A_{k}(t)=\frac{\omega_{k}^{1/2}(t_{i})a_{i}^{3/2}}{\omega_{k}^{1/2}(t)a^{3/2}}\approx\frac{(k/a_{i})^{1/2}a_{i}^{3/2}}{m^{1/2}a^{3/2}}.$
(4.15)

By definition, the backreaction will be important in the non-linear regime so
the amplification factors $N_{k}$ are calculated via the workflow outlined in
Section 3.2 of Eroncel:2022vjg . By plugging these terms into the expression
for the power spectrum (4.11) we obtain

$P_{\delta}^{\rm
quad}(k)=\pi\quantity(\int\differential{q}N_{q}^{2})^{-2}\int\differential[3]{q}\frac{N_{q}^{2}N_{\absolutevalue{\mathbf{q}-\mathbf{k}}}^{2}\cos^{2}\quantity(\varphi_{q}-\varphi_{\absolutevalue{\mathbf{q}-\mathbf{k}}})}{q^{2}\absolutevalue{\mathbf{q}-\mathbf{k}}^{2}},$
(4.16)

where momentum integrals are over the modes which are amplified during the
fragmentation, hence they are finite. The amplification factors $N_{k}$ are
interpreted as their asymptotic limits, so they and the whole power spectrum
is independent of time101010This is true only when the potential is
approximated by a quadratic, in other words when the self-interactions are
neglected.. By averaging over the phase factors we can replace the cosine
square term by $1/2.$ Also, without loss of generality, we can choose
$\mathbf{q}$ to be aligned with the $\hat{\mathbf{z}}$-direction. Then
$\absolutevalue{\mathbf{k}-\mathbf{q}}=\sqrt{k^{2}+q^{2}-2kqu}$ where $u$ is
the cosine of the angle between $\mathbf{k}$ and $\mathbf{q}$. After
integrating over the azimuthal angle we find

$P_{\delta}^{\rm
quad}(k)=\pi^{2}\quantity(\int\differential{q}N_{q}^{2})^{-2}\int\differential{q}N_{q}^{2}\int_{-1}^{1}\differential{u}\frac{N_{\absolutevalue{\mathbf{q}-\mathbf{k}}}^{2}}{\absolutevalue{\mathbf{q}-\mathbf{k}}^{2}}.$
(4.17)

It is more convenient to make a change of parameter from $u$ to
$p\equiv\absolutevalue{\mathbf{q}-\mathbf{k}}$. This way we obtain

$P_{\delta}^{\rm
quad}(k)=\frac{1}{2}\frac{2\pi^{2}}{k}\quantity(\int\differential{q}N_{q}^{2})^{-2}\int\frac{\differential{q}}{q}N_{q}^{2}\int_{\absolutevalue{q-k}}^{q+k}\frac{\differential{p}}{p}N_{p}^{2}.$
(4.18)

The dimensionless power spectrum is  $\mathcal{P}_{\delta}^{\rm
quad}(k)=\frac{k^{3}}{2\pi^{2}}P_{\delta}^{\rm
quad}(k)=\frac{1}{2}k^{2}\quantity(\int\differential{q}N_{q}^{2})^{-2}\int\frac{\differential{q}}{q}N_{q}^{2}\int_{\absolutevalue{q-k}}^{q+k}\frac{\differential{p}}{p}N_{p}^{2}.$
(4.19) Let $\mathfrak{q}\equiv q/(a_{\ast}m)$. Then, one can show that at
small momenta  $\mathcal{P}_{\delta}^{\rm quad}(\kappa\ll
1)\approx\kappa^{3}\quantity(\int\differential{\mathfrak{q}}N_{\mathfrak{q}}^{2})^{-2}\int\frac{\differential{\mathfrak{q}}}{\mathfrak{q}^{2}}N_{\mathfrak{q}}^{4}\propto\kappa^{3},$
(4.20) where $\kappa\equiv k/(a_{\ast}m)$. Therefore we found that the power
spectrum approaches to that of a white noise at low momenta. On the other
hand, the large momentum behavior is  $\mathcal{P}_{\delta}^{\rm
quad}(\kappa\gg
1)\approx\kappa\quantity(\int\differential{\mathfrak{q}}N_{\mathfrak{q}})^{-1}$
(4.21) These modes collapse too late to contribute to structure formation and
will thus be not too relevant. We note that this result assumes that for all
the modes the initial power spectrum is given by (4.14) which is true as long
as the mode is super-horizon at the onset of the kination-like scaling of the
homogeneous mode. Therefore, the initial conditions are not valid for the
modes $k>a_{\rm kin}H_{\rm kin}$. As a result, the high-momenta behavior will
also be modified. When and how this modification happens depends on the UV
completion, and is irrelevant for the gravitational signatures that we will
study in the next section. Thus, we will not comment on the high-momenta
behavior any further.

Figure 9: The dimensionless power spectrum $\mathcal{P}_{\delta}^{\rm
quad}(k)$ for an ALP model with constant mass $m=10^{-15}\;$\mathrm{eV}$$ as a
function of the dimensionless momentum $\kappa=k/(ma_{\ast})$ for different
values of the decay constant. The power spectra is plotted on the linear scale
in the left figure, and on the log scale in the right figure. The dashed lines
show the low-momenta limit of the power spectrum
$\mathcal{P}\propto\kappa^{3}$ calculated via (4.20), while the dot-dashed
lines show the high-momenta limit $\mathcal{P}\propto\kappa$ calculated via
(4.21). In all cases, ALP make all of dark matter.

We show the plots of the quadratic contribution to the power spectrum in
Figure 9 for a couple of benchmark points. We assumed a constant ALP mass
$m=10^{-15}\,$\mathrm{eV}$$, and varied the ALP decay constant $f$. For all
the curves, the ALP makes all of dark matter. We can see that all power
spectra are peaked at $\kappa\sim\mathcal{O}(1)$ while the peak position moves
to slightly higher momenta for smaller decay constants. The reason behind this
is that the fragmentation starts earlier for smaller decay constants which
permits exponential growth of higher modes. However, this also causes
amplification of a larger range of modes which results in a more broad
spectrum that has a smaller peak value. We can also see that the low/high
momentum behavior of all spectra is consistent with the expressions (4.20) and
(4.21) respectively.

###  4.2 Linear contribution

For fragmented axions to explain all of dark matter, the resulting matter
power spectrum should be consistent with the observations around the CMB
scales $k\sim 0.1\,\rm{Mpc}^{-1}$. The relevant modes for fragmentation are
those with $k\sim k_{\rm peak}\sim a_{\ast}m$. To express this in
dimensionless units, we use the following relation derived in Eroncel:2022vjg
:

$\frac{m_{\ast}}{m_{0}}\quantity(\frac{a_{\ast}}{a_{0}})^{3}=\frac{\pi}{8}\frac{\rho_{\Theta,0}}{\Lambda_{b,0}^{4}},$
(4.22)

where $m_{\ast}$ is the ALP mass at trapping, $\Lambda_{b,0}=\sqrt{m_{0}f}$ is
the today’s value of the barrier height, and $\rho_{\Theta,0}$ is the ALP
energy density today without taking fragmentation into account111111The
fragmentation might slightly modify the relic density today. We made a semi-
analytical estimate in Eroncel:2022vjg , and found that this modification is
not significant.. From (4.22) and assuming a constant mass we get

$k_{\rm peak}\sim
10^{5}\,\mathrm{Mpc}^{-1}\times\quantity(\frac{m}{10^{-15}\,$\mathrm{eV}$})^{1/3}\quantity(\frac{f}{10^{14}\,$\mathrm{GeV}$})^{-2/3}.$
(4.23)

Therefore, the fragmentation scales are much shorter compared to the CMB
ones121212The fragmentation scales can be close to the CMB scales at the very
low ALP mass and very high ALP decay constant region of the parameter space.
However, this region is excluded by the BBN constraints Eroncel:2022vjg .. As
a result, the quadratic contribution will be very suppressed due to the
$k^{3}$ decay at low momenta. Thus, the matter power spectrum at CMB scales
will be determined mainly by the linear contribution.

We cannot evaluate the linear contribution since $\theta_{B}$ is not known.
However, from the discussion of the linear case in Section 3, we expect that
the modes which are super-horizon during fragmentation will lock on the CDM
prediction once they enter the horizon. Therefore, we approximate the linear
contribution of the low momentum modes by the power spectrum of CDM:

$P_{\delta}^{\rm lin}(k;z)\approx P_{\delta,\rm CDM}(k;z).$ (4.24)

This approximation has been used in Lee:2020wfn to study the gravitational
signatures of axion miniclusters when Peccei-Quinn breaking happens after
inflation, and also in Graham:2015rva ; Alonso-Alvarez:2018tus where dark
matter is generated from inflationary perturbations.

Figure 10: The combined dimensionless power spectrum
$\mathcal{P}_{\delta}=\mathcal{P}_{\delta,\rm CDM}+\mathcal{P}_{\delta}^{\rm
quad}$ for an ALP model with constant mass $m=10^{-15}\,$\mathrm{eV}$$ as a
function of the comoving momentum $k$ for different values of the decay
constant. The dashed lines show the low-momenta behavior without the CDM
contribution, while the dotted lines show the extrapolation to high momenta
using (4.21). Note that the CDM contribution is added only to the low-momenta
part of the power spectrum, and we assumed $a/a_{0}=10^{-6}$ when calculating
the CDM spectrum which is much later than the trapping scale factor $a_{\ast}$
for all benchmark points.

We show the full power spectrum including the CDM contribution in Figure 10.
The benchmark points are the same in Figure 9. Even though the quadratic
contribution $\mathcal{P}_{\delta}^{\rm quad}$ is time-independent, the CDM
contribution has a logarithmic time-dependence during the radiation era. We
have evaluated the CDM contribution at $a/a_{0}=10^{-6}$ which is much later
than the trapping scale factor $a_{\ast}$ for all benchmark points. The high
momenta behaviors are shown by the dotted lines since the position of their
cutoff is model-dependent. We can see that the spectra is peaked at higher
comoving momenta for smaller decay constants. This happens because of two
reasons: First, as we have already seen in Figure 9, the power spectrum is
peaked at higher dimensionless momenta for smaller decay constants. On top of
this effect, the trapping happens later for smaller decay constants which
implies that the amplified modes have larger comoving momenta. This shift to
larger momenta will have a significant effect on the gravitational signatures
as we shall see in the next section.

In order to study the gravitational signatures we need to evolve the power
spectrum $\mathcal{P}^{\rm quad}$ until today. For this we will use the
linearized evolution equation (3.27) even though the fluctuations are clearly
non-linear. The reasoning behind this approximation is that the linearized
analysis is all we need to have an estimate of the gravitational signatures
which becomes more clear when we discuss the Press-Schechter formalism in
Section 5.2.

###  4.3 Comparison with the post-inflationary scenario

In our discussion so far, we have always considered the _pre-inflationary_
scenario where the spontaneous symmetry breaking that creates ALPs as Nambu-
Goldstone bosons happens before inflation. In the absence of an initial
kinetic energy, the ALP angle $\theta$ takes random values
$\theta\in[-\pi,\pi)$ in different patches of the universe after symmetry
breaking. If this process happens before inflation, each patch is inflated to
huge sizes so that the initial angle is the same everywhere in the observable
universe which makes the ALP field homogeneous at early times.

However, the symmetry breaking can also happen after inflation. This is the
so-called _post-inflationary_ scenario. In this case the observable universe
will contain many patches all having a different initial angle. This leads to
the emergence of topological defects like cosmic strings and domain walls,
also known as the string-wall network Kibble_1976 ; PhysRevD.26.435 ;
PhysRevLett.48.1867 . As a result, the ALP field becomes very inhomogeneous in
contrast to the pre-inflationary scenario.

It is interesting to compare our predictions for the power spectrum with the
one predicted in the post-inflationary scenario. A simplified linear analysis
of the post-inflationary scenario for the QCD axion is given in
Enander:2017ogx by assuming a quadratic potential. Full non-linear analysis of
ALP models with different temperature dependencies including the constant mass
scenario has been recently presented in OHare:2021zrq . In this section, first
we calculate the power spectrum of a constant mass ALP model by adapting the
method used in Enander:2017ogx , and then make a comparison with the numerical
results of OHare:2021zrq . Finally we compare the predictions in the post-
inflationary scenario with the predictions of KMM and SMM/LMM.

#### Linearized analysis

Following Enander:2017ogx , we separate the ALP mode functions as

$\theta_{k}(t)=\hat{\theta}_{k}f_{k}(t),$ (4.25)

where the statistical properties and the initial conditions are encoded in
$\hat{\theta}_{k}$, while the functions $f_{k}(t)$ contain the time
evolutions, which are normalized such that $f_{k}(t_{i})=1$. For a quadratic
potential, the equations of motion for $f_{k}$’s are

$\ddot{f}_{k}+3H\dot{f}_{k}+\quantity(\frac{k^{2}}{a^{2}}+m^{2})f_{k}=0,$
(4.26)

where $m$ is the ALP mass which is assumed to be constant. The oscillations
start around $t=t_{1}=(2H_{1})^{-1}$, where in the constant mass case $H_{1}$
is defined by Blinov:2019rhb ; OHare:2021zrq

$m=\frac{8}{5}H_{1}.$ (4.27)

We also define $\mathcal{H}_{1}=a_{1}H_{1}=L_{1}^{-1}$ where $L_{1}$ is the
comoving horizon at $t_{1}$ OHare:2021zrq . Then, in terms of the
dimensionless time $t_{m}=mt$, (4.26) takes the form

$f_{k}^{\prime\prime}(t_{m})+\frac{3}{2t_{m}}f_{k}^{\prime}(t_{m})+\quantity[\frac{5}{16}\frac{K^{2}}{t_{m}}+1]f_{k}(t_{m})=0,$
(4.28)

where we have also defined $K\equiv kL_{1}$.

The evolution of the mode functions after the onset of oscillations can well
be described by the WKB approximation:

$f_{k}(t_{m})=\frac{c_{k}}{t_{m}^{3/4}\sqrt{\widetilde{\omega}_{k}(t_{m})}}\cos\quantity(\int_{t_{m}^{\rm
WKB}}^{t_{m}}\differential{t_{m}^{\prime}\widetilde{\omega}_{k}(t_{m}^{\prime})}+\varphi_{k}),\quad\widetilde{\omega}_{k}(t_{m})\equiv\sqrt{\frac{5}{16}\frac{K^{2}}{t_{m}}+1},$
(4.29)

where $c_{k}$ and $\varphi_{k}$ are numerical factors which we have determined
by numerically solving the equation of motion (4.28) until $t_{m}^{\rm WKB}\gg
1$. We have verified that by choosing $t_{m}^{\rm WKB}=50$, the WKB
approximation reproduces the numerical solution very well.

The power spectrum is given by Enander:2017ogx

$P(q)=2\quantity(2\pi)^{3}\frac{\int\differential[3]{k}P_{\theta}(k)P_{\theta}(\absolutevalue{\mathbf{k}-\mathbf{q}})F(k,k-q)^{2}}{\quantity[\int\differential[3]{k}P_{\theta}(k)F(k,k)]^{2}},$
(4.30)

where

$F(k,k^{\prime})=\dot{f}_{k}(t)\dot{f}_{k^{\prime}}(t)+\quantity(\frac{\mathbf{k}\cdot\mathbf{k}^{\prime}}{a^{2}}+m^{2})f_{k}(t)f_{k^{\prime}}(t),$
(4.31)

and $P_{\theta}(k)$ is the field power spectrum at some initial temperature
$T_{i}$ defined such that

$\expectationvalue{\hat{\theta}_{k}\hat{\theta}_{k^{\prime}}^{\ast}}=\quantity(2\pi)^{3}\delta^{(3)}(\mathbf{k}-\mathbf{k}^{\prime})P_{\theta}(k).$
(4.32)

A reasonable choice for $P_{\theta}(k)$ is Enander:2017ogx ; OHare:2021zrq

$\boxed{P_{\theta}(k)=\frac{8\pi^{4}}{3\sqrt{\pi}k_{\rm
cr}^{3}}\exp(-\frac{k^{2}}{k_{\rm cr}^{2}}),}$ (4.33)

where the cutoff is given by $k_{\rm cr}=a_{i}H_{i}$, and the normalization is
fixed such that $\expectationvalue{\theta(\mathbf{x})^{2}}=\pi^{2}/3$. Note
that at late times when all the relevant modes become non-relativistic we have
$\widetilde{\omega}_{k}\approx 1$, and the power spectrum becomes independent
of time131313As we have stated in Section 4.1, this is true only for the
linearized analysis where the potential is approximated by a quadratic, and
therefore the self-interactions are neglected. .

Last thing we need to specify is the initial time. This needs to be before the
onset of oscillations but not too early since at those times the string-wall
network plays an important role Vaquero:2018tib ; OHare:2021zrq . Again
following Enander:2017ogx we have chosen $T_{i}=3T_{1}$ as the initial
temperature where $T_{1}$ is defined by $H(T_{1})=H_{1}$. This corresponds to
an initial time of $t_{m,i}=4/45$ and a momentum cutoff of $k_{\rm
cr}L_{1}=9$. After this point the calculation is identical to the one we have
presented in Section 4.1 so we do not repeat it here.

#### Results and comparison with lattice simulations and other mechanisms

Figure 11: Late-time dimensionless power spectrum of the ALP field in the
post-inflationary scenario with constant mass obtained via the linearized
analysis of Enander:2017ogx (blue lines). We compare this spectrum with the
lattice result of OHare:2021zrq (orange lines). The left and right plots show
the same spectra except for the scalings of the vertical axis. The orange band
on the right plot show the low-momenta behavior reported in OHare:2021zrq such
that $\mathcal{P}(kL_{1}\ll 1)\approx C\quantity(kL_{1})^{3}$ with $C\simeq
0.03\text{ -- }0.05$. The lattice result is obtained at $t_{m}=2\times 10^{3}$
in our units, however we have checked that the linearized result does not
change significantly after $t_{m}\sim 10^{3}$.

We now present our results and make a comparison with the lattice result
obtained in OHare:2021zrq in Figure 11. The blue lines show the spectra
obtained via the linearized analysis, while the orange lines show the lattice
results. The orange band on the right plot is given by $\mathcal{P}(kL_{1}\ll
1)\approx C\quantity(kL_{1})^{3}$ with $C\simeq 0.03\text{ -- }0.05$, which is
the low-momentum behavior reported in OHare:2021zrq . We can see that the non-
linear effects shifts the spectrum to larger momentum values and make the
spectrum broader. As a result, the peak value gets reduced by a factor of
$\approx 3.5$.

Figure 12: Left plot: The line on the ALP parameter space where the correct
relic density is obtained in the post-inflationary scenario as determined by
(4.34). This line agrees perfectly with the one obtained in the Standard/Large
Misalignment Mechanism by setting the initial angle to $\Theta_{i}\approx
3.06$. Right plot: The comparison of the power spectra in the post-
inflationary scenario obtained via the linearized analysis of Enander:2017ogx
, and the lattice calculation of OHare:2021zrq with the spectra obtained in
the Kinetic and Large Misalignment Mechanism for the same benchmark point.

Before concluding this section, it would be instructive to compare the power
spectra of all the dark matter production mechanisms that we have seen so far
on a single plot. For a given ALP mass, the decay constant is fixed in the
post-inflationary scenario. The relic density for this scenario is given by
OHare:2021zrq

$\boxed{h^{2}\Omega_{\theta}\simeq
0.019\quantity(\frac{g_{\rho}(T_{1})}{70})^{3/4}\quantity(\frac{g_{s}(T_{1})}{70})^{-1}\quantity(\frac{m}{$\mathrm{\SIUnitSymbolMicro
eV}$})^{1/2}\quantity(\frac{f}{10^{12}\,$\mathrm{GeV}$})^{2}.}$ (4.34)

One can then solve for the ALP decay constant $f$ to get the correct relic
density. We show this line on the ALP parameter space in the left plot of
Figure 12 with the label “post-inflation”. This line matches perfectly with
the line one gets in the Standard/Large Misalignment Mechanism by setting the
initial angle to $\Theta_{i}\approx 3.06$. However, the power spectrum of the
fluctuations is completely different due to the highly inhomogeneous initial
conditions. We demonstrate this on the right plot of Figure 12. In this plot,
we set the ALP mass to $m=10^{-15}\,$\mathrm{eV}$$, and solve (4.34) to get
the correct relic density and found $f\approx 3.53\times
10^{14}\,$\mathrm{GeV}$$. We then calculated the power spectra predicted by
the Kinetic and the Standard/Large Misalignment Mechanisms and make the
comparison with the post-inflationary spectra obtained both with the
linearized analysis (solid line), and the lattice results (dashed line). We
see that the spectrum is greatly enhanced and broaden in the post-inflationary
scenario.

##  5 Gravitational collapse of fluctuations and the halo spectrum

Once the matter fluctuations become sufficiently dense, they decouple from the
ambient Hubble flow, and form gravitationally bound structures known as
_halos_. This process is called _gravitational collapse_. Studying this
process precisely is quite difficult, and requires N-body simulations.
However, qualitative results can be derived by exploiting the approximate
spherical symmetry. This will be our goal in this section. We start by
deriving the critical density beyond which a density fluctuation will
experience gravitational collapse in Section 5.1. Then we introduce the Press-
Schechter formalism in Section 5.2 which we shall use in Section 5.3 to
compute the halo mass function, the number density of halos per logarithmic
mass bin. We conclude the section by calculating the halo spectrum in Section
5.4 that contains information about the mass-density relation of the halos at
a given benchmark.

###  5.1 Critical density for collapse

We consider a universe which can either be matter or radiation dominated, and
take a spherical shell of _physical radius_ $r$. The time evolution of this
radius obeys Kolb:1994fi

$\frac{\ddot{r}}{r}=-\frac{4\pi G}{3}(\rho+3p)=-\frac{4\pi
G}{3}\quantity(\rho_{m}+\rho_{r}+3p_{r})=-\frac{4\pi G}{3}\rho_{m}-\frac{8\pi
G}{3}\rho_{r},$ (5.1)

where $\rho_{m}$ is the matter density141414Since we are neglecting the
baryons in this work, the matter density is the same as the ALP density,
$\rho_{m}=\rho_{\theta}$.  including the over-densities, and $\rho_{r}$ is the
radiation density which we assume to be homogeneous everywhere. The radiation
pressure is $p_{r}=\rho_{r}/3$, and we assumed that the matter is pressureless
$p_{m}=0$. The latter assumption is accurate for CDM, however it breaks down
for ALPs at small scales due to the quantum pressure discussed in the previous
section. We will comment later on the consequences. For now, we continue our
discussion by assuming $p_{m}=0$.

Let $M$ denote the mass inside the spherical region, i.e.

$M=\frac{4\pi}{3}\rho_{m}r^{3}=\frac{4\pi}{3}\quantity(\overline{\rho}_{m}+\Delta\rho_{m})r^{3}=\frac{4\pi}{3}\overline{\rho}_{m}\quantity(1+\delta)r^{3},$
(5.2)

where $\overline{\rho}_{m}$ is the background matter density, and
$\delta\equiv\Delta\rho_{m}/\overline{\rho}_{m}$ is the size of the
fluctuation. We will assume that this mass remains constant during the
gravitational collapse. In terms of $M$, the equation of motion for $r$
becomes

$\ddot{r}=-\frac{8\pi G}{3}\rho_{r}r-\frac{MG}{r^{2}}.$ (5.3)

Let us parametrize $r(t)$ by

$r(t)=a(t)b_{R}(t)R,$ (5.4)

where $R$ is the _comoving_ radial coordinate, and $b_{R}$ denotes the
deviation from the ambient Hubble flow. The $R$-subscript indicates that the
evolution is valid for the spherical shell with comoving label $R$. We will
drop the subscript below to simplify the notation.

Next, we will switch the conformal time, and make a change of variables to
$y\equiv a/a_{\rm eq}$. After some algebra, one can obtain a differential
equation for the evolution of $b(y)$: Kolb:1994fi

$y(1+y)\derivative[2]{b}{y}+\quantity(1+\frac{3}{2}y)\derivative{b}{y}+\frac{1}{2}\quantity(\frac{1+\delta_{i}}{b^{2}}-b)=0,$
(5.5)

where $\delta_{i}$ is the _initial_ over-density in the region
$\delta_{i}=\delta\rho_{m}(y_{i})/\overline{\rho}_{m}(y_{i})$ set at some
early time when $b(y_{i})=1$. By deriving this result, we have used

$M=\frac{4\pi}{3}\rho_{m,i}r_{i}^{3}=\frac{4\pi}{3}\overline{\rho}_{m,i}a_{i}^{3}\quantity(1+\delta_{i})R^{3}=\frac{4\pi}{3}\frac{\overline{\rho}_{\rm
eq}}{2}a_{\rm eq}^{3}\quantity(1+\delta_{i})R^{3},$ (5.6)

where $\overline{\rho}_{\rm eq}=2\overline{\rho}_{m,\rm
eq}=2\overline{\rho}_{r,\rm eq}$ is the _total_ average energy density at the
matter-radiation equality.

For a given value of the initial over-density $\delta_{i}$, one can solve
(5.5) numerically with the initial conditions $b(y_{i})=1$, and
$b^{\prime}(y_{i})=0$ where $y_{i}\ll 1$. Initially, the radius will expand
with the Hubble flow. Then it will start to deviate, and at $y_{\rm ta}$ it
_turns around_ , i.e. the expansion will stop, and the radius will start
shrinking. This happens when $\dot{r}=0$, or $b^{\prime}(y_{\rm ta})+b(y_{\rm
ta})/y_{\rm ta}=0$. Numerical solution shows that the turnaround happens at
Enander:2017ogx

$y_{\rm ta}\approx\frac{0.7}{\delta_{i}}.$ (5.7)

Without a pressure term, the collapse will continue without a restriction, and
the overdensity will eventually collapse to a single point at $y=y_{c}$ where
the overdensity $\delta$ becomes infinite. Again, by numerical solution we
observed that

$y_{c}\approx\frac{1.1}{\delta_{i}}.$ (5.8)

Numerical calculation of $\delta_{i}y_{\rm ta}$, and $\delta_{i}y_{c}$ for a
bunch of initial overdensities is shown on the left plot of Figure 13.

Figure 13: Left: Plot of $\delta_{i}y_{\rm ta}$ and $\delta_{i}y_{c}$ for
various values of the initial density $\delta_{i}$. The solid lines denote the
numerical results, while dashed lines are the mean values for the initial
densities considered in the plot. Right: The critical density $\delta_{c}(x)$
predicted by the linear theory calculated numerically via (5.13) together with
approximate expressions in the deep matter (dashed orange line), and deep
radiation eta (dashed green line).

After having calculated the time of collapse using the _non-linear_ theory, we
can ask the following question: What would be the prediction of the _linear_
theory if we had extrapolated it until the collapse $y_{c}$? Ellis:2020gtq
This is the necessary ingredient for using the Press-Schechter formalism which
we will review in the next section. The conservation of the mass $M$ implies

$\quantity[1+\delta(y)]b^{3}(y)=1+\delta_{i}.$ (5.9)

In the linear approximation we define $\delta(y)=\delta_{i}+\Delta\delta(y)$,
and evaluate all expressions at linear order in $\Delta\delta(y)$. Then (5.9)
implies $b(y)=1-\Delta\delta(y)/3(1+\delta_{i})$, and (5.5) becomes

$y(1+y)\derivative[2]{\Delta\delta}{y}+\quantity(1+\frac{3}{2}y)\derivative{\Delta\delta}{y}-\frac{3}{2}\quantity[\frac{2}{3}+\frac{1}{3(1+\delta_{i})}]\Delta\delta=-\frac{3}{2}\delta_{i}.$
(5.10)

Note that in the $\delta_{i}\ll 1$ limit, this equation reduces to the
_Meszaros_ equation as expected. The solution of this equation with the
initial condition $\Delta\delta(y\rightarrow 0)=0$ is

$\Delta\delta(y)=\frac{3\delta_{i}(1+\delta_{i})}{3+2\delta_{i}}\quantity[{}_{2}F_{1}\quantity(\delta_{-},\delta_{+},1,-y)-1],\quad\delta_{\pm}\equiv\frac{1}{4}\pm\frac{\sqrt{17\delta_{i}^{2}+42\delta_{i}+25}}{4+4\delta_{i}},$
(5.11)

where ${}_{2}F_{1}$ is the hypergeometric function. So in the linear theory,
$\delta(y)$ is given by

$\delta^{\rm
lin}(y)=\delta_{i}\quantity[\frac{3+3\delta_{i}}{3+2\delta_{i}}\,_{2}F_{1}\quantity(\delta_{-},\delta_{+},1,-y)-\frac{\delta_{i}}{3+2\delta_{i}}].$
(5.12)

Since the non-linear theory predicts that the initial fluctuation $\delta_{i}$
will collapse at $y_{c}$, we can calculate the prediction of the linear theory
for the overdensity at collapse by replacing $\delta_{i}$ with the
corresponding $y_{c}$, and setting $y=y_{c}$ in (5.12). So the prediction of
the linear theory at time of the collapse is
$\delta_{c}(y_{c})=\evaluated{\delta_{i}\quantity[\frac{3+3\delta_{i}}{3+2\delta_{i}}\,_{2}F_{1}\quantity(\delta_{-},\delta_{+},1,-y_{c})-\frac{\delta_{i}}{3+2\delta_{i}}]}_{\delta_{i}\approx
1.1/y_{c}}.$ (5.13) The quantity $\delta_{c}$ is called the _critical density
at collapse_. Small fluctuations $\delta_{i}\ll 1$ collapse during the deep
matter era where $y_{c}\gg 1$. In this case we recover the well-known result:

$\delta_{c}(y_{c}\gg
1)\approx\frac{1.1}{y_{c}}\quantity(1+\frac{3}{2}y_{c})\approx\frac{3}{20}\quantity(12\pi)^{2/3}\approx
1.686,\quad\delta_{i}\ll 1.$ (5.14)

On the other hand, for large fluctuations we can consider the $y_{c}\ll 1$
limit to approximate the hypergeometric function at $1+y_{c}$. Then we find

$\delta_{c}(y_{c}\ll
1)\approx\delta_{i}\quantity(1+\frac{3}{2}y_{c})\approx\frac{1.1}{y_{c}}\quantity(1+\frac{3}{2}y_{c})\approx\frac{1.1}{y_{c}},\quad\delta_{i}\gg
1.$ (5.15)

Suprisingly, the approximate expression for small and large fluctuations are
the same. So for all regimes, one can use
$\delta_{c}\approx(1.1/x_{c})(1+3x_{c}/2)$ as the expression for the critical
density at collapse Ellis:2020gtq . The plot of the critical density
$\delta_{c}$ as a function of $y_{c}$ is plotted on the right plot of Figure
13.

###  5.2 Press-Schechter Formalism

The formation of the dark matter halos can be studied analytically via the
Press-Schechter (PS) formalism Press:1973iz . It relies on the following
postulate: The fraction of matter $\mathscr{F}$ which is inside collapsed
structures of comoving size larger than $R$ at any given scale factor $a$ is
the same as the probability $\mathscr{P}$ of finding an overdensity
$\delta_{R}(\mathbf{x},a)>\delta_{c}(a)$ where $\delta_{R}(\mathbf{x},a)$ is
the overdensity _smoothed_ at the scale $R$, and $\delta_{c}(a)$ is the
critical density for collapse at scale factor $a$. If the initial fluctuations
prior to collapse are small, then there is a one-to-one correspondence between
the mass $M$ of the halo and the comoving radius $R$. Then the PS postulate
can be expressed as

$\mathscr{F}(>M;a)=\mathscr{P}(\delta_{R}(a)>\delta_{c}(a)).$ (5.16)

The relation between $M$ and $R$ depends on the choice of the _window
function_ $W_{R}$ which is used to smooth out the density contrast. For any
window function, the smoothed density contrast is given by

$\delta_{R}(\mathbf{x},a)=\int\differential[3]{\mathbf{x}^{\prime}}W_{R}(\absolutevalue{\mathbf{x}-\mathbf{x^{\prime}}})\delta(\mathbf{x}^{\prime},a),\quad\int\differential[3]{\mathbf{x}}W_{R}(\mathbf{x})=1,$
(5.17)

where $\delta$ is the energy density contrast. By a straightforward
calculation we can obtain the variance of the smoothed density contrast as

$\sigma_{R}^{2}(a)=\expectationvalue{\absolutevalue{\delta_{R}(\mathbf{x},a)}^{2}}=\int_{0}^{\infty}\frac{\differential{k}}{k}\absolutevalue{\widetilde{W}_{R}(k)}^{2}\mathcal{P}_{\delta}(k;a),$
(5.18)

where $\widetilde{W}$ is the Fourier transform of the window function and
$\mathcal{P}_{\delta}$ is the dimensionless power spectrum corresponding to
the density contrast.

The most common and natural choice for the window function is the _spherical
top-hat (STH)_ window function which is defined by

$W_{\rm
STH}(\mathbf{x})=\quantity(\frac{4\pi}{3}R^{3})^{-1}\times\begin{cases}1,&\mathbf{x}\leq
R\\\ 0,&\mathbf{x}>R\end{cases}.$ (5.19)

Its Fourier transform is

$\widetilde{W}_{\rm
STH}(kR)=\frac{3}{(kR)^{3}}\quantity[\sin(kR)-kR\cos(kR)].$ (5.20)

With this window function, the smoothed density contrast
$\delta_{R}(\mathbf{x})$ is just the integral of the density contrast over a
spherical region of comoving size $R$ centered at $\mathbf{x}$. Then the mass
enclosed in the region is unambiguously given by

$M(R)=\frac{4\pi}{3}\overline{\rho}a^{3}R^{3}\approx\frac{4\pi}{3}\overline{\rho}_{m,0}a_{0}^{3}R^{3},$
(5.21)

where the approximation holds for small fluctuations. We will use this
relation even if the fluctuations are large151515The Press-Schechter Formalism
relies on the assumption that there is a one-to-one correspondence between the
mass $M$ of the halo and the comoving radius $R$ which breaks down if the
fluctuations are large. A modified version of the formalism which does not
rely on this one-to-one correspondence is proposed in Enander:2017ogx .
However, we found that this formalism is not very useful to construct the halo
spectrum, which is our main goal in this section. Therefore, we will present
our estimates using the original Press-Schechter formalism..

Note that in order to calculate the variance via (5.18), one needs to perform
an integration over all momentum values, not just the ones we have simulated
in the previous sections. The contribution of higher modes can be neglected
since the power spectrum is suppressed, however the contribution of the lower
momentum modes should be included. Since we have shown that the power spectrum
converges to the CDM result for low momentum modes, we calculate the variance
as161616In the case of complete fragmentation we took
$\mathcal{P}_{\delta}=\mathcal{P}_{\delta}^{\rm quad}$ in the second term of
(5.22) since for all the simulated modes the linear contribution is
negligible.

$\sigma_{R}^{2}(a)=\int_{0}^{k_{\rm
min}}\frac{\differential{k}}{k}\absolutevalue{\widetilde{W}_{R}(k)}^{2}\mathcal{P}_{\delta,\rm
CDM}(k;a)+\int_{k_{\rm min}}^{k_{\rm
max}}\frac{\differential{k}}{k}\absolutevalue{\widetilde{W}_{R}(k)}^{2}\mathcal{P}_{\delta}(k;a),$
(5.22)

where $k_{\rm min}$ and $k_{\rm max}$ are the minimum and the maximum modes
which are simulated. To calculate the CDM contribution during the radiation
era, we have used the _gauge-invariant_ density contrast $\delta_{\rm
CDM}^{\rm GI}$ given by Zhang:2017flu

$\delta_{\rm CDM}^{\rm GI}(k;a)=9\Phi_{k}(0)\quantity[\frac{\sin
t_{k}}{t_{k}}+\frac{\cos t_{k}-1}{t_{k}^{2}}-\operatorname{Ci}(t_{k})+\ln
t_{k}+\gamma_{E}-\frac{1}{2}],\quad t_{k}=\frac{k/a}{\sqrt{3}H}$ (5.23)

so that the dimensionless power spectrum becomes171717The difference between
the density contrast calculated in the Newtonian gauge and the gauge-invariant
one decays rapidly on sub-horizon scales. However, on super-horizon scales the
Newtonian one becomes a constant while the gauge-invariant one scales with
$t_{k}^{2}$. We are using the gauge-invariant one in this calculation in order
to have a well-defined integral for the variance.

$\boxed{\mathcal{P}_{\delta,\rm CDM}^{\rm rad}(k;a)=\frac{4}{9}A_{s}\times
81\quantity[\frac{\sin t_{k}}{t_{k}}+\frac{\cos
t_{k}-1}{t_{k}^{2}}-\operatorname{Ci}(t_{k})+\ln
t_{k}+\gamma_{E}-\frac{1}{2}]^{2}.}$ (5.24)

During the matter era, we use

$\boxed{\mathcal{P}_{\delta,\rm CDM}^{\rm
mat}(k;z)=A_{s}\quantity(\frac{k_{\rm
eq}^{2}}{5\Omega_{m}H_{0}^{2}})^{2}\frac{K^{4}\mathcal{T}_{\rm
BBKS}^{2}(K)}{(1+z)^{2}},}$ (5.25)

where $\Omega_{m}\approx 0.32$, $K\equiv\sqrt{2}k/k_{\rm eq}$ with $k_{\rm
eq}$ is the wavenumber which enters the horizon at matter-radiation equality,
and $\mathcal{T}_{\rm BBKS}$ is the BBKS transfer function Bardeen:1985tr :

$\mathcal{T}_{\rm
BBKS}(K)=\frac{\ln(1+0.121K)}{0.121K}\quantity[1+0.201K+\quantity(0.831K)^{2}+\quantity(0.282K)^{3}+\quantity(0.346K)^{4}]^{-1/4},$
(5.26)

We review the derivation of (5.25) in Appendix B.

Figure 14: Evolution of the variance $\sigma_{R}$ for the benchmark point
$m=10^{-15}\,$\mathrm{eV}$$ and $f=1.3\times 10^{14}\,$\mathrm{GeV}$$ shown
with solid lines. The red and blue lines denote redshifts in radiation and
matter eras respectively. The dashed lines show the corresponding spectra for
CDM. We see that at large scales, i.e. large comoving radius, the variance
converges to the CDM prediction at all times. During the radiation era, the
variance grows logarithmically on large scales like CDM, while at small scales
it is frozen since the variance integral (5.22) is dominated by the modes that
are above the Jeans scale. Slightly after the matter-radiation equality, all
relevant modes drop below the Jeans scale, so the variance grows linearly with
the scale factor at all scales. We also show a horizontal gray line at which
the variance becomes equal to the critical density at collapse during the
matter era.

In Figure 14, we show the plots of the variance calculated at different
redshifts for a benchmark ALP model with $m=10^{-15}\,$\mathrm{eV}$$ and
$f=1.3\times 10^{14}\,$\mathrm{GeV}$$ with solid lines, while the dashed lines
show the CDM predictions. The lines with red and blue color tones denote the
redshifts in radiation and matter eras respectively. We can observe that the
variance evolution differs between small scales (smaller comoving radius) and
large scales (larger comoving radius). On large scales, all lines do converge
to the CDM predictions which grow logarithmically/linearly with the scale
factor during the radiation/matter eras. On the other hand, the small scales
are frozen initially since the variance integral (5.22) at those scales are
dominated by the modes which are above the Jeans scale. Only slightly after
the matter-radiation equality these scales drop below the Jeans scale, and as
a result the variance grows linearly at all scales.

Figure 15: The variance of the density power spectum today (5.22) calculated
with the spherical top-hat window function (5.19) for different ALP benchmarks
with the mass $m=10^{-15}\,$\mathrm{eV}$$, and CDM. The curves with the red
tones denote the benchmarks with incomplete fragmentation (linear regime),
while the ones with the green tones show the benchmarks where the
fragmentation is complete (non-linear regime).

We compare the variance evaluated at today for different ALP benchmarks, and
CDM in Figure 15. We fixed the ALP mass to $m=10^{-15}\,$\mathrm{eV}$$, and
gradually decreased the ALP decay constant starting from a value which is
close to the Kinetic-Standard Misalignment boundary. The curves with the red
and green tones denote the benchmarks where the fragmentation is incomplete
and complete respectively, corresponding to linear and non-linear regimes. In
the linear regime, decreasing the axion decay constant yields larger
enhancement of the small scales. The variance gets peaked for the benchmark
which is closest to the incomplete-complete fragmentation boundary. Decreasing
the ALP decay constant further reverses the trend, yielding to less
enhancement. This feature is closely related to the shape of the power
spectrum that we show in Figure 9. As the ALP decay constant decreases in the
non-linear regime, the power spectrum becomes more broader, and its peak is
also shifted to higher momentum modes. These modes experience less growth in
the matter era due to the Jeans suppression, and as a result the variance gets
suppressed.

###  5.3 Mass distribution of the halos

We now continue our discussion of the PS formalism in order to obtain the mass
distribution of the halos. We have seen that the relevant quantity to
calculate is $\mathscr{P}(\delta_{R}(a)>\delta_{c}(a))$. The easiest approach
is to assume that $\delta(\mathbf{x},a)$, hence $\delta_{R}(\mathbf{x},a)$ is
a Gaussian random field. Then the probability density function (PDF) for
$\delta_{R}$ is given by

$f(\delta_{R};a)=\frac{1}{\sqrt{2\pi\sigma_{R}^{2}(a)}}\exp(-\frac{\delta_{R}^{2}(a)}{2\sigma_{R}^{2}(a)}).$
(5.27)

Then by using the PS postulate we can write

$\mathscr{F}\quantity(>M;a)=\mathscr{P}(\delta_{R}(a)>\delta_{c}(a))=\int_{\delta_{c}}^{\infty}\differential\delta
f(\delta;a)=\frac{1}{2}\operatorname{erfc}\quantity(\frac{\delta_{c}(a)}{\sqrt{2\sigma_{R}^{2}(a)}}).$
(5.28)

In the case of CDM, the variance $\sigma_{R}^{2}$ goes to infinity as
$R\rightarrow 0$. This implies that all the matter in the universe must
finally be in virialized objects so we expect $\mathscr{F}(>M;a)\rightarrow 1$
as $M\rightarrow 0$ Maggiore:2009rv . However, the PS result (5.28) has the
limit $\mathscr{F}(>M;a)\rightarrow 1/2$ as $M\rightarrow 0$, meaning that the
PS formalism predicts that only half of the matter reside in collapsed
structures. The reason behind this disagreement is that in PS formalism only
overdense regions do collapse. In reality, underdense regions can also
collapse if they are enclosed within a larger overdense region mo2010 . This
effect can be included by using the _excursion set formalism_ Bond:1990iw ,
for a nice review see Zentner:2006vw . The correct expression is181818This
result is derived under the assumption that the variance is smoothed via a
_Fourier-space top-hat (FTH)_ window function $\widetilde{W}_{\rm
FTH}(kR)=\Theta(1-kR)$, with $\Theta$ being the Heaviside function. This is
simply the analog of the spherical top-hat in momentum space. The major
drawback of this window function is that one cannot define a mass for the
overdensity as we have done in (5.21) since its volume is not well-defined
Maggiore:2009rv . The reason for choosing the FTH window function is that only
with this choice the excursion set problem has an analytical solution
Zentner:2006vw .

$\mathscr{F}\quantity(>M;a)=\operatorname{erfc}\quantity(\frac{\delta_{c}(a)}{\sqrt{2\sigma_{R}^{2}(a)}}),$
(5.29)

which has the correct limit as $M\rightarrow 0$.

With this final result we can write the comoving number density of collapsed
objects within the mass range $(M,M+\differential{M})$ as

$\begin{split}\derivative{n(M;a)}{M}\differential{M}&=\frac{\overline{\rho}_{m,0}}{M}\absolutevalue{\derivative{\mathscr{F}(>M;a)}{M}}\differential{M}\\\
&=\sqrt{\frac{2}{\pi}}\frac{\overline{\rho}_{m,0}}{M^{2}}\frac{\delta_{c}(a)}{\sigma_{M}(a)}\absolutevalue{\derivative{\ln\sigma_{M}(a)}{\ln
M}}\exp\quantity(-\frac{\delta_{c}^{2}(a)}{2\sigma_{M}^{2}(a)})\differential{M}.\end{split}$
(5.30)

This is known as the _Press-Schechter mass function_. In the literature, one
commonly uses the so-called _halo mass function (HMF)_ which gives the number
density of halos per logarithmic mass bin:

$\boxed{\derivative{n(M;a)}{\ln
M}=\frac{1}{2}\frac{\overline{\rho}_{m,0}}{M}f_{\rm
PS}(\nu)\absolutevalue{\derivative{\ln\sigma_{M}^{2}(a)}{\ln M}},}$ (5.31)

where $\nu=\delta_{c}(a)/\sigma_{M}(a)$, and

$f_{\rm PS}(\nu)=\sqrt{\frac{2}{\pi}}\,\nu\exp(-\frac{\nu^{2}}{2}).$ (5.32)

The advantage of using this form is that one can replace $f_{\rm PS}(\nu)$ by
other fit functions so that the HMF agrees more accurately with the results of
the N-body simulations. An often adopted fit function is the Sheth-Tormen fit
Sheth:1999su which uses an ellipsoidal collapse model rather than a spherical
one as in Press-Schechter. However, this fit function contains numerical
factors which are obtained by solving the excursion set problem explicitly for
CDM. So the numerical factors need to be recalculated for each model
Du:2016zcv . In this work we will use the PS mass function (5.32) since it is
more general and simpler.

Figure 16: Left plot: The evolution of the halo mass function (5.31) for the
benchmark point with $m=10^{-15}\,$\mathrm{eV}$$ and $f=1.3\times
10^{14}\,$\mathrm{GeV}$$. We see that the number density of lighter halos
decrease with redshift, while the number density of heavier ones increase.
This is a consequence of merging of smaller halos with time forming larger
halos. Right plot:  The halo mass function today for different benchmark
points and CDM. The ALP mass is set to $m=10^{-15}\,$\mathrm{eV}$$. The color
codes are the same as in Figure 15.

We show our results about the halo mass function in Figure 16. On the left
plot we show the time evolution of the HMF for the benchmark point
$m=10^{-15}\,$\mathrm{eV}$$ and $f=1.3\times 10^{14}\,$\mathrm{GeV}$$. We see
that the structure formation starts slightly after the matter-radiation
equality with light halos, and with time heavier halos begin to form. The
number density of light halos increase initially, however they start to
decrease after some time. Meanwhile, the number density of heavy halos is
always increasing meaning that the number density shifts from lighter halos to
heavier halos. This feature arises because light halos merge with each other
forming heavier halos.

The halo mass function today for different benchmark points is presented on
the right plot of Figure 16. The number density of heavy halos are not
affected by the fragmentation, while the number density of light halos are
suppressed compared to CDM due to the existence of the Jeans scale. On
intermediate scales, there is an enhancement compared to CDM due to the
fragmentation

Figure 17: Analog of Figure 16, except that the vertical axis is in terms of
the dimensionless halo mass function defined in (5.33). This dimensionless
quantity gives roughly the fraction of dark matter that resides in halos with
masses in the logarithmic mass bin $\quantity[\ln M,\ln M+\differential{\ln
M}]$. The stars on the right plot show the local maxima of $X(M)$ which
approximately gives the mass of the halos that have the largest fraction of
dark matter, excluding very heavy halos.

We also present our results of the halo mass function in terms of the
dimensionless variable Enander:2017ogx

$\boxed{X(M)\equiv\frac{M}{\overline{\rho}_{m,0}}\derivative{n(M;a)}{\ln M}.}$
(5.33)

in Figure 17. This approximately gives the fraction of dark matter that reside
in halos with masses between $\quantity[\ln M,\ln M+\differential{\ln M}]$. We
also denote the local maxima of these curves with stars on the right plot.
Note that for large $M$, the curves for all benchmarks converge to the CDM
prediction for which $X(M)$ grows with $M$191919The CDM prediction for $X(M)$
has a global maximum around $M\sim 10^{14}M_{\odot}$, see Figure 18..
Therefore heavier halos contain a larger fraction of dark matter even with
fragmentation. However, these halos collapse very late, so they are not
sufficiently dense to be interesting for observational prospects as we shall
see in the next section.

|  $m=10^{-5}\,$\mathrm{eV}$$ red, left |  $m=10^{-10}\,$\mathrm{eV}$$ blue, middle |  $m=10^{-15}\,$\mathrm{eV}$$ green, right
---|---|---|---
KMM solid | $f=4.0\times 10^{11}\,$\mathrm{GeV}$$ | $f=6.7\times 10^{12}\,$\mathrm{GeV}$$ | $f=1.2\times 10^{14}\,$\mathrm{GeV}$$
LMM dot-dashed |  $f=4.0\times 10^{11}\,$\mathrm{GeV}$$ $\absolutevalue{\pi-\Theta_{i}}\approx 1.2\times 10^{-9}$ |  $f=6.7\times 10^{12}\,$\mathrm{GeV}$$ $\absolutevalue{\pi-\Theta_{i}}\approx 7.0\times 10^{-9}$ |  $f=1.2\times 10^{14}\,$\mathrm{GeV}$$ $\absolutevalue{\pi-\Theta_{i}}\approx 3.7\times 10^{-7}$
Post-inf dotted | $f=1.5\times 10^{12}\,$\mathrm{GeV}$$ | $f=2.3\times 10^{13}\,$\mathrm{GeV}$$ | $f=3.5\times 10^{14}\,$\mathrm{GeV}$$
SMM dashed |  $f=3.0\times 10^{12}\,$\mathrm{GeV}$$ $\absolutevalue{\pi-\Theta_{i}}\approx 0.6$ |  $f=4.0\times 10^{13}\,$\mathrm{GeV}$$ $\absolutevalue{\pi-\Theta_{i}}\approx 0.4$ |  $f=6.0\times 10^{14}\,$\mathrm{GeV}$$ $\absolutevalue{\pi-\Theta_{i}}\approx 0.4$
Table 1: The list of benchmarks used in constructing the halo mass functions
in Figure 18 and the halo spectra shown in Figure 20.  Figure 18:
Dimensionless halo mass functions predicted by the production mechanisms
Kinetic Misalignment (solid), Large Misalignment (dot-dashed), post-
inflationary scenario (dotted), and Standard Misalignment (dashed). Different
colors denote different ALP masses; $m=10^{-5}\,$\mathrm{eV}$$ (red, left),
$m=10^{-10}\,$\mathrm{eV}$$ (blue, middle), $m=10^{-15}\,$\mathrm{eV}$$
(green, right). The ALP decay constants are given in Table 1. The CDM
prediction is shown via the brown line, and the curves do converge to the CDM
prediction for larger halo masses.

In Figure 18, we show a comparison of the dimensionless halo mass functions
corresponding to different ALP masses and production mechanisms. The masses
are $m=10^{-5}\,$\mathrm{eV}$$ (red, left), $m=10^{-10}\,$\mathrm{eV}$$ (blue,
middle), $m=10^{-15}\,$\mathrm{eV}$$ (green, right). Different linestyles
denote different production mechanisms: Kinetic Misalignment (solid), Large
Misalignment (dash-dotted), post-inflationary scenario (dotted), and Standard
Misalignment (dashed). When comparing Kinetic and Large Misalignment
mechanisms we have used the same value for the axion decay constant. For the
post-inflationary scenario we set the decay constant such that the correct
relic density is obtained. The list of benchmarks that we used to produce this
plot are tabulated in Table 1. The CDM prediction is shown by the brown line,
and all the curves converge to it at large halo masses.

We end this section with a short discussion of the low-mass behavior of the
halo mass function. As we can see from the right plot of Figure 16, even
though the halo mass function for ALP is suppressed compared to CDM it still
continues to grow with decreasing mass. However, this is not physical since
the existence of the Jeans scale should heavily suppress the structure
formation below a certain mass. This is a well-known problem of the Press-
Schechter formalism applied to ALP dark matter. In Arvanitaki:2019rax and
Kulkarni:2020pnb , it has been argued that this isssue stems from the fact
that even for a small comoving radius $R$, the spherical top-hat window
function still collects contributions from larger scales, i.e. small
$k$-modes. Other works claim that the issue is not the spherical top-hat
window function, but the critical density for collapse that we have computed
in Section 5.1 should be modified to account for the Jeans scale Marsh:2013ywa
; Bozek:2014uqa ; Du:2016zcv . We will not investigate this issue further in
this work, and continue our calculations by using the spherical top-hat window
function, and scale-independent critical density.

###  5.4 Halo spectrum

In this section we derive the halo spectrum predicted by the ALP
fragmentation, compare it with the predictions of Standard Misalignment and
the post-inflationary scenario. For this we need to know the density profile
of the halos. We start by discussing the profiles of the CDM halos, and later
comment on the modifications due to the wavelike nature of ALPs.

#### CDM Halos

The N-body simulations of CDM show that for a halo identified, or measured, at
redshift $z$, the density profile is well-approximated by the _Navarro-Frenk-
White (NFW)_ profile given by Navarro:1996gj

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

where $r_{s}$ is called the _scale radius_ which is defined by the condition

$\evaluated{\derivative{\ln\rho(r)}{\ln r}}_{r_{s}}=-2,$ (5.35)

$\delta_{\rm char}$ is the dimensionless _characteristic density_ , and
$\rho_{c}(z)$ is the critical density of the universe at redshift $z$. An
alternative but equivalent definition is

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

where $\rho_{s}=\rho(r=r_{s})$ is called the _scale density_. From this
defition, one can also define the _scale mass_ $M_{s}$ as the mass enclosed
inside the scale radius:

$M_{s}=4\pi\int_{0}^{r_{s}}\differential{r}r^{2}\rho(r)=16\pi\rho_{s}r_{s}^{3}f(1),\quad
f(x)\equiv\ln(1+x)-\frac{x}{1+x}.$ (5.37)

It is well known that the energy density inside a virialized halo is
approximately $200$ times the mean energy density when the halo was formed.
This motivates the definition of a quantity called the _virialized radius_
$r_{200}$ such that the mean energy density enclosed in a spherical region of
radius $r_{200}$ is equal to $200$ times the critical density:

$200\rho_{c}=\quantity(\frac{4\pi}{3}r_{200}^{3})^{-1}4\pi\int_{0}^{r_{200}}\differential{r}r^{2}\rho(r).$
(5.38)

One also defines the _concentration parameter_ $c_{200}$ by

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

Then we can use (5.38) to relate the characteristic density $\delta_{\rm
char}$ to the concentration parameter $c_{200}$:

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

Even though we have not written down the redshift dependence explicitly, both
$\delta_{\rm char}$ and $c_{200}$ are redshift-dependent.

The quantities $\delta_{\rm char}$ and $c_{200}$ need to be determined by the
N-body simulations. However there is an emprical method originally developed
by NFW for CDM Navarro:1996gj , which can also be used for axion miniclusters
Ellis:2020gtq . Consider a halo identified, or measured, at redshift $z=z_{\rm
ind}$. Assign this halo a mass $M_{200}$ which is equal to the mass enclosed
in the spherical region of radius $r_{200}$. By definition we have

$M_{200}=\frac{4\pi}{3}\quantity(200\times\rho_{c}(z_{\rm ind}))r_{200}^{3}.$
(5.41)

This mass can be identified with the mass $M$ of the collapsing region defined
in (5.21), i.e. $M\leftrightarrow M_{200}$202020The reasoning behing this
identification is the assumption that the mass $M$ inside the collapsing
region as defined by (5.21) remains approximately constant during the collapse
and the virialization Blinov:2021axd .. Then assign this halo a _collapse
redshift_ $z_{\rm col}(M_{200},\mathfrak{f})$212121The parameter $z_{\rm col}$
has nothing to do with the redshift of the linear collapse $y_{c}=(1+z_{\rm
eq})/(1+z_{\rm c})$ that we have introduced in Section 5.1.  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}$. This
can be estimated using the Press-Schechter formalism: Lacey:1993iv ;
Lacey:1994su

$\boxed{\operatorname{erfc}\quantity[X(z_{\rm col})-X(z_{\rm
ind})]=\frac{1}{2},}$ (5.42)

where

$\boxed{X(z)=\frac{\delta_{c}(z)}{\sqrt{2\quantity[\sigma^{2}(\mathfrak{f}M_{\rm
200};z)-\sigma^{2}(M_{\rm 200};z)]}}.}$ (5.43)

In this equation, $\delta_{c}(z)$ is the critical density for gravitational
collapse, and $\sigma^{2}(M_{200};z)$ is the variance at the redshift $z$
evaluated with a spherical tophat window function with the smoothing scale

$R=\quantity(\frac{3M_{\rm 200}}{4\pi\rho_{m,0}a_{0}^{3}})^{1/3},$ (5.44)

where $\rho_{m,0}$ is the CDM/axion density today222222We are assuming that
CDM/axion make all of matter., and $a_{0}$ is the scale factor today.

Now assume that $\delta_{\rm char}(z_{\rm ind})\rho_{c}(z_{\rm ind})$, and
hence $\rho_{s}$, is proportional to the energy density of the universe at the
time of collapse:

$\delta_{\rm char}(z_{\rm ind})\rho_{c}(z_{\rm ind})=4\rho_{s}(z_{\rm
ind})=\mathfrak{C}(z_{\rm col})\rho_{c}(z_{\rm col}),$ (5.45)

where $\mathfrak{C}(z_{\rm col})$ is in general $z_{\rm col}$-dependent
proportionality factor. Both $\mathfrak{C}$ and $\mathfrak{f}$ are free
parameters which need to be fit with the N-body simulations. For both CDM
Navarro:1996gj and axion miniclusters Ellis:2020gtq , it is found that
$\mathfrak{f}=10^{-2}$ provides a good fit to the simulations. Hence, we will
also use this value in our calculations. To determine $\mathfrak{C}(z_{\rm
col})$ we have used the CDM simulations compiled in Sanchez-Conde:2013yxa ,
and axion minicluster simulations in Eggemeier:2019khm . The details of the
fitting procedure is explained in the Appendix C. Once this fitting is done,
all the halo properties, most importantly the scale mass $M_{s}$ and the scale
density $\rho_{s}$, can be calculated.

#### ALP Halos

At small scales, the halo density profile for ALP dark matter cannot be
described by the NFW profile due to the wavelike nature of ALPs. For recent
reviews, see for example Marsh:2015xka and Niemeyer:2019aqm . Numerical
simulations show that the central profiles of ALP halos do not obey the
$r^{-1}$ scaling, the so-called “cusp”, as in the NFW halos Hu:2000ke ;
Schive:2014dra ; Schive:2014hza . Instead, the density profiles are described
by _solitons_ 232323Precisely, these are _pseudo-solitons_ since the ALP field
is time-dependent, and they are not absolutely stable Marsh:2015xka . True
soliton solutions require a conserved Noether charge Diez-Tejedor:2013sza .
Therefore they occur for complex scalar field dark matter where the conserved
$U(1)$ charge guarantees stability. These are known as the _boson stars_ in
the literature Liddle:1992fmk . which are the eigenstates of the time-
independent Schrödinger-Poisson equations PhysRev.187.1767 ;
PhysRevLett.66.1659 . The soliton profile can be fitted by Schive:2014dra

$\rho_{\rm sol}(r)=\rho_{\rm
sol}(0)\quantity[1+\quantity(2^{1/8}-1)\quantity(\frac{r}{r_{1/2}})^{2}]^{-8},$
(5.46)

where

$\rho_{\rm sol}(0)\approx
1.9\times\quantity(\frac{m}{10^{-17}\,$\mathrm{eV}$})^{-2}\quantity(\frac{r_{1/2}}{\rm{pc}})^{-4}M_{\odot}\rm{pc}^{-3}$
(5.47)

is the density at the core, and $r_{1/2}$ is the radius at which the density
is half of the core density. By using this form of the profile one can
calculate all the scale parameters $r_{s},\rho_{s}$ and $M_{s}$. Then one can
show that the scale density and the scale mass are related to each other by

$\boxed{\rho_{s}^{\rm sol}\approx 1.2\times
10^{-4}\quantity(\frac{m}{10^{-17}\,$\mathrm{eV}$})^{6}\quantity(\frac{M_{s}^{\rm
sol}}{M_{\odot}})^{4}\,M_{\odot}\rm{pc}^{-3}.}$ (5.48)

This ground state configuration is independent of the power spectrum, and is
reached rapidly by a process called _gravitational cooling_ Seidel:1993zk ;
Guzman:2006yc .

#### Results for the halo spectrum

Far away from the center of the halo, the density profile of the ALP halos can
again be explained by the NFW profile. Therefore, a realistic description
requires a profile which smoothly interpolates between the soliton profile
(5.46) at small scales, and the NFW profile (5.36) on large scales. A precise
way to obtain this matching is currently not known. In this work, we
calculated the halo spectrum by assuming the NFW profile at all scales,
however used the $M_{s}^{\rm sol}$–$\rho_{s}^{\rm sol}$ relation (5.48) as an
upper bound on the scale density for a given scale mass.

Figure 19: Halo spectrum for some benchmark points predicted by kinetic
fragmentation assuming the NFW profile for the halos. The ALP mass is assumed
to be constant and is fixed to $m=10^{-15}\,$\mathrm{eV}$$. Labels on top of
the curves denote the ALP decay constants. The gray line labeled as the
“soliton profile” is the spectrum if all ALP halos have the soliton profile
(5.46), and it is given by (5.48). The curves with the red tones denote the
benchmarks with incomplete fragmentation (linear regime), while the ones with
the green tones are the benchmarks with complete fragmentation (non-linear
regime). The dashed lines on the left of the soliton line are the halo spectra
if the ALP halos are described by the NFW profile at all scales. The stars
denote the halos for which $X(M)$ as defined by (5.33) has a local maxima,
meaning that the probability of these halos is larger compared to nearby
halos. Above the light/dark thin blue lines, weak gravitational lensing of
compact halos can induce localized distortions in the proper motion $\mu$ of
the background sources that can be observed by SKA/GAIA. Above the thin pink
line, the weak lensing can cause correlations in stellar proper accelerations
$\alpha$ that can be probable by Theia. In the region enclosed by the dark
yellow line, the compact halos can cause distortions in the waveform of the
gravitational waves emitted from a merger of two black holes that can be
observed with aLIGO. See Section 6.1 for more details.

We show the results of the halo spectrum in Figure 19. We assumed a constant
ALP mass, fixed it to $m=10^{-15}\,$\mathrm{eV}$$, and varied the ALP decay
constant. The color codes of the curves are the same as in Figure 15; the ones
with the red tones denote the benchmarks with incomplete fragmentation, and
the ones with the green tones are the cases with complete fragmentation. By
comparing the results with the Figure 19 we can clearly see that there is a
direct relationship between the variance today and the peak scale density. The
benchmark which has the largest variance today, is also the benchmark which
has the largest scale density. This is a direct consequence of the fact that a
larger variance implies larger fluctuations, and they collapse earlier
yielding to denser halos. We also show the soliton relation (5.48) by the gray
line on the left. The dashed lines on the left of it are the spectra if the
ALP halos have the NFW profile at all scales.

Finally, we compare the halo spectra corresponding to the different ALP masses
and production mechanisms in Figure 20. The ALP masses are
$m=10^{-5}\,$\mathrm{eV}$$ (red, left), $m=10^{-10}\,$\mathrm{eV}$$ (blue,
middle), and $m=10^{-15}\,$\mathrm{eV}$$ (green, right). For all masses we
compare the production mechanisms Kinetic misalignment with fragmentation
(solid lines), Large misalignment (dash-dotted lines), post-inflationary
scenario (dotted lines), and Standard misalignment (dashed lines)242424Kinetic
misalingment without fragmentation line is indistinguishable from the Standard
misalignment line so we don’t show in the plot. . The straight faint lines
labeled via the ALP mass show the soliton spectrum corresponding to the given
ALP mass. When comparing Kinetic and Large misalignment mechanisms we have
used the same value for the ALP decay constant. For the post-inflationary
scenario we set the decay constant such that ALPs make up all of the dark
matter. The list of benchmarks that we used when constructing the halo spectra
can be found in Table 1. We also show the region of the $M_{s}$–$\rho_{s}$
plane which can be probed by future experiments by thin lines. We see that
low-mass axions provide much more optimistic discovery prospects since the
halo spectra are peaked at larger masses.

Figure 20: The halo spectra corresponding to the benchmarks listed in Table 1
together with the regions observable by future lensing probes that we briefly
summarized in Section 6.1. Different colors show the different ALP masses;
$m=10^{-5}\,$\mathrm{eV}$$ (red, left), $m=10^{-10}\,$\mathrm{eV}$$ (blue,
middle), $m=10^{-15}\,$\mathrm{eV}$$ (green, right). Different linestyles show
different production mechanisms; Kinetic misalignment with fragmentation
(solid), Large misalignment (dot-dashed), post-inflationary scenario (dotted),
and Standard misalignment (dashed). The straight faint lines labeled via the
ALP mass show the soliton spectrum corresponding to the given ALP mass.

##  6 Observational prospects

In this section, we briefly comment on the phenomenological consequences of
the halo spectra that we derived in the previous section. In Section 6.1 we
discuss various experiments that have a potential to probe the halo spectrum
at small scales. In Section 6.2 we discuss the consequences of the compact ALP
halos for the terrestrial ALP detection experiments, such as holoscopes.

###  6.1 Probes of the halo spectrum

The compact dark matter halos which are denser than the CDM ones can be probed
by future gravitational surveys via their direct gravitational interactions. A
detailed study of the discovery prospects is beyond the scope of this paper.
Here we give a quick overview by using the sensitivity curves in the
$M_{s}$–$\rho_{s}$ plane presented in Arvanitaki:2019rax .

  * •

Weak gravitational lensing: When a compact DM halo passes near the line of
sight of a background light source such as a star or any other luminuous
object, it can induce apparent motions of these light sources through
gravitational lensing. In the case of _strong_ lensing, these distortions are
easily visible and form Einstein rings, arcs, and multiple images. If the
distortions are weak so that they can only be detected by analyzing a large
number of sources statistically, then one speaks of _weak_ gravitational
lensing. The discovery prospects of compact DM halos via weak gravitational
lensing has been studied in VanTilburg:2018ykj . The weak lensing of a compact
halo can induce localized distortions in the proper motion $\mu$ of the
background light sources that can be observed by GAIA 2018A&A…616A…1G and SKA
2004NewAR..48.1473F . We call this _$\mu$ lensing_, and show the observable
region via thin blue lines in Figures 19 and 20. Similarly, the weak lensing
can cause correlations in stellar proper accelerations $\alpha$ that can be