ALP dark matter mini-clusters from kinetic fragmentation

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.


Introduction
Axion-like-particles (ALPs) are prime candidates for dark matter, and they have been subject to extensive theoretical research in recent years (see [1][2][3][4] 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 Λ b f via their interactions with a strong gauge sector that breaks the continuous shift symmetry explicitly. This explicit breaking creates a periodic potential V (θ, T ) = Λ 4 b (T )[1 − cos(θ)] for the ALP field θ which leads to an ALP mass m(T ) ∝ Λ 2 b (T )/f . By far the best-known ALP model is the QCD axion, see [5] for a recent review, that arises as a by-product of the Peccei-Quinn solution to the Strong CP problem [6][7][8]. In this case, the strong sector is QCD so that Λ b ∼ Λ QCD , thus the QCD axion mass m 0 ≡ 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 [9][10][11]. 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 θ i = [−π, π) on each Hubble patch. The later evolution is governed byθ + 3H(t)θ + m 2 (t) sin θ = 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 θ = 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 [12,13]. 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 [14]. 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 [12,[15][16][17][18][19]. 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 θ = ±π, then move to the low-temperature minimum θ = 0 after some critical temperature T c [20,21]. A large kinetic energy can also be realized in models with more than one ALP fields [22,23].
In [14], 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 [24]. It was demonstrated in [14] 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 O(1) density fluctuations in the postinflationary scenario. These fluctuations gravitationally collapse very early, and form dense objects called ALP miniclusters [25][26][27][28]. In recent years, this process has been studied in detail both semi-analytically [29][30][31][32], and numerically [33][34][35][36][37][38][39][40][41]. For a recent review see [42]. A recent discussion of how the properties of the QCD axion miniclusters in the postinflationary scenario are modified in the presence of kinetic misalignment can be found in [43].
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 δθ H I /(2πf I ) where f I is the effective decay constant during inflation which might be different from its value today [44]. The ratio of the isocurvature fluctuations to the adiabatic ones are highly constrained by the Planck 2018 data [45,46]. 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 [47]. A well-known example related to the periodic potential is when the initial angle θ i is very close to the top of the potential 1 , i.e. |π − θ i | 1, referred as the large misalignment mechanism (LMM) [51] or extreme axion [52]. In this case, the onset of oscillations is delayed from the conventional condition m ∼ 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 [51][52][53][54][55][56]. This feature also arises in ALP models with non-periodic potentials for sufficiently large initial angles [57][58][59][60][61][62]. 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 [63].
The enhancement of fluctuations does also occur in KMM [14,24]. 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 ≈ 0 in the early time limit, and dark energy has no adiabatic perturbations [1]. On the other hand, the equation of state of the homogeneous mode in KMM is w ≈ 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 ∼ m osc a osc that get amplified most strongly. Here m osc and a 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 [64]. An important consequence of fragmentation is that the peak of the matter power spectrum can reach 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 [51]. 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 [65].
In [14] we demonstrated that fragmentation occurs both in temperature-dependent and independent potentials. In this work, we concentrate on ALP models with temperatureindependent 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.
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 [14] which are relevant for this work. The most important quantity for this paper is the trapping temperature T * 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.

Brief review of ALP fragmentation in kinetic misalignment
In this section, we summarize the main findings of [14]. We consider an ALP field θ with the Lagrangian In KMM, the evolution of the ALP homogeneous mode Θ can be divided into two regimes. When the ALP kinetic energy f 2Θ2 /2 dominates over the size of the potential barrier 2m 2 f 2 , then the ALP field rolls with the velocityΘ ∝ a −3 so that its energy density redshifts as ρ Θ ∝ 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 where s(T ) is the entropy density of the universe. Assuming entropy conservation, this quantity is conserved at early times when ρ Θ ∝ a −6 . The ALP relic density today can be expressed as where we took Ω DM = 0.12 [45]. 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 * : 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 [14] T 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 * , which can be obtained after solving for T * in (2.5). These regimes can be summarized as follows: 1. Standard misalignment (m/H * < 3): The onset of oscillations is not delayed from its conventional value m(T osc ) = 3H(T osc ), and there is no kinetic misalignment.
2. Kinetic misalignment with incomplete fragmentation (3 < m/H * 4 × 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. Complete fragmentation after trapping (4 × 10 1 m/H * 5 × 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 end < T * .
4. Complete fragmentation before trapping (5 × 10 2 m/H * 5 × 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 end > T * . 5. Non-perturbative (5 × 10 3 m/H * ): 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 [14]. We do this beacuse in this regime m/H * is not very large so the analytical method of [14] 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 nonperturbative analysis such as a lattice simulation, we will rely on the analytical theory in [14] 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 [30]. This is the topic of Section 4.

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.

Early evolution of fluctuations: before and during fragmentation
Our starting point is the Friedmann-Lemaitre-Robertson-Walker metric including curvature perturbations Φ, Ψ in the Newtonian gauge 2 : The anisotropic stress vanishes for scalar fields at the linear order in perturbation theory which sets Ψ = −Φ. At early times, we assume that the dominant component of the energy density is radiation, so Φ is an external field for the ALP.
We expand the ALP field as and keep only the terms that are linear in fluctuations δθ and Φ. Since the fluctuations are small, the average energy density is given by the energy density of the homogeneous mode: where denotes derivative with respect to conformal time. The density contrast is The homogeneous mode satisfies the following equation of motion in physical time t: Θ + 3HΘ + m 2 sin Θ = 0, (3.5) while the equation of motion for the mode functions θ k reads where Φ k are the Fourier modes of the curvature perturbations. Following [51], we introduce the following dimensionless quantities: where the approximate expression for t m is valid during radiation domination. The momentum variable 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 Φ k in radiation era: where Φ k (0) is a stochastic variable that depends on the amplitude of the scalar primordial fluctuations A s by Here R k is the comoving curvature perturbation, k = 0.05 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 × 10 −9 consistent with the Planck 2018 measurements (TT,TE,EE+lowE+lensing 68%) [45]. 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 [51,52,55] 3 : where we have used H = (2t) −1 . After choosing the correct initial conditions for the homogeneous mode Θ and the mode functions θ k , the evolution of the ALP field and its fluctuations can be determined by numerically solving (3.10) and (3.11).

Initial conditions
To fix the initial conditions for the homogeneous mode, we used the parameter that we have introduced in [14]: The homogeneous mode is rolling when > 1, and oscillating when < 1. In the former case, satisfies the following relation [14]: where E is the complete elliptic integral of the second kind 4 , and is an another dimensionless time variable defined such that τ = 1 coincides with the time of trapping. At any τ , this equation can be solved numerically to get (τ ). At early times when 1, we can approximate By plugging this result into (3.12) we get Here Θ i is the initial angle at the start of the simulation which is irrelevant in the case of kinetic misalignment mechanism 5 . We have chosen it as Θ i = π/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. ρ Θ ∝ a −6 . We describe the calculation in Appendix A.1, and in much more detail in [14]. For t k < 1, we can express the solution as a power series in t k as (see (A.7)) 4 We use the following definition for the elliptic integral: Note that most software packages such as Mathematica and scipy [67] use m = k 2 instead of k as their argument when defining elliptic integrals. 5 For some choices of tm,i and Θ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.   The conversion between t k and t m can be done via 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.

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 eV and f = 1.3 × 10 14 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, * = m/2H * . 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 [55] where Ci is the cosine integral, and γ E is the Euler-Mascheroni constant. Note that all modes approach to 3Φ k (0) in the super-horizon limit. This is becuase in this limit the adiabatic perturbations satisfy for all species i and j, where δ 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 δ r ≈ 2Φ 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 θ ≈ 1 before trapping gives δ ≈ 3Φ k (0) at early times in the super-horizon limit.
Later behavior can be interpreted as follows: • k = 0.1: This mode enters the horizon at t m = 3/4 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.
• 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.
• k = 1.7: This mode experiences significant growth due to the parametric resonance effect explained in [14], and its amplitude is enhanced significantly compared to CDM.

Power spectrum
By solving for the mode functions, we obtain the dimensionless power spectrum where δ k = δ k Φ 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 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.20) at t m = m/2H = 2 × 10 3 . For all benchmarks points, the ALP mass is m = 10 −15 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.

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 a * , i.e. t m 1, the ALP energy density redshifts like matter with an average equation of state w ≈ 0. Therefore, one can use the WKB approximation for both the homogeneous mode Θ, and the mode functions θ k (t). One makes the ansatz [68] where Θ ± are constants. With this approximation, one can directly derive the evolution of the density constrast as [68] δ k + 2Hδ k + c 2 s,eff where ρ r is the energy density in the radiation, δ r is the density constrast of the radiation, and c 2 s,eff is the effective sound speed of the ALP: 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 [69]. 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πGρ r term on the LHS is negligible during the matter era since c 2 s,eff ρ r ρ Θ 6 . Also, in radiation era we have 16πGρ r ≈ 6H 2 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 bÿ The terms 4πGρ Θ is responsible for the graviational collapse and the structure formation during the matter era, while c 2 s,eff 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 2 s,eff k 2 /a 2 4πGρ Θ , 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 [1] 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. Following [51], we define y ≡ a/a eq . (3.29) In terms of this variable (3.25) becomes 7 Therefore, the momentum mode k drops below the Jeans scale at y = 2 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 match m which is much later 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.
than the onset of oscillations t match m 1, but long before the matter-radiation equality y match 1. By using the fact that y = 2 1/4 H eq t m /m in radiation era, we find the matching conditions as Once it drops below the Jeans scale during the matter era, it starts to grow linearly. The evolution of the mode k = 5.0 is similar, except it does not experience exponential growth during trapping, and drops below the Jeans scale much later.

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 [51,52] is a special case of the Standard Misalignment Mechanism where the initial angle Θ i is very close to the top of the cosine potential, i.e. |π − Θ i | 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 mt as the dynamical variable. The initial conditions are given by If the initial angle is very close to the bottom of the potential |Θ i | 1, then the cosine potential can be approximated by a quadratic, which implies sin Θ ≈ Θ. In this case the homogeneous mode equation of motion (3.10) with initial conditions (3.32) has an exact solution: In this case, the energy density is given by In the late time limit t m 1, in other words much after the onset of oscillations, one can use asymptotic forms of the Bessel functions to show that so the ALP field scales as cold dark matter at late times. Redshifting this energy density until today we obtain where T m is the temperature at which t m = 1, i.e. 2H(T m ) = m, and g ρ (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 Θ 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 In most of the parameter space, the relative error is less than 10%. density today is given by where F anh is the anharmonicity correction which should behave as F anh (Θ i ) → 1 as Θ i → 0. By solving the equation of motion (3.5) numerically with different initial angles, we found that the fit function 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.
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 |π − Θ i | 1, the dimensionless time at which the oscillations start can be approximated by [51] t osc m = ln 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 [1]. 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 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: The results and comparison with KMM We first show the dimensionless power spectra at t m = m/2H = 2×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 k 1 behave as CDM, while larger modes k 1 are suppressed with respect to CDM. Only when the initial angle is significantly close to the top of the potential |π − Θ i | 10 −5 , there is an enhancement in the modes k ∼ O(1). The power spectrum reaches to O(1) values when |π − Θ i | 10 −10 , after which the linear perturbation theory becomes unreliable.
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 eV. We have chosen two benchmark values for the ALP decay constant; f = 6.0 × 10 14 GeV, and f = 1.3 × 10 14 GeV. The former is close to the boundary of the SMM-KMM transition in the language of Section 2 8 , 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 Θ i via (3.37) and (3.38) such that all of dark matter is made of ALPs. For f = 6.0 × 10 14 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 × 10 14 GeV (green lines) both models predict an enhancement of the modes with k ∼ O(1), however the enhancement in the case of KMM is much larger compared to LMM. 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 τ = 2H * t by definition. In terms of t m this is On the other hand, the onset of oscillations in LMM is given by (3.39). Then for f = 1.3 × 10 14 GeV one gets t osc m ≈ 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 where v is the velocity potential which for the ALP field takes the form [14] v Thus, in Fourier space the gauge-invariant density contrast is 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 In the SMM/LMM case, we can use (A.17) and (A.25) to obtain at leading order By roughly comparing the density contrast for a mode k ∼ 1 at t m ∼ 1 we find The enhancement in LMM happens when |π − Θ i | 1. In this case, on top of the two orders of magnitude numerical suppression, there is also a drastic suppression from the cosine term: 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 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 ω 2 k = k 2 /t m − cos Θ. During the period when the homogeneous mode is stuck at the top we have cos Θ ≈ cos Θ i ≈ −1. So the effective frequency of the mode becomes tachyonic when 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 Therefore, higher momentum modes are suppressed with 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.

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 Even though this approximation clearly breaks down during the fragmentation, it will eventually be a good one once the self-interactions become negligible 9 . 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 where ρ is the average energy density. For a quadratic potential (4.1) it takes the form (4.3) 9 One 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 [70][71][72]. For more information about the oscillons and their consequences for dark matter see [59] and the references therein.
In deriving this expression we have neglected the O(Φ 2 ) terms, but the ALP field θ(t, x) is kept non-perturbative. We now expand the ALP field as Hereθ k 's are stochastic variables that carry the statistical properties, while θ k 's are cnumber functions which depend only on the magnitude of the momentum k ≡ k. We assume thatθ k 's satisfy The function θ B describes the evolution of the long wavelength modes k < k * for which ∇θ ≈ 0. We took k * = a * H * , but this choice does not affect our results as long as k * is much smaller than the relevant modes for the fragmentation for which k ∼ ma * [14]. 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 [14]. This method does not work for the modes that are super-horizon during fragmentation, i.e. k * < a * H * . By plugging the expansion (4.4) into (4.3) we find By approximating 1 ± 2Φ ≈ 1, the density contrast reduces to 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 θ B → Θ, and ρ ≈ ρ Θ . Therefore, we refer this contribution to the density contrast as δ lin . Since the second term is quadratic in fluctuations, we refer it as δ quad , and write δ = δ lin + δ quad . On the scales k ∼ O(1) × ma * where the fragmentation is efficient, the quadratic contribution δ quad dominates over the linear contribution δ lin . However, at scales k ma * the quadratic contribution is small as we shall show soon, so the linear contribution dominates.
After this splitting, the power spectrum takes the form where 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 Φθθ whereΦ andθ 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θ ∼ Φ. Therefore, by assuming that the three-point function for the curvature perturbations vanishes ΦΦΦ = 0, the last term in (4.8) becomes zero. Then the power spectrum becomes where P lin The linear contribution is given by and it reduces to the power spectrum in the linear case when θ b → Θ, and ρ → ρ Θ . The quadratic contribution is derived in [30], and it reads The relevant contribution for the late time observables is the quadratic one. So we discuss it first, and then comment on the linear contribution.

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 [14]: where ϕ k is a momentum-dependent phase factor, P θ (k) is the initial power spectrum before the fragmentation and the redshift factor A k (t) is given by 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 [14]. By plugging these terms into the expression for the power spectrum (4.11) we obtain 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 time 10 . By averaging over the phase factors we can replace the cosine square term by 1/2. Also, without loss of generality, we can choose q to be aligned with theẑ-direction. Then |k − q| = k 2 + q 2 − 2kqu where u is the cosine of the angle between k and q. After integrating over the azimuthal angle we find It is more convenient to make a change of parameter from u to p ≡ |q − k|. This way we obtain The dimensionless power spectrum is Let q ≡ q/(a * m). Then, one can show that at small momenta where κ ≡ k/(a * 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 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 kin H 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.
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 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 κ ∼ 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.

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 ∼ 0.1 Mpc −1 . The relevant modes for fragmentation are those with k ∼ k peak ∼ a * m. To express this in dimensionless units, we use the following relation derived in [14]: where m * is the ALP mass at trapping, Λ b,0 = √ m 0 f is the today's value of the barrier height, and ρ Θ,0 is the ALP energy density today without taking fragmentation into account 11 . From (4.22) and assuming a constant mass we get (4.23) Therefore, the fragmentation scales are much shorter compared to the CMB ones 12 . 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 θ B is not known. However, from the discussion of the linear case in Section 3, we expect that the modes which are superhorizon 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: This approximation has been used in [73] to study the gravitational signatures of axion miniclusters when Peccei-Quinn breaking happens after inflation, and also in [58,74] where dark matter is generated from inflationary perturbations. 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 P 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 * 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 P quad until today. For this we will use the linearized evolution equation (3.27) even though 11 The fragmentation might slightly modify the relic density today. We made a semi-analytical estimate in [14], and found that this modification is not significant. 12 The 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 [14]. . 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 * for all benchmark points. 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.

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 θ takes random values θ ∈ [−π, π) 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 [75][76][77]. 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 [30] 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 [41]. In this section, first we calculate the power spectrum of a constant mass ALP model by adapting the method used in [30], and then make a comparison with the numerical results of [41]. Finally we compare the predictions in the post-inflationary scenario with the predictions of KMM and SMM/LMM.

Linearized analysis
Following [30], we separate the ALP mode functions as where the statistical properties and the initial conditions are encoded inθ 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 arë where m is the ALP mass which is assumed to be constant. The oscillations start around We also define H 1 = a 1 H 1 = L −1 1 where L 1 is the comoving horizon at t 1 [41]. Then, in terms of the dimensionless time t m = mt, (4.26) takes the form where we have also defined K ≡ kL 1 . The evolution of the mode functions after the onset of oscillations can well be described by the WKB approximation: where c k and ϕ k are numerical factors which we have determined by numerically solving the equation of motion (4.28) until t WKB m 1. We have verified that by choosing t WKB m = 50, the WKB approximation reproduces the numerical solution very well.
The power spectrum is given by [30] where and P θ (k) is the field power spectrum at some initial temperature T i defined such that where the cutoff is given by k cr = a i H i , and the normalization is fixed such that θ(x) 2 = π 2 /3. Note that at late times when all the relevant modes become non-relativistic we have ω k ≈ 1, and the power spectrum becomes independent of time 13 . 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 [38,41]. Again following [30] 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 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
We now present our results and make a comparison with the lattice result obtained in [41] 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 P(kL 1 1) ≈ C(kL 1 ) 3 with C 0.03 -0.05, which is the low-momentum behavior reported in [41]. 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 ≈ 3.5.
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 [41] h 2 Ω θ 0.019 g ρ (T 1 ) 70 (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 Θ i ≈ 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 eV, and solve (4.34) to get the correct relic density and found  [30] (blue lines). We compare this spectrum with the lattice result of [41] (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 lowmomenta behavior reported in [41] such that P(kL 1 1) ≈ C(kL 1 ) 3 with C 0.03 -0.05. The lattice result is obtained at t m = 2 × 10 3 in our units, however we have checked that the linearized result does not change significantly after t m ∼ 10 3 .    [30], and the lattice calculation of [41] with the spectra obtained in the Kinetic and Large Misalignment Mechanism for the same benchmark point.

S t a n d a r d M i s a l i g n m e n t K i n e t i c M i s a l i g n m e n t F r a g m e n t a t i o n a f t e r t r a p p i n g F r a g m e n t a t i o n b e f o r e t r a p p i n g N o n -p e r t u r b a t i v e
f ≈ 3.53 × 10 14 GeV. We then calculated the power spectra predicted by the Kinetic and the Standard/Large Misalignment Mechanisms and make the comparison with the postinflationary 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.

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.

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 [28] where ρ m is the matter density 14 including the over-densities, and ρ r is the radiation density which we assume to be homogeneous everywhere. The radiation pressure is p r = ρ 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.
where ρ m is the background matter density, and δ ≡ ∆ρ m /ρ 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 Let us parametrize r(t) by 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 ≡ a/a eq . After some algebra, one can obtain a differential equation for the evolution of b(y): [28] y(1 + y) where δ i is the initial over-density in the region δ i = δρ m (y i )/ρ m (y i ) set at some early time when b(y i ) = 1. By deriving this result, we have used where ρ eq = 2ρ m,eq = 2ρ r,eq is the total average energy density at the matter-radiation equality.
For a given value of the initial over-density δ i , one can solve (5.5) numerically with the initial conditions b(y i ) = 1, and b (y i ) = 0 where y i 1. Initially, the radius will expand with the Hubble flow. Then it will start to deviate, and at y ta it turns around, i.e. the expansion will stop, and the radius will start shrinking. This happens whenṙ = 0, or b (y ta ) + b(y ta )/y ta = 0. Numerical solution shows that the turnaround happens at [30] 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 δ becomes infinite. Again, by numerical solution we observed that Numerical calculation of δ i y ta , and δ i y c for a bunch of initial overdensities is shown on the left plot of Figure 13. 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 ? [79] 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 In the linear approximation we define δ(y) = δ i + ∆δ(y), and evaluate all expressions at linear order in ∆δ(y). Then (5.9) implies b(y) = 1 − ∆δ(y)/3(1 + δ i ), and (5.5) becomes where 2 F 1 is the hypergeometric function. So in the linear theory, δ(y) is given by Since the non-linear theory predicts that the initial fluctuation δ i will collapse at y c , we can calculate the prediction of the linear theory for the overdensity at collapse by replacing δ 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 The quantity δ c is called the critical density at collapse. Small fluctuations δ i 1 collapse during the deep matter era where y c 1. In this case we recover the well-known result: On the other hand, for large fluctuations we can consider the y c 1 limit to approximate the hypergeometric function at 1 + y c . Then we find Suprisingly, the approximate expression for small and large fluctuations are the same. So for all regimes, one can use δ c ≈ (1.1/x c )(1+3x c /2) as the expression for the critical density at collapse [79]. The plot of the critical density δ c as a function of y c is plotted on the right plot of Figure 13.

Press-Schechter Formalism
The formation of the dark matter halos can be studied analytically via the Press-Schechter (PS) formalism [80]. It relies on the following postulate: The fraction of matter F which is inside collapsed structures of comoving size larger than R at any given scale factor a is the same as the probability P of finding an overdensity δ R (x, a) > δ c (a) where δ R (x, a) is the overdensity smoothed at the scale R, and δ 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 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 where δ is the energy density contrast. By a straightforward calculation we can obtain the variance of the smoothed density contrast as where W is the Fourier transform of the window function and P δ 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 Its Fourier transform is (5.20) With this window function, the smoothed density contrast δ R (x) is just the integral of the density contrast over a spherical region of comoving size R centered at x. Then the mass enclosed in the region is unambiguously given by where the approximation holds for small fluctuations. We will use this relation even if the fluctuations are large 15 . 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 as 16 where k min and k 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 δ GI CDM given by [55] δ GI CDM (k; a) = 9Φ k (0) so that the dimensionless power spectrum becomes 17 During the matter era, we use where Ω m ≈ 0.32, K ≡ √ 2k/k eq with k eq is the wavenumber which enters the horizon at matter-radiation equality, and T BBKS is the BBKS transfer function [81]: 15 The 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 [30]. 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. 16 In the case of complete fragmentation we took P δ = P quad δ in the second term of (5.22) since for all the simulated modes the linear contribution is negligible. 17 The 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 2 k . We are using the gauge-invariant one in this calculation in order to have a well-defined integral for the variance. 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 eV and f = 1.3 × 10 14 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.
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 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 nonlinear 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.

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 P(δ R (a) > δ c (a)). The easiest approach is to assume that δ(x, a), hence δ R (x, a) is a Gaussian random field. Then the probability density function (PDF) for δ R is given by . In the case of CDM, the variance σ 2 R goes to infinity as R → 0. This implies that all the matter in the universe must finally be in virialized objects so we expect F (> M ; a) → 1 as M → 0 [82]. However, the PS result (5.28) has the limit F (> M ; a) → 1/2 as M → 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 [83]. This effect can be included by using the excursion set formalism [84], for a nice review see [85]. The correct expression is 18 The advantage of using this form is that one can replace f PS (ν) 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 [86] 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 [87]. In this work we will use the PS mass function (5.32) since it is more general and simpler. 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 eV and f = 1.3 × 10 14 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 We also present our results of the halo mass function in terms of the dimensionless variable [30]  in Figure 17. This approximately gives the fraction of dark matter that reside in halos with masses between [ln M, ln M + d 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  Table 1: The list of benchmarks used in constructing the halo mass functions in Figure 18 and the halo spectra shown in Figure 20.
to the CDM prediction for which X(M ) grows with M 19 . 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.
In Figure 18, we show a comparison of the dimensionless halo mass functions cor- 19 The CDM prediction for X(M ) has a global maximum around M ∼ 10 14 M , see Figure 18.  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 [51] and [88], 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 [87,89,90]. 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.

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 [91] ρ(r; z) = δ char ρ c (z) (r/r s )(1 + r/r s ) 2 , The quantities δ 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 [91], which can also be used for axion miniclusters [79]. Consider a halo identified, or measured, at redshift z = z 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 where ρ m,0 is the CDM/axion density today 22 , and a 0 is the scale factor today. Now assume that δ char (z ind )ρ c (z ind ), and hence ρ s , is proportional to the energy density of the universe at the time of collapse: δ char (z ind )ρ c (z ind ) = 4ρ s (z ind ) = C(z col )ρ c (z col ), (5.45) where C(z col ) is in general z col -dependent proportionality factor. Both C and f are free parameters which need to be fit with the N-body simulations. For both CDM [91] and axion miniclusters [79], it is found that f = 10 −2 provides a good fit to the simulations.
Hence, we will also use this value in our calculations. To determine C(z col ) we have used the CDM simulations compiled in [95], and axion minicluster simulations in [40]. 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 ρ s , can be calculated. 20 The 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 [92]. 21 The parameter z col has nothing to do with the redshift of the linear collapse yc = (1 + zeq)/(1 + zc) that we have introduced in Section 5.1. 22 We are assuming that CDM/axion make all of matter.

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 [1] and [42]. 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 [96][97][98]. Instead, the density profiles are described by solitons 23 which are the eigenstates of the time-independent Schrödinger-Poisson equations [101,102]. The soliton profile can be fitted by [97] ρ sol (r) = ρ sol (0) 1 + 2 1/8 − 1 r r 1/2 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 , ρ s and M s . Then one can show that the scale density and the scale mass are related to each other by This ground state configuration is independent of the power spectrum, and is reached rapidly by a process called gravitational cooling [103,104].

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 sol s -ρ sol s relation (5.48) as an upper bound on the scale density for a given scale mass.
We show the results of the halo spectrum in Figure 19. We assumed a constant ALP mass, fixed it to m = 10 −15 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  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 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 (nonlinear 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 µ 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 α 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.
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 eV (red, left), m = 10 −10 eV (blue, middle), and m = 10 −15 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) 24 . 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 -ρ 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.

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.

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 -ρ s plane presented in [51].
• 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 [105]. The weak lensing of a compact halo can induce localized distortions in the proper motion µ of the background light sources that can be observed by GAIA [106] and SKA [107]. We call this µ 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 α that can be probable Theia [108]. We call this α lensing and show the observable region via thin pink line in Figures 19 and 20.
• Photometric microlensing: Another way of probing the halo spectrum is the photometric microlensing where the brightness of a background light source changes due to the passage of a compact halo [109,110]. The potentially observable region as obtained in [51] is shown by the thin orange line in Figure 20.
• Pulsar Timing Arrays: The Pulsar Timing Arrays (PTAs) are powerful probes of DM substructure at small scales [73,111,112]. A transit of a compact object near the timing system can modify the frequency of the pulsar in two ways. First, it can modify the photon geodesic along the line of sight causing a Shapiro time delay [113]. Second, the compact objects can lead to an acceleration of the pulsar or the Earth modifying the period of the pulsar via the Doppler effect [114]. We show the potentially observable region obtained in [112] via the thin green line in Figure 20.
• Diffraction of Gravitational Waves: Final observable that we comment on is the diffraction of gravitational waves. A distribution of compact DM halos can also induce distortions in the amplitude and the phase of the gravitational waves emitted from a merger of two black holes [115]. The potentially observable region via aLIGO as obtained in [51] is shown via dark yellow lines in Figures 19 and 20.
All the observable prospects mentioned above assume a monochromatic halo mass distribution, in other words X(M ) is taken to be proportional to a δ-function. This means that, for a given M s , the lines are derived under the assumption that some fraction of dark matter resides in halos with scale mass M s and scale density ρ s . In [51] this fraction is chosen to be 0.3 which we also use. The monochromatic assumption is necessary to draw the observable regions on the M s -ρ s plane. The observable prospects will be modified with more realistic mass distributions like we have derived in Section 5.3. For the halometry observations, such as µ-and α-lensing it is expected that this will not modify the qualitative features too much [105]. However, taking into account realistic sub-halo mass distributions reduces the prospects of the PTAs significantly, at least for CDM [73,112].

Consequences for terrestrial experiments
If the ALPs have couplings to the Standard Model, they can be detected via direct detection terrestrial experiments, most notably via Haloscopes 25 . These experiments search for a monochromatic signal with frequency set by the ALP mass m 0 , and amplitude set by the local DM density ρ DM ≈ 1.1 × 10 −2 M pc −3 . The coherence time of the DM signal is given by τ ∼ 2π m 0 v 2 where v ∼ 10 −3 is the galactic virial velocity of dark matter [116]. If a significant fraction of DM in the Milky Way is bound to minihalos that are much denser that the local DM density, they can have a big impact on the terrestrial DM searches if they pass through the Earth [51,92]. During the crossing of a minihalo, the DM signal is boosted by a factor of B ∼ ρ s /ρ DM , and is coherent for a much longer time. However, detecting such a signal would also require some modifications in the experimental setups. See [51] for more details on this point.
To assess the relevance of such an effect, we need to quantify how likely a microhalo will pass through the Earth. This information is encoded in the differential encounter rate given by dΓ enc d ln M = d n d ln M σ enc v vir ≡ ρ DM M X(M )σ enc v vir , (6.1) 25 For a detailed list of haloscope experiments and the corresponding references we refer the reader to our companion paper [14].
where σ enc ∼ πr 2 s is the encounter cross-section, v vir ∼ 200 km/s is the virial velocity of DM inside the Milky-Way, and d n/d ln M is the sub-halo mass function (sHMF). The latter tells us the distribution of sub-halos of the Milky Way and it is different than the HMF that we have calculated in Section 5.3 that tells us about the global distribution of halos in the whole universe. Therefore, a precise calculation of the encounter rate requires the calculation of the sHMF which we have left for future work. Here we assume a monochromatic sHMF by taking X(M ) = Ω sub δ(ln M − ln M s ), (6.2) where Ω sub is the fraction of DM in the Milky Way that is bound to minihalos. With this assumption we can estimate the encounter rate as where we have assumed that the minihalos have the NFW profile and used (5.37). We see that this effect is relevant mostly for very light sub-halos, and therefore for heavy ALP masses m 0 10 −5 eV. The crossing time of a minihalo can be estimated as We show our results in Figure 21. The solid and dashed black contours show the encounter rate and crossing time respectively. We took Ω sub = 0.3 as the fraction of DM bound to sub-halos, consistent with the previous section. We also show the halo spectra for three benchmark points given in the legend. The solid and dot-dashed lines show the spectra in the case of Kinetic and Large Misalignment respectively with the same m 0 and f values. The dashed lines show the soliton profiles. The halos that have densities ρ s 10 3 ρ DM are expected to be tidally distrupted before the encounter according to the estimate in [51].
We see that the probability of a minihalo encounter during the lifetime of an experiment is rather low, even with highly optimistic assumptions. Repeating the calculation with a more realistic sHMF will likely give a lower encounter rate according to the analysis done in [92] for a different model. We plan to revisit this in a future work.

Conclusion
In summary, we have analysed the implications of ALP fragmentation in the Kinetic Misalignment Mechanism (KMM) for the formation of ALP mini-clusters, and the observational consequences. KMM is particularly appealing as it makes the low f /large coupling region, which is accessible to terrestrial experiments, viable for dark matter. This work is a follow-up to our companion article [14]. The main logic and results of this paper are summarised in Figure 1, quoting the key equations. We solved the equations for the mode functions of the ALP fluctuations (equations (3.10) and (3.11)), obtained the initial conditions (3.17), and solved for the late time evolution (3.27). A key result is the comparison of the density contrasts in KMM versus the Standard Misalignement Mechanism (SMM) and Large Misalignment Mechanism (LMM) before fragmentation starts (equations (3.47) and (3.48)). This shows the large suppression in SMM and LMM in comparison with KMM, initial fluctuations being much larger in KMM.
An important figure is Figure 4, showing the full evolution of the modes, comparing CDM with KMM. The period just before fragmentation is illustrated in Figure 8, showing clearly that the large initial fluctuations in KMM is what makes the difference.
Our two main final predictions are the halo mass functions ( Figure 18) and the halo spectra ( Figure 20). This leads to promising prospects for experiments such as GAIA, SKA, Theia, PTA, aLIGO. On the other hand, we do not expect that the formation of dense ALP halos will impact direct detection experiments such as haloscopes as we found that the probability of a minihalo encounter during the lifetime of an experiment is rather low.
We will provide a more precise study of the predictions for the QCD axion, relying on lattice simulations, in a future work. A key question for the relevance of this work is the motivation for the kinetic misalignment mechanism. For this, the general class of UV completions realising this framework should be scrutinised. This is what is done in [65] where we analyse in detail the parameter space in the (m, f ) plane.
In conclusion, there are great experimental prospects for probing the early evolution of the ALP field and thus access high energy scales through a variety of detection approaches, combining haloscopes, helioscopes, light-shining-through-the-wall experiments as well as astronomical observations as the ones we have discussed in this paper. Here the scale factors a MD and a in are irrelevant as long as they are chosen at deep matter and radiation domination respectively. The normalization is chosen such that T (0) = 1. The main idea for introducing a transfer function is to separate the evolution of the modes from their initial conditions. The evolution at later times a > a MD can be described by where D(z) is called the growth function. It is defined such that D(z) = (1 + z) −1 in a matter-dominated universe without a cosmological constant Ω Λ = 0. A non-zero Ω Λ suppresses the growth at very late times so that todays value of the growth factor is D(0) ≈ 0.79 instead of 1. Since this suppression happens very late, and most of the structures that we are interested in collapse deep in the matter era, we neglect the suppression from the vacuum energy and assume D(z) = (1 + z) −1 .
To get the density constrast, we make use of the Poisson equation given by which is valid for all modes not just the sub-horizon ones. From this we can solve for the Fourier modes of the density constrast as In this work we neglect the spectral tilt by taking n s = 1, and we also neglect the baryons so the matter power spectrum is the same as the CDM power spectrum. Finally by changing to the dimensionless power spectrum P = k 3 P (k)/2π 2 , and defining K ≡ √ 2k/k eq we get (5.25).
C Fitting procedure to obtain C(z col ) To determine C(z col ) we have used two sets of simulations: • A set of CDM simulations compiled by Conde and Prada [95].  Prada [95] for CDM, and Eggemeier et al. [40] for axion miniclusters.
• An axion minicluster simulation by Eggemeier et al. [40] based on the lattice calculation of the post-inflationary scenario for the QCD axion that is performed in [38].
Both sets of simulations give the results for the concentration parameter c 200 as a function of the halo mass M 200 . To perform the matching, first we did calculate the collapse redshifts for the tabulated halo masses. For CDM halos we use the power spectrum (5.25), while for the axion miniclusters we took the latetime power spectrum obtained in [38], calculated the variance using the spherical top-hat window function, and finally evaluated until today as [73,79]  where σ i is the initial variance obtained from the power spectrum. The reason that the Jeans scale plays no role in this evolution is because all relevant modes are well below the Jeans scale at the matter-radiation equality.
After the collapse redshifts are calculated, it becomes a trivial task to determine C(z col ) using the given concentration parameter data. For CDM halos that collapse around z ∼ 10 we found that C ∼ O(1) × 10 3 , and it is an increasing function of z col . On the other hand, axion miniclusters collapse around z col ∼ 10 3 for which C ≈ 9.76 × 10 4 did provide the best fit. Finally by combining the CDM and the axion minicluster data we obtained the following fit function for C(z col ): C(z col ) =