aainstitutetext: Deutsches Elektronen-Synchrotron DESY, Notkestr. 85, 22607
Hamburg, Germanybbinstitutetext: Department of Physics, Osaka University,
Toyonaka, Osaka 560-0043, Japanccinstitutetext: II. Institute of Theoretical
Physics, Universität Hamburg D-22761, Germany

# ALP Dark Matter from Kinetic Fragmentation: Opening up the Parameter Window

Cem Eröncel b Ryosuke Sato a,c Géraldine Servant a,c Philip Sørensen
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The main mechanism responsible for Axion-Like-Particle (ALP) production in the
early universe is the so-called misalignment mechanism. Three regimes have
been investigated in this context: standard misalignment, large misalignment
and kinetic misalignment. The latter applies if the axion inherits a large
initial velocity in the early universe, such that the field rolls through many
wiggles during its evolution, before it gets trapped in one minimum. This
largely opens the region of parameter space for ALP dark matter towards higher
values for the axion-photon coupling, which can be probed by the whole set of
next decade’s upcoming experiments. In fact, almost the entire parameter space
in the [mass, decay constant] plane can now accommodate dark matter. In this
paper, we show that in kinetic misalignment, the axion field is almost always
entirely fragmented, meaning that the energy density of the homogeneous field
is redistributed over higher-mode axions. We present a general model-
independent analytical description of kinetic fragmentation, including
discussion of the modified initial conditions for the mode functions due to
the axion’s initial velocity, and how they impact the growth of the adiabatic
fluctuations. We calculate precisely the parameter regions corresponding
respectively to standard misalignment, kinetic misalignment with weak
fragmentation, fragmentation after trapping and fragmentation before trapping.
While axion fragmentation can impact the precise determination of the relic
abundance, another main observational implication is the formation of much
denser compact axion halos, that is described in a companion paper. We also
point out a new gravitational-wave signature that arises in the large
misalignment regime with complete fragmentation and could be seen in
measurements of $\mu$ distortions in the Cosmic Microwave Background.

DESY 22-106
OU-HET-1148

##  1 Introduction

Axion-Like-Particles (ALPs) appear in many extensions of the Standard Model.
Their physics case has been extensively scrutinised in the last decades,
especially in the last few years which have seen an increased interest in
novel detection strategies and new experimental proposals Irastorza:2018dyq ;
Adams:2022pbo . ALPs are pseudo-Nambu-Goldstone bosons, arising from the
spontaneous breaking of a global $U(1)$ symmetry. As such, they are naturally
light particles, with a mass $m$ much smaller than the abelian symmetry
breaking scale $f$, also referred to as the ALP decay constant. ALPs are
characterised by their small mass $m\ll f$, and by their small coupling to
Standard Model particles which scales as $1/f$. The axion mass is generated at
some lower scale $\Lambda_{b}\ll f$ typically related to some new strong
dynamics that break the $U(1)$ symmetry explicitly, leading to
$m\propto\Lambda_{b}^{2}/f$. For a generic ALP, $m$ and $f$ are independent
parameters. In the case of the QCD axion, they are related as
$mf\propto\Lambda_{QCD}^{2}$. The best-known example of an axion is in fact
the QCD axion, which arises as a consequence of the Peccei-Quinn solution to
the strong-CP problem Peccei:1977hh ; Weinberg:1977ma ; Wilczek:1977pj . ALPs
have a characteristic two-photon vertex. This coupling is at the basis of most
search strategies which rely on axion-photon conversion in external magnetic
fields.

Particularly well-motivated is the case that ALPs compose the dark matter of
the Universe Preskill:1982cy ; Abbott:1982af ; Dine:1982ah . Because of their
small mass and small couplings, ALPs are cosmologically stable. In most of the
literature, ALPs are produced as extremely cold relics in the early universe
through the misalignment mechanism. The energy density of the axion field is
stored into the oscillations of the zero-mode homogeneous axion field
$f\theta(t)$. The scalar field equation takes the form
$\ddot{\theta}+3H(t)\dot{\theta}+m^{2}(t)\sin\theta=0$. The initial conditions
are usually taken to be $\theta(0)=\theta_{i}$ where $\theta_{i}\in[-\pi,\pi]$
is referred to as the initial misalignment angle, and $\dot{\theta}(0)=0$. At
early times, the field is frozen due to Hubble friction. Oscillations start as
soon as the expansion rate of the universe $H(t)$ drops below the axion mass.
From that time, the field behaves effectively as pressureless cold matter. The
axion contribution to the dark matter density today crucially depends on the
initial displacement of the axion field with respect to the minimum of the
axion potential. We show the corresponding predictions for different
assumptions of this initial angle in Fig. 1. Constraints can be represented in
the $[m,1/f]$ plane. The testable region of parameter space of both the QCD
axion and of generic ALPs is still restricted. However, there are many plans
for testing further the parameter space. Another ALP production mechanism is
from the decay of cosmic strings or domain walls into axion particles, which
roughly predicts comparable regions of parameter space as those from the
misalignment mechanism. Experiments such as ADMX can probe some parts of the
parameter space of the QCD axion, although not in the region where it is Dark
Matter, while other regions of generic ALP DM are out of reach of most future
experiments, see Fig. 1.

In this paper, we consider a new mechanism for ALP production in the early
universe, so far overlooked, in which ALPs are first produced relativistically
while the energy density initially stored in the homogeneous zero mode of the
axion field gets transferred into relativistic axion particles. This happens
generically if the axion field has initially sufficiently large velocity to
roll over many wiggles of the periodic axion potential until the expansion
rate of the universe is sufficiently small to enable axion fragmentation. This
can happen at relatively low temperatures and leads to a viable scenario of
Cold Dark Matter if the fragmentation temperature is below a TeV. These axion
particles quickly cool down and dominate the energy density. The corresponding
axion periodic potential emerges from some new confining dynamics happening in
a dark sector at energies between a keV and a TeV.

Recently, the usual assumption that the axion field is initially static in the
early universe has been questioned Co:2019jts ; Chang:2019tvx . In particular,
it was shown that such initial large velocity can delay the onset of the usual
misalignment mechanism for axion DM production. The so-called axion kinetic
misalignment has been put forward to predict axion dark matter with a larger
coupling than in the original setup, which is excellent news for experiments
which now have stronger motivations to search for ALP DM. We show that there
is an additional effect which takes place in the axion kinetic misalignment
framework, called axion fragmentation Fonseca:2019ypl . Axion fragmentation
was studied in the context of the relaxion in Fonseca:2019lmc and for dark
matter in axion monodromy models in Jaeckel:2016qjp ; Berges:2019dgr .

While most axion cosmology literature usually assumes that the axion field
does not have any initial velocity, such initially large velocity is actually
a natural outcome in a class of UV completions extensively discussed in
Co:2019jts ; Co:2019wyp ; Co:2020jtv ; Co:2020xlh ; Gouttenoire:2021jhk .
Considering the early evolution of the full complex Peccei-Quinn scalar field,
if the radial mode has a large initial VEV, it can induce a kick for the
angular (axion) mode that starts rotating for a long time (similar to what
happens in Affleck-Dine baryogenesis Affleck:1984fy ; Dine:1995kz ). In this
paper, we remain agnostic about the UV completion. We work in full generality,
only considering the late-stage evolution when the complex scalar field has
reached the bottom of the potential. The radial mode of the complex scalar
field can be ignored, we only consider the axion degree of freedom after it
has acquired a mass. Because of its large initial velocity, it can overshoot
the potential barrier and keep rotating until it eventually fragments
entirely.

It is the main goal of this paper to investigate in detail and in full
generality (model-independently) the properties of ALP DM produced by
fragmentation in the kinetic misalignment mechanism. We generalise the
analysis of Ref.Fonseca:2019ypl in two ways. First, we focus on the axion
being DM (while Fonseca:2019ypl did not impose that the axion is DM). Second,
we not only investigate axion fragmentation before the axion is trapped by the
barriers but also show that fragmentation can take place after trapping. We
provide a discussion of the initial conditions for the mode functions and
derive the field power spectrum. The distinctive feature in this scenario is
the momentum distribution of the axion particles. Ultimately, our analysis
sets the basis to derive the phenomenological implications of this modified
power spectrum, which will be presented in Eroncel:2022c . We will provide
predictions in the explicit UV completions in another companion publication
Eroncel:2022a , deriving the precise model parameter space regions for ALP
dark matter through kinetic fragmentation.

The plan of this article is as follows. In Section 2 we provide an analytical
derivation of the ALP relic abundance and trapping temperature in the kinetic
misalignment framework. The latter is crucial when calculating the power
spectrum. The analytical theory of parametric resonance in kinetic
misalignment is derived in section 3. The determination of the various regions
in the [axion mass, decay constant] parameter space are determined in Section
4 where the impact of fragmentation on the relic abundance prediction is also
discussed. All phenomenological constraints are compiled in Section 5. Section
6 presents the prediction of the gravitational-wave spectrum induced by
fragmentation. We summarise our results by listing the important equations and
figures in Section 7. Technical details are gathered in appendices A, B and C.
Appendix D lists all the experimental bounds and projections that are reported
in our plots together with the corresponding references.

Figure 1:  The experimental landscape in the hunt for ALPs assuming a KSVZ-
like axion-photon coupling $g_{\theta\gamma\gamma}^{\rm KSVZ}$ given in
(5.10). Coloured regions are excluded. The thin lines indicate the
sensitivities of future experiments. Used data is listed in Appendix D. Orange
constraints apply to any ALP while the green ones assume the ALP is DM. The
yellow thick line corresponds to the QCD axion. The four other parallel
straight lines indicate the correct dark matter relic abundance contours for
different assumptions of the initial misalignment angle. Above the thick
orange line, the axion produced from the standard misalignment mechanism is
under-abundant to explain DM.

##  2 ALP dark matter from kinetic misalignment

We consider the cosmological evolution of an ALP field $\theta$ whose
Lagrangian is given by

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

where $f$ is the vacuum expectation value of the complex scalar field radial
component. The metric is taken to be the flat Friedmann-Lemaitre-Robertson-
Walker metric111In general, the metric does also have curvature perturbation
terms. These will determine the initial conditions for the mode functions as
we will demonstrade in Section 3.3.

$\differential{s}^{2}=-\differential{t}^{2}+a^{2}(t)\delta_{ij}\differential{x}^{i}\differential{x}^{j}.$
(2.2)

We decompose the ALP field $\theta(t,\mathbf{x})$ into a homogeneous mode
$\Theta(t)$, and small fluctuations $\delta\theta(t,\mathbf{x})$, where the
latter can be expanded into Fourier modes as

$\delta\theta(t,\mathbf{x})=\int\frac{\differential[3]{k}}{(2\pi)^{3}}\theta_{\mathbf{k}}(t)e^{-i\mathbf{k}\cdot\mathbf{x}}.$
(2.3)

In this work we assume that the Fourier modes $\theta_{\mathbf{k}}(t)$ evolve
independently. We also separate their time evolution from their statistical
properties by defining

$\theta_{\mathbf{k}}(t)=\hat{\theta}_{\mathbf{k}}\theta_{k}(t),$ (2.4)

where $\hat{\theta}_{\mathbf{k}}$’s carry the statistical properties, and
$\theta_{k}$’s are $c$-number functions which depend only on the magnitude of
the momentum $k\equiv\absolutevalue{\mathbf{k}}$ due to the isotropy of the
equations of motions in the linear approximation. We also assume that
$\hat{\theta}_{\mathbf{k}}$’s obey statistical homogeneity and isotropy:

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

The initial conditions for the mode function $\theta_{k}$ will be determined
by the initial field power spectrum $P_{\theta}$ which depends on the physical
process that creates the initial conditions for the Fourier modes. We study
the initial conditions in Section 3.3.

If we neglect the back-reaction, the equations of motion of the homogeneous
mode $\Theta$ and of the fluctuations $\theta_{k}$ are decoupled from each
other. These are respectively given by

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

and

$\ddot{\theta}_{k}+3H\dot{\theta}_{k}+\quantity[\frac{k^{2}}{a^{2}}+m^{2}(T)\cos\Theta]\theta_{k}=0.$
(2.7)

We start by studying the evolution of the homogeneous mode, and derive today’s
value of the relic density in Sections 2.1 and 2.2. Then we discuss the
implications of the parametric resonance in Section 3.1. We conclude this
section by providing a semi-analytical estimate of the back-reaction of the
fluctuations on the homogeneous mode. We comment on the implications in
Section 3.2.

###  2.1 Derivation of the relic density today from the adiabatic invariant

In the Kinetic Misalignment mechanism Co:2019jts , we work under the
assumption that the ALP field has a large initial velocity such that its
initial kinetic energy is much larger than its potential energy. Therefore it
travels over many barriers without being trapped. As the universe expands, the
kinetic energy will redshift while the potential energy will grow since it is
proportional to the axion mass. At some temperature $T_{\ast}$, the energy in
the homogeneous mode $\rho_{\Theta}$ becomes smaller than the height of the
barriers. After this time the ALP field cannot continue rolling, as it is
_trapped_ by the barrier. The trapping temperature $T_{\ast}$ is defined by
the condition222From now on the $\ast$-subscript denotes the quantities
evaluated at $T=T_{\ast}$.

$\rho_{\Theta}(T_{\ast})=\frac{1}{2}f^{2}\dot{\Theta}^{2}(T_{\ast})+m^{2}(T_{\ast})f^{2}\quantity[1-\cos(\Theta(T_{\ast}))]=2m^{2}(T_{\ast})f^{2}\quantity[1-\cos(\Theta(T_{\ast}))].$
(2.8)

We introduce the dimensionless parameter $\varrho(t)$ as the ratio between the
total ALP energy density and its potential energy333We will use the time and
temperature interchangebly depending on the nature of the quantity that we
want to calculate. We always assume radiation domination so $t=(2H)^{-1}$, and
$H(T)\propto g_{\ast}(T)^{1/2}T^{2}$ where $g_{\ast}$ is the effective number
of relativistic degrees of freedom. We assume that $g_{\ast}$ is constant
throughout the parametric resonance so we can change from temperature to time
and vice versa by $T^{2}\propto t^{-1}$.

$\varrho(t)\equiv\frac{\rho_{\Theta}}{2m^{2}(t)f^{2}}=\frac{1}{4}\frac{\dot{\Theta}^{2}}{m^{2}(t)}+\sin^{2}\quantity(\frac{\Theta}{2}).$
(2.9)

The field is trapped at $\varrho(t_{\ast})=1$. If $\varrho>1$ it keeps
rolling, while for $\varrho<1$ it is oscillating around its minimum.

The evolution of $\varrho(t)$ can be studied analytically by using _action-
angle formalism_ goldsteinClassicalMechanics2001 ; landauMechanicsVolume1976 ,
see also Masso:2005zg . Neglecting the back-reaction, the evolution of the
homogeneous mode is governed by the Lagrangian density

$\mathcal{L}_{0}=f^{2}a^{3}(t)\quantity[\frac{1}{2}\dot{\Theta}^{2}-m^{2}(t)\quantity(1-\cos\Theta)].$
(2.10)

This Lagrangian has a discrete shift symmetry given by
$\Theta\rightarrow\Theta+2\pi$, so the motion in $\Theta$-space can be
considered periodic. In the following, we define the period to be the time it
takes for the field to roll from one maximum to the other.

Now we assume that the parameters carrying explicit time dependence, namely
$a$ and $m$, change _adiabatically_ , in other words they remain approximately
constant during a single oscillation. This assumption is justified if the
Hubble time scale is larger compared to the period of the motion. We will
justify this assumption later. In the case of periodic motion where the
Lagrangian changes adiabatically, the _action variable_ $J$ defined by

$J\equiv\oint\Pi_{q}\differential{q},$ (2.11)

is a constant of motion, where the integral is over a single period, and

$q=f\Theta\mbox{\quad
and\quad}\Pi_{q}=\partialderivative{\mathcal{L}_{0}}{\dot{q}}=a^{3}(t)f\dot{\Theta}.$
(2.12)

Assuming that the scale factor and the axion mass are constant during one
period, $J$ becomes

$J\approx
f^{2}a^{3}(t)\oint\differential{\Theta}\dot{\Theta}=2m(t)f^{2}a^{3}\oint\sqrt{\varrho-\sin^{2}\quantity(\frac{\Theta}{2})}.$
(2.13)

The integral can be performed analytically and one gets

$J\approx
m(t)f^{2}a^{3}(t)\times\begin{cases}8\sqrt{\varrho}\,\textrm{E}\quantity(1/\sqrt{\varrho}),&\varrho>1\\\
8\quantity[(\varrho-1)\textrm{K}(\sqrt{\varrho})+\textrm{E}(\sqrt{\varrho})],&\varrho<1\end{cases},$
(2.14)

where K and E are complete elliptic integrals of the first and second kind
respectively444We use the following definitions for the elliptic integrals:
$\textrm{K}(k)=\int_{0}^{\pi/2}\frac{\differential{\varphi}}{\sqrt{1-k^{2}\sin^{2}\varphi}}\mbox{\quad
and\quad}\textrm{E}(k)=\int_{0}^{\pi/2}\differential{\varphi}\sqrt{1-k^{2}\sin^{2}\varphi}$
Note that most software packages such as Mathematica and scipy 2020SciPy-NMeth
uses $m=k^{2}$ instead of $k$ as their argument when defining elliptic
integrals. . By using the limiting behavior of these as $\varrho\rightarrow 1$
we can show that

$\lim_{\varrho\rightarrow 1^{+}}J(\varrho)=\lim_{\varrho\rightarrow
1^{-}}J(\varrho)=8m(t)f^{2}a^{3}(t).$ (2.15)

The adiabatic invariant $J$ remains approximately constant throughout the
whole evolution. This makes $J(\varrho)$ a perfect quantity to compare the
early and late time behaviors. At early times when $\varrho\gg 1$ we have

$J(\varrho\gg 1)\approx 2\pi f^{2}\dot{\Theta}a^{3},$ (2.16)

whereas at late times we get

$J(\varrho\ll 1)\approx\pi\frac{\rho_{\Theta}a^{3}}{m_{0}},$ (2.17)

where $m_{0}$ is the zero-temperature axion mass. Let us introduce the _yield_
quantity defined by

$Y\equiv n_{\Theta}/s,$ (2.18)

where $s$ is the entropy density of the universe, and

$n_{\Theta}=f^{2}\dot{\Theta}$ (2.19)

is the Noether charge of the axion shift symmetry. This quantity is conserved
during the early evolution when the ALP mass is negligible so the equation of
motion of the ALP field is just $\ddot{\Theta}+3H\dot{\Theta}\approx 0$. The
conservation of $J(\varrho)$ implies that (2.16) and (2.17) should be equal
which gives

$\frac{\rho_{\Theta,0}}{s_{0}}\approx 2m_{0}Y.$ (2.20)

where $\rho_{\Theta,0}$ and $s_{0}$ are respectively the relic density and
entropy density today. This equation is also derived in Co:2019jts , however
the numerical factor there was obtained numerically. Here we did present an
analytical derivation.

From this result we can obtain a simple formula for the fractional energy
density of the ALP field today as

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

where we took $h^{2}\Omega_{\rm DM}=0.12$ Aghanim:2018eyx . Note that the
yield does only depend on the zero temperature axion mass, not on the high-
temperature behavior. The latter will be important in determining the trapping
temperature $T_{\ast}$. For the QCD axion the zero-temperature mass and the
axion decay constant are related by $m_{0}f\approx(75.6$\mathrm{MeV}$)^{2}$
Borsanyi:2016ksw . Then the above relation implies

$h^{2}\Omega^{\rm QCDaxion}_{\phi,0}\approx h^{2}\Omega_{\rm
DM}\quantity(\frac{10^{9}\,$\mathrm{GeV}$}{f})\quantity(\frac{Y}{40}).$ (2.22)

Before closing this subsection, we mention in which conditions the adiabacity
assumption is justified. To do this we need to calculate the periods of the
motion before and after trapping. By neglecting the energy loss due to the
Hubble expansion during a single oscillation, the periods $\mathcal{T}_{>}$
before and $\mathcal{T}_{<}$ after trapping can be derived as

$\mathcal{T}_{>}(\varrho)=\frac{2}{m\sqrt{\varrho}}\textrm{K}(1/\sqrt{\varrho})\mbox{\quad
and\quad}\mathcal{T}_{<}(\varrho)=\frac{2}{m}\textrm{K}(\sqrt{\varrho}).$
(2.23)

At early times when $\varrho\gg 1$ we have

$\mathcal{T}_{>}(\varrho\gg
1)\approx\frac{\pi}{m}\sqrt{\frac{1}{\varrho}}\approx\frac{2\pi}{\dot{\Theta}},$
(2.24)

so the adiabacity condition reads $\dot{\Theta}>2\pi H$. At early times
$\dot{\Theta}\propto a^{-3}$ and $H\propto a^{-2}$ so it is sufficient if this
condition is satisfied at trapping $\dot{\Theta}=2m_{\ast}$. This yields
$m_{\ast}>\pi H_{\ast}$ as the adiabacity condition. On the other hand, at
late times when $\varrho\ll 1$ we have

$\mathcal{T}_{<}(\varrho\ll 1)\approx\frac{\pi}{m},$ (2.25)

so the adiabacity condition becomes $m>\pi H$. Again, it is sufficient that
this condition is satisfied at trapping which yields

$\frac{m_{\ast}}{H_{\ast}}\gtrsim\pi,$ (2.26)

as the adiabacity condition. As we will show explicity in Section 4, the most
interesting region of the parameter space is where $m_{\ast}/H_{\ast}\gg 1$
which justifies the adiabacity assumption. We note that even for
$m_{\ast}/H_{\ast}\gg 1$, the adiabacity is broken for a period around
$\varrho=1$, but this period is shorter for larger $m_{\ast}/H_{\ast}$.

###  2.2 Trapping temperature

The trapping temperature is an important quantity since fragmentation happens
around the same time as trapping. We recall that the trapping happens at
$\varrho=1$ for which

$J(\varrho=1)=8m_{\ast}f^{2}a_{\ast}^{3}.$ (2.27)

By comparing this with (2.17), and using conservation of $J(\varrho)$ we
obtain

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

where $\Lambda_{b,0}^{4}=m_{0}^{2}f^{2}$ is the zero-temperature barrier
height. In order to calculate the trapping temperature we need to specify the
scaling of the axion mass with respect to the temperature. We consider the
following ansatz:

$m^{2}(T)=m_{0}^{2}\times\begin{cases}\quantity(T/T_{c})^{-\gamma},&T>T_{c}\\\
1,&T\leq T_{c}\end{cases},\qquad T_{c}\equiv 2.12\times\Lambda_{b,0}.$ (2.29)

The idea behind this choice is for $\Lambda_{b,0}=75.6\,$\mathrm{MeV}$$ and
$\gamma=8.16$, this reproduces the lattice result Borsanyi:2016ksw for the QCD
axion quite well at high temperatures $T\gtrsim 1\,$\mathrm{GeV}$$. By using
this ansatz we can calculate the trapping temperature via (2.28). The result
is

$\boxed{\frac{T_{\ast}}{\Lambda_{b,0}}\approx(2.12)^{\frac{\gamma}{6+\gamma}}\quantity(2\times
10^{8})^{\frac{2}{6+\gamma}}\quantity(\frac{g_{s}(T_{\ast})}{72})^{-\frac{2}{6+\gamma}}\quantity(\frac{\Lambda_{b,0}}{\rm{GeV}})^{\frac{2}{6+\gamma}}\quantity(\frac{h^{2}\Omega_{\phi,0}}{h^{2}\Omega_{\rm
DM}})^{-\frac{2}{6+\gamma}}.}$ (2.30)

We can observe that the trapping temperature does depend only on the zero
temperature barrier height, and on the scaling of the axion mass at high
temperature. For the QCD axion, the trapping temperature does not depend on
the axion mass, and it is given by

$T_{\ast}^{\rm
QCD}\approx\quantity(1.21\,$\mathrm{GeV}$)\quantity(\frac{h^{2}\Omega_{\phi,0}}{h^{2}\Omega_{\rm
DM}})^{-0.141}.$ (2.31)

We show a plot of the trapping temperature as a function of the zero-
temperature barrier height for various choices of $\gamma$ in Figure 2.
Equation (2.30) is a key quantity that determines the different fragmentation
regimes.

Figure 2: Trapping temperature $T_{\ast}$ given in Eq. 2.30 (left figure), and
the ratio between $T_{\ast}$ and the zero-temperature barrier height
$\Lambda_{b,0}$ (right plot) as a function of $\Lambda_{b,0}$ for different
choices of $\gamma$ which parametrizes the early-time scaling of the barrier
height, and thus of the axion mass, with temperature as defined in (2.29). We
see that this ratio increases for larger $\Lambda_{b,0}$, and the trapping
temperature becomes more sensitive to the early time behavior of the axion
potential. The dashed vertical purple line shows the zero-temperature QCD
barrier height $75.6\,$\mathrm{MeV}$$.

##  3 Analytical theory of parametric resonance in kinetic misalignment

###  3.1 Parametric resonance

After the discussion of the evolution of the homogeneous mode $\Theta$, now we
turn our attention to fluctuations. We generalise the results of
Fonseca:2019ypl by including the regime of fragmentation after trapping. We
also use different initial conditions for the fluctuations. If one neglects
the expansion of the universe, the equation of motion for axion mode functions
(2.7) takes the form

$\ddot{\theta}_{k}+\quantity[\frac{k^{2}}{a^{2}}+m^{2}\cos\Theta]\theta_{k}=0,$
(3.1)

where $a$ and $m$ are constants. This is a second-order differential equation
with periodic coefficients, so it has the form of the _Hill’s equation_
Hill1966 . According to the Floquet theorem, the solutions should be of the
form

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

where $\theta_{\pm}(t;k)$ are periodic functions in time, and $\mu_{k}$’s are
in general complex numbers known as _Floquet exponents_. If $\Re{\mu_{k}}>0$,
then the mode does grow exponentially during one oscillation, and one says
that the mode experiences instability via parametric resonance. The primary
goal of the _Floquet analysis_ is to determine the _instability bands_ of the
system, which are the regions where $\Re{\mu_{k}}>0$.

In the case of oscillations after trapping, a closed-form expression for the
Floquet exponents has been obtained in Greene:1998pb . We have found that this
method can easily be extended to derive an analogous expression for Floquet
exponents during the evolution before trapping. In this section we will just
state the result, and discuss the consequences. The detailed derivation in
presented in Appendix A.

Let $\varrho$ be the dimensionless energy of the homogeneous mode as defined
in (2.9), and let us introduce a dimensionless momentum

$\kappa\equiv k/(ma).$ (3.3)

Then the Floquet exponents have the following analytical solution:

$\mu_{k}=\sqrt{8\kappa^{2}\quantity(\varrho-\kappa)\quantity(1-\varrho+\kappa^{2})}\times\begin{cases}\mathcal{T}_{>}^{-1}(\varrho)\,\mathcal{I}\quantity(\arcsin\quantity(1/\sqrt{2\varrho-1})),&\varrho>1\\\
\mathcal{T}_{<}^{-1}(\varrho)\,\mathcal{I}\quantity(\pi/2),&\varrho<1\end{cases},$
(3.4)

The first line in (3.4) is new while the second line was derived in the
literature. $\mathcal{T}_{>}$ and $\mathcal{T}_{<}$ are the periods before and
after trapping respectively, and $\mathcal{I}(\varphi)$ is the following
integral:

$\mathcal{I}(\varphi)=\int_{0}^{\varphi}\frac{\differential{\vartheta^{\prime}}}{1+\quantity(1-2\varrho+2\kappa^{2})\sin^{2}\vartheta^{\prime}}\frac{2\sin^{2}\vartheta^{\prime}}{\sqrt{\quantity(1+\sin^{2}\vartheta^{\prime})\quantity[1+(1-2\varrho)\sin^{2}\vartheta^{\prime}]}}.$
(3.5)

The periods $\mathcal{T}$ can be derived from the conservation of energy. They
read

$\mathcal{T}_{>}(\varrho)=\frac{2}{m\sqrt{\varrho}}\textrm{K}(1/\sqrt{\varrho})\mbox{\quad
and\quad}\mathcal{T}_{<}(\varrho)=\frac{2}{m}\textrm{K}(\sqrt{\varrho}).$
(3.6)

From the square root term in (3.4) we can immediately obtain the instability
band. Before trapping it is given by

$\varrho-1<\kappa^{2}<\varrho,$ (3.7)

while its form after trapping is

$0<\kappa^{2}<\varrho.$ (3.8)

If one neglects the Hubble expansion, both $\varrho$ and $\kappa$ do not
change, so these modes keep growing. Therefore, in order to get realistic
results we need to incorporate the Hubble expansion, which we do next.

The Hubble expansion and the temperature dependence of the axion mass can be
incorporated by restoring the time dependence of $\varrho$ and $\kappa$ in the
expression for the Floquet exponent (3.4):

$\kappa\rightarrow\kappa(t)=\frac{k/a(t)}{m(t)}\mbox{\quad
and\quad}\varrho\rightarrow\varrho(t),$ (3.9)

where $\varrho(t)$ is given by (2.9). As a result, the Floquet exponents now
become time-dependent:

$\mu_{k}\rightarrow\mu_{k}\quantity(\kappa(t),\varrho(t)).$ (3.10)

Then we can obtain the total amplification factor $N_{k}$ of a given mode by
integrating the time-dependent Floquet exponents over time:

$N_{k}(t)=\exp(\int_{t_{i}}^{t}\differential{t^{\prime}}\mu\quantity(\kappa(t^{\prime}),\varrho(t^{\prime}))),$
(3.11)

where $t_{i}$ is some initial time before which the parametric resonance is
not effective. At this point it is convenient to introduce the dimensionless
quantities via

$\kappa_{\ast}\equiv\frac{k/a_{\ast}}{m_{\ast}}\mbox{\quad,\quad}\tau\equiv
2H_{\ast}t\mbox{\quad,\quad}\widetilde{\mu}\equiv\frac{\mu}{m_{\ast}},$ (3.12)

where $\ast$-subscript implies that the quantities are evaluated at $T_{\ast}$
defined by (2.8). In terms of these quantities (3.11) becomes

$\boxed{N_{k}\quantity(\tau)=\exp(\frac{m_{\ast}}{2H_{\ast}}\int_{\tau_{i}}^{\tau}\differential{\tau^{\prime}}\widetilde{\mu}\quantity(\kappa(\tau^{\prime}),\varrho(\tau^{\prime})))\equiv\exp(\frac{m_{\ast}}{2H_{\ast}}\mathcal{B}_{k}(\tau)).}$
(3.13)

The growth factor $\mathcal{B}_{k}$ does depend on the temperature scaling
$\gamma$ of the axion mass, but not on the model parameters such as $m_{\ast}$
and $H_{\ast}$. For the modes which amplify most efficiently its value is
$\sim\mathcal{O}(0.5)$. Therefore the efficiency of the fragmentation is
effectively determined by the hiearchy between the axion mass and Hubble rate
at trapping. This result can be understood physically. If the Hubble is much
smaller than the axion mass at the beginning of oscillations, the redshifting
of the homogenous mode becomes very slow allowing the axion to probe non-
quadratic parts of its potential for a longer time. In the Standard
Misalignment Mechanism, $m_{\ast}/H_{\ast}\sim 3$ so the expansion quickly
redshifts the amplitude of the oscillations which makes the parametric
resonance ineffective. In a nutshell, efficient parametric resonance requires
a mechanism which delays the onset of oscillations. The Kinetic Misalignment
Mechanism provides this via large initial kinetic energy. In the Large
Misalignment Mechanism Zhang:2017dpp ; Arvanitaki:2019rax this is achieved by
tuning the initial angle to the top of the axion potential such that the onset
of oscillations is delayed due to the small potential gradient at the top.

It is also instructive to study how the shape of the instability bands changes
with time. For this, let us introduce the following quantity:

$\widetilde{J}\quantity(t)=\frac{m_{\ast}}{m(t)}\quantity(\frac{a_{\ast}}{a(t)})^{3}.$
(3.14)

Then we can write $\kappa(\tau)$ as

$\kappa(\tau)=\kappa_{\ast}\widetilde{J}(\tau)\tau.$ (3.15)

To obtain $\varrho(\tau)$ we note that the adiabatic invariant (2.14) can be
expressed in terms of $\widetilde{J}$ by

$J(\varrho(\tau))=8f^{2}m(t)a^{3}(t)\widetilde{J}(\tau),$ (3.16)

so that

$\widetilde{J}(\varrho)=\begin{cases}\sqrt{\varrho}\,\textrm{E}(1/\sqrt{\varrho}),&\varrho>1\\\
\quantity(\varrho-1)\textrm{K}(\sqrt{\varrho})+\textrm{E}(\sqrt{\varrho}),&\varrho<1\end{cases}.$
(3.17)

At a given time this function can be inverted to get $\varrho(\tau)$. In
general this needs to be done numerically, however in asymptotic regimes we
can use the following approximations:

$\widetilde{J}(\varrho\gg 1)\approx\frac{\pi}{2}\sqrt{\varrho}\mbox{\quad
and\quad}\widetilde{J}(\varrho\ll 1)\approx\frac{\pi}{4}\varrho.$ (3.18)

We are now ready to discuss the behavior of the instability bands. Before
trapping, the expression for the instability band (3.7) can be written as

$\frac{\varrho-1}{\quantity(\widetilde{J}(\tau)\tau)^{2}}<\kappa_{\ast}^{2}<\frac{\varrho}{\quantity(\widetilde{J}(\tau)\tau)^{2}}.$
(3.19)

Initially $\varrho$ is large, so the width of the instability band is narrow.
As the energy of the homogeneous mode decreases, the instability band
gradually widens. The behavior after trapping is slightly more involved. The
expression for the instability band (3.8) becomes

$0<\kappa_{\ast}^{2}<\frac{\varrho}{\quantity(\widetilde{J}(\tau)\tau)^{2}}.$
(3.20)

At late times when $\varrho\ll 1$ we have $\varrho\approx 4\widetilde{J}/\pi$.
Therefore the upper limit of the instability band can be approximated as

$\frac{\varrho}{\quantity(\widetilde{J}(\tau)\tau)^{2}}\approx\frac{4}{\pi}\tau^{-1/2}\times\begin{cases}\tau^{\gamma/4},&T>T_{c}\\\
m_{0}/m_{\ast},&T<T_{c}\end{cases}.$ (3.21)

We see that for $\gamma>2$, the upper limit of the instability band grows with
time until the axion reaches its zero-temperature mass, and then decreases
with time. This means that the modes with $\kappa_{\ast}>1$ which did enter
and exit the instability band at the rolling stage will re-enter the
instability band after trapping. On the other hand, for $\gamma<2$, and also
for constant axion mass, the width of the instability band will shrink after
trapping. We show the evolution of the instability bands, Equations (3.19) and
(3.20), as a function of the scale factor for these two cases in Figure 3.

Figure 3: Time evolution of the instability bands, Eq. (3.19) and (3.20), as a
function of the scale factor for two benchmark points with $\gamma=8$ (upper
plot), and constant mass $\gamma=0$ (lower plot). In both cases we set
$m_{0}=10^{-11}\,\rm{eV}$, and choose the decay constant $f$ such that
$m_{\ast}/H_{\ast}=10^{2}$.

Even though the Floquet analysis predicts a very wide instability band, the
fragmentation is most efficient around trapping as can be observed in Figure
4. In this plot, we have used the parameters of Figure 3 with $\gamma=8$, but
we also show the value of the Floquet exponents, and zoom-in into the region
around trapping. The white lines denote the boundaries of the instability
band.

Figure 4: Evolution of the instability bands together with the value of the
Floquet exponents by zooming into the region around trapping. The white lines
denote the boundaries of the instability band. Fragmentation is most efficient
around trapping $a\approx a_{\ast}$. The model parameters are identical to the
ones used in Figure 3 with $\gamma=8$.

Our approximation so far sizably overpredicts the amplification factor. The
reason is that the Floquet analysis gives the amplification factor during one
oscillation by neglecting the damping of the mode due to the redshift. More
specifically, if the Floquet exponents for a mode are imaginary at all times,
our approximation would predict that the amplitude of this mode will stay
constant which is of course not correct. Instead, the amplitude of this mode
will get redshifted, and becomes smaller. To account for this, we assume that
on top of the amplification given by (3.13), all modes redshift like a free
particle. This decay factor can be calculated using WKB approximation and it
reads

$A_{k}(t)\propto\omega_{k}^{-1/2}(t)a^{-3/2}(t)\mbox{\quad
where\quad}\omega_{k}(t)=\sqrt{\frac{k^{2}}{a^{2}(t)}+m^{2}(T)}.$ (3.22)

So our final ansatz for the mode functions is

$\absolutevalue{\theta_{k}(t)}=\theta_{k,i}A_{k}(t)N_{k}(t),$ (3.23)

where $\theta_{k,i}$ is related to the initial field power spectrum (2.5) by
$\theta_{k,i}=\sqrt{P_{\theta}(k)}$, $A_{k}(t)$ is normalized such that it is
unity initially, and we have omitted the oscillatory term. We provide a
comparison between this ansatz and the full numerical solutions of the mode
functions in Figure 5. The parameters are the same as in Figures 3 and 4. The
thin solid lines show the numerical solutions, while the thick dashed lines
are calculated via (3.23). We started the numerical solution at $\tau=0.1$
with an initial amplitude $\theta_{k,i}=1$ for all the modes. We confirm that
our ansatz (3.23) provides a very reasonable approximation to the numerical
solution.

Figure 5: Comparison between the numerical solution (thin solid lines) and our
analytical ansatz (3.23) (thick dashed lines) for three benchmark momentum
modes, showing a very reasonable agreement. The parameters are the same as in
Figures 3 and 4.

###  3.2 An estimate of the back-reaction

In our discussion so far we have neglected the back-reaction of the
fluctuations on the homogeneous mode. However this approximation breaks down
when the fragmentation is efficient so that the fluctuations cannot be
considered small compared to the homogeneous mode Fonseca:2019ypl . If the
only thing we want to know is whether at a given point on the ALP parameter
space the relic density is dominated by fragmented axions or by the
homogeneous mode, then the back-reaction effects might not be important.
However, ultimately we are interested in observational consequences of
fragmentation such as miniclusters Eroncel:2022c . For this we need to
calculate the density power spectrum of the ALP field after fragmentation, and
in order to get accurate results we need to take the back-reaction into the
account.

Precise study of the back-reaction should be performed by using non-
perturbative methods such as lattice. However, we can obtain a simple estimate
semi-analytically. For this, we only need to assume that fragmentation steals
energy from the homogeneous mode, and gives it to the fluctuations. This
implies that the energy gained by the fluctuations during the fragmentation is
equal to the energy lost by the homogeneous mode during the same period:

$\Delta\rho_{\rm fluct}=-\Delta\rho_{\Theta},$ (3.24)

where $\rho_{\Theta}$ is the energy density in the homogeneous mode. This idea
has also been used in Fonseca:2019ypl for an extensive discussion of the
fragmentation before trapping. As we will show later, our calculations
reproduce these results under appropriate limits.

Figure 6: This plot compares the estimates of the amplification factor $N_{k}$
performed by including/not including the back-reaction effects through the
procedure (3.30). It shows the importance of back-reaction effects as $f$ is
decreased. This plot is made for a constant mass $m_{*}=m_{0}=10^{-15}$ eV.

The energy density in the fluctuations is

$\rho_{\rm
fluct}=\frac{f^{2}}{2}\int\frac{\differential[3]{k}}{(2\pi)^{3}}\quantity[\absolutevalue{\dot{\theta}_{k}(t)}^{2}+\quantity(\frac{k^{2}}{a^{2}}+m^{2}\cos\Theta)\absolutevalue{\theta_{k}(t)}^{2}],$
(3.25)

where we have neglected higher-order terms in $\theta_{k}$. Since the mode
functions oscillate very rapidly, we can assume that the kinetic and potential
terms contribute equally. We also average over the oscillations of the
homogeneous mode. With these assumptions and our ansatz (3.23), we approximate
the energy density in fluctuations as

$\rho_{\rm
fluct}(t)\approx\frac{f^{2}}{4\pi^{2}}\int\differential{k}k^{2}\quantity(\frac{k^{2}}{a^{2}}+m^{2}\overline{\cos\Theta})P_{\theta}(k)A_{k}^{2}(t)N_{k}^{2}(t),$
(3.26)

where

$\overline{\cos\Theta}\equiv\frac{1}{\Theta_{\rm max}-\Theta_{\rm
min}}\int_{\Theta_{\rm min}}^{\Theta_{\rm
max}}\differential{\Theta}\cos\Theta=\begin{cases}0,&\varrho>1,\\\
\operatorname{sinc}(2\arcsin\sqrt{\varrho}),&\varrho<1\end{cases}.$ (3.27)

The extra $1/2$ factor in (3.26) arises from averaging over the mode function
oscillations. From this result, we can estimate the energy lost by the
homogeneous mode during a short time period $\Delta t$ as

$-\frac{\Delta\rho_{\Theta}}{\Delta t}=\frac{\Delta\rho_{\rm fluct}}{\Delta
t}\approx\frac{f^{2}}{4\pi^{2}}\int\differential{k}k^{2}\quantity(\frac{k^{2}}{a^{2}}+m^{2}\overline{\cos\Theta})P_{\theta}(k)A_{k}^{2}(t)\frac{N_{k}^{2}(t)-N_{k}^{2}(t-\Delta
t)}{\Delta t}.$ (3.28)

Note that we have only varied the amplification factor $\mathcal{N}_{k}$ since
we are only interested in the change in the energy due to the fragmentation.
In the $\Delta t\rightarrow 0$ limit we find

$-\derivative{\rho_{\Theta}}{t}=\derivative{\rho_{\rm
fluct}}{t}\approx\frac{f^{2}}{4\pi^{2}}\int\differential{k}k^{2}\quantity(\frac{k^{2}}{a^{2}}+m^{2}\overline{\cos\Theta})P_{\theta}(k)A_{k}^{2}(t)\quantity[2N_{k}^{2}(t)\mu\quantity(\kappa(t),\varrho(t))],$
(3.29)

where we have used (3.11) when taking the time derivative.

With this result, we can employ the following procedure to calculate the
amplification factors including back-reaction. Choose an initial time
$\tau_{i}$ at which $P_{\theta}(k)$ is known and $N_{k}(\tau_{i})=1$ for all
modes. Pick a sufficiently small time-step $\Delta\tau$, and at each time
interval calculate the amplification factors by

$\ln N_{k}(\tau+\Delta\tau)=\ln
N_{k}(\tau)+\frac{m_{\ast}}{2H_{\ast}}\int_{\tau}^{\tau+\Delta\tau}\differential{\tau}^{\prime}\widetilde{\mu}_{k}(\tau^{\prime})\approx\ln
N_{k}(\tau)+\frac{m_{\ast}}{2H_{\ast}}\Delta\tau\,\widetilde{\mu}_{k}(\tau).$
(3.30)

For each time step, find $\varrho$ via (3.14) and (3.17) and therefore the
time-dependent Floquet exponent $\mu_{k}$. Take into account the decrease in
$\varrho$ due to the back-reaction via (3.29). Then repeat the steps either
until fragmentation is no longer efficient, or when $\varrho$ drops to zero
indicating complete fragmentation. At this stage, the dynamics of the system
becomes non-perturbative, so our calculation should be considered an estimate.
We illustrate the importance of back-reaction effects in Figure 6. We will use
this procedure in the following, especially to determine the boundary between
the ‘after-trapping’ and ‘before-trapping’ fragmentation regimes defined in
Section 6.

###  3.3 Initial conditions for the mode functions

In order to study the consequences of fragmentation for ALP dark matter, we
need to specify the initial conditions for the mode functions prior to
fragmentation. A scalar field, like any component in the universe, can have
adiabatic and isocurvature fluctuations. Adiabatic ones arise solely due to
the temperature fluctuations in the universe, and are unavoidable for any
cosmological fluid including scalar fields. On the other hand, the
isocurvature fluctuations are the fluctuations on constant density slices, and
they describe any other kind of fluctuation which is not adiabatic.

In the case of standard misalignment mechanism, adiabatic fluctuations of an
ALP field are negligible during the early evolution when it is Hubble frozen
Marsh:2015xka . They start to grow once the homogeneous mode starts
oscillating. If the ALP field is present during inflation, for example if the
Peccei-Quinn (PQ) symmetry is broken during the inflation in the case of QCD
axion, then it will pick up quantum fluctuations given by
$\delta\phi=H_{I}/2\pi$, where $H_{I}$ is the inflation scale. This causes
fluctuations in the initial axion angle $\Theta_{i}$ given by
$\delta\Theta_{i}\simeq H_{I}/(2\pi S_{I})$ where $S_{I}$ is the effective
axion decay constant during inflation which can be different than its value at
low temperature Linde:1991km .

The situation changes drastically in the case of rotating axions. In this
case, the ALP field has a much more complicated cosmological history which is
highly model-dependent. A common property in all of these histories is that
the ALP field starts moving much earlier compared to the standard case, and it
has a large velocity long before the ALP potential turns on. As we show below,
this velocity together with curvature perturbations acts as a source term for
adiabatic fluctuations and make them grow. In addition to these, the kick
mechanism can induce isocurvature fluctuations from the quantum fluctuations
of the field responsible for producing the kick. These can constrain the
parameter space due to their isocurvature nature, and can also lead to domain
wall problems in some cases Co:2020dya .

An extensive discussion of model realizations of the kinetic misalignment
mechanism including the possible cosmological histories as well as the
implications of isocurvature fluctuations was presented in Gouttenoire:2021jhk
. The detailed consequences of ALP dark matter in these UV completions will be
presented in Eroncel:2022a . For simplicity, and to remain agnostic about the
UV completion, we consider a standard cosmological history in this work, and
assume that the modes that are relevant for fragmentation do enter the horizon
when the ALP field scales as kination. We will present the necessary condition
for the validity of this assumption towards the end of this sub-section.

To get the evolution, we start by writing the FRLW metric with curvature
perturbations included. In conformal time $\eta$ and Newtonian (conformal)
gauge, the metric reads555A brief review of cosmological perturbation theory
and our notation convention can be found in Appendix B.

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

Here $\Psi$ and $\Phi$ are Bardeen potentials Bardeen:1980kt , and they
represent the curvature perturbations. In Appendix B.2, we show that the mode
functions have the following equations of motion at early times:

$\theta_{k}^{\prime\prime}+2\mathcal{H}\theta_{k}^{\prime}+k^{2}\theta_{k}=-4\Phi_{k}^{\prime}\Theta^{\prime},$
(3.32)

where $\prime$ denotes the derivative with respect to conformal time, and
$\mathcal{H}=aH$ is the conformal Hubble parameter. The right side of this
equation represents the source term due to the curvature fluctuations. This
term clearly represents the difference between the standard and kinetic
misalignment mechanisms. In the standard case, the ALP field is frozen due to
the Hubble friction, so $\Theta^{\prime}=0$. Therefore, the source term is
absent, and the adiabatic fluctuations remain zero until the oscillations
start. However, in the case of kinetic misalignment mechanism, the ALP field
receives a kick at a much earlier time, so the right side is non-zero for a
much longer time. As a result, the standard and kinetic misalignment
mechanisms predict different initial conditions for mode functions. To our
knowledge, this is the first time that this is pointed out. We will discuss
more thoroughly further implications of this source term in Eroncel:2022b .

By an explicit calculation of the mode functions which we describe in Appendix
B.2, see (Eq. B.53), we obtain the following result for the power spectrum:

$\boxed{P_{\theta}(k;a)=\absolutevalue{\theta_{k}(a)}^{2}\approx\frac{2\pi^{2}}{k^{3}}\quantity(\frac{1}{3})^{2}A_{s}\quantity(\frac{\dot{\Theta}}{H})^{2}}$
(3.33)

where $A_{s}=2.1\times 10^{-9}$ is the amplitude of the primordial power
spectrum at the pivot scale $k=0.05\,\textrm{Mpc}^{-1}$ as measured by Planck
2018 (TT,TE,EE+lowE+lensing 68%) Aghanim:2018eyx . This result is valid for
both super- and sub-horizon modes, but assumes that the mode is super-horizon
when the ALP field starts its kination-like scaling which we denote by $a_{\rm
kin}$. The behavior of the modes which were sub-horizon at $a_{\rm kin}$
cannot be determined without specifying the cosmological history before
$a_{\rm kin}$.

From this result alone, we can put a bound on the duration of the kination-
like scaling. For this, we evaluate the variance of the axion velocity:

$\expectationvalue{\quantity(\delta\dot{\theta})^{2}}=\frac{1}{2\pi^{2}}\int_{0}^{\infty}\differential{k}k^{2}\absolutevalue{\dot{\theta}_{k}}^{2}=\frac{1}{2\pi^{2}a^{2}}\int_{0}^{\infty}\differential{k}k^{2}\absolutevalue{\theta_{k}^{\prime}(\eta)}^{2}.$
(3.34)

The integral is dominated by the modes which are sub-horizon, but were super-
horizon at $a_{\rm kin}$. Then by approximating
$\absolutevalue{\theta_{k}^{\prime}(\eta)}^{2}\approx
k^{2}\absolutevalue{\theta_{k}(\eta)}^{2}$, and using (3.33) we find at late
times

$\expectationvalue{\quantity(\delta\dot{\theta})^{2}}\approx\frac{1}{2}\quantity(\frac{1}{3})^{2}A_{s}\dot{\Theta}^{2}\quantity(\frac{a}{a_{\rm
kin}})^{2}.$ (3.35)

Then we can estimate the density contrast by

$\delta_{\phi}\sim\frac{\delta\quantity(\dot{\theta}^{2}/2)}{\dot{\Theta}^{2}/2}\sim\frac{2}{\dot{\Theta}}\sqrt{\expectationvalue{\quantity(\delta\dot{\theta})^{2}}}\approx\sqrt{2A_{s}}\quantity(\frac{1}{3})\frac{a}{a_{\rm
kin}}.$ (3.36)

This becomes $\mathcal{O}(1)$ when

$\boxed{a/a_{\rm kin}\gtrsim 10^{5}.}$ (3.37)

This implies that if the ALP field scales as kination more than
$\ln(10^{5})\sim 10$ number of e-folds, then the ALP field cannot be
considered homogeneous anymore. A similar bound applies if the ALP drives a
kination era. The bound (3.37) is an important constraint when considering UV
completions Eroncel:2022a , in particular models where the ALP drives
temporarily a kination era in the early universe, enhancing primordial GW
signals. Such bound on the total duration of a kination era was overlooked in
previous literature on kination, as discussed in Gouttenoire:2021jhk . We will
investigate this issue further in an upcoming work Eroncel:2022b . In this
paper, we do assume that this bound is not violated, so that the ALP field can
be considered homogeneous at the onset of fragmentation.

As we have stated previously, in this work we will assume that all the modes
relevant for fragmentation are super-horizon at $a_{\rm kin}$. This assumption
requires that

$\frac{k}{a_{\rm kin}H_{\rm
kin}}<1\quad\Rightarrow\quad\kappa_{\ast}\frac{m_{\ast}}{H_{\ast}}\frac{a_{\rm
kin}}{a_{\ast}}<1.$ (3.38)

If we demand that this condition is satisfied for the relevant modes for
fragmentation, i.e. $\kappa_{\ast}\sim\mathcal{O}(1)$, then we need to demand
that

$\frac{m_{\ast}}{H_{\ast}}\lesssim\frac{a_{\ast}}{a_{\rm kin}}.$ (3.39)

The quantity $a_{\ast}/a_{\rm kin}$ depends on the specifics of the UV
completion, however we have shown that the homogeneity of the ALP field prior
to fragmentation requires $a_{\ast}/a_{\rm kin}\gtrsim 10^{5}$ which puts the
bound $m_{\ast}/H_{\ast}\lesssim 10^{5}$. In the next section, we will show
that this bound is satisfied in the region of the ALP parameter space where
our analytical approximation is under control.

With all these assumptions above, we can fix the field power spectrum at the
onset of fragmentation by

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

As we will see explicitly in the coming sections, the choice of $a_{i}$ is not
relevant for the final result, as long as it is early enough that
$\dot{\Theta}\propto a^{-3}$, and any fragmentation prior to $a_{i}$ can be
neglected.

##  4 ALP dark matter from fragmentation

###  4.1 Overview of the fragmentation regions

One can split the ALP parameter space into different regions depending whether
the fragmentation happens and when. There are four different scenarios:

  1. 1.

Standard misalignment: In this case the onset of oscillations is not delayed
from its conventional value $m(T_{\rm osc})=3H(T_{\rm osc})$, and the standard
misalignment mechanism is at play.

  2. 2.

Kinetic misalignment with weak fragmentation: The ALP field has a non-zero
initial velocity such that the onset of oscillations is delayed, but the
particle production is not strong enough so that the energy density in
fluctuations is always subdominant compared to the zero mode.

  3. 3.

Fragmentation after trapping: The ALP field is completely fragmented, but the
fragmentation ends after it would have been trapped by the potential in the
absence of fragmentation, i.e. $T_{\rm end}<T_{\ast}$.

  4. 4.

Fragmentation before trapping: The ALP field is completely fragmented, and the
fragmentation ends before it would have been trapped by the potential in the
absence of fragmentation. It other words the fragmentation is complete at
$T_{\rm end}>T_{\ast}$ where $T_{\ast}$ is given by (2.31).

The boundaries between these regions depend strongly on the hierarchy between
the axion mass and Hubble at trapping, and have a very mild dependence on the
other model parameters. In the rest of this section, we will calculate these
boundaries while giving details on the properties of the fragmentation in each
region.

An overview of these regions on the $m_{0}$–$f$ plane along with various model
parameters can be seen in Figures 7 and 8.

Figure 7: Overview of the fragmentation regions for temperature-dependent
axion mass with $\gamma=8$ (left plot), and for constant axion mass (right
plot). The solid contours denote the zero-temperature barrier heights, while
the dashed ones are the $m_{\ast}/3H_{\ast}$ contours, where $T_{\ast}$ is
given in Eq. 2.30. Figure 8: Same as Figure 7, but the dashed contours now
show the trapping temperature $T_{*}$ in GeV, as given in Eq. 2.30.

####  4.1.1 Boundary between regions $1$ and $2$:

In order to have kinetic misalignment, the onset of oscillations needs to be
delayed so $m(T_{\ast})>3H(T_{\ast})$. Therefore the boundary between the
regions $1$ and $2$ is given by the condition

$m(T_{\ast})=3H(T_{\ast}).$ (4.1)

In a generic ALP model, the critical value of the axion decay constant $f_{\rm
crit}$ for a given zero-temperature axion mass $m_{0}$ can be calculated via
(2.28) and (2.29). It can approximately be expressed as

$\boxed{f_{\rm crit}^{1\text{ -- }2}\approx
10^{12}\,$\mathrm{GeV}$\quantity[\quantity(3\pi)^{-12-2\gamma}\quantity(\frac{g_{s}(T_{\ast})}{72})^{2+\gamma}\quantity(\frac{m_{0}}{$\mathrm{eV}$})^{-4-\gamma}\quantity(\frac{h^{2}\Omega_{\Theta,0}}{h^{2}\Omega_{\rm
DM}})^{8+2\gamma}]^{\frac{1}{16+3\gamma}}}$ (4.2)

For the QCD axion this value is

$f_{\rm crit}^{1\text{ -- }2}\approx 2.15\times
10^{11}\;$\mathrm{GeV}$,\mbox{\quad QCD Axion.\quad}$ (4.3)

This is in agreement with the result obtained in Co:2019jts .

####  4.1.2 Boundary between regions $2$ and $3$:

To calculate this boundary, we need to find the region where the energy
density in the fluctuation remains subdominant compared to the energy density
in the homogeneous mode. In this region, the back-reaction can be neglected.
The energy density in the homogeneous mode is

$\rho_{\Theta}=2f^{2}m^{2}\varrho,$ (4.4)

while the energy density in the fluctuations is (3.26)

$\rho_{\rm
fluct}=\frac{f^{2}}{4\pi^{2}}\int\differential{k}k^{2}\quantity(\frac{k^{2}}{a^{2}}+m^{2}\overline{\cos\Theta})P_{\theta}(k)A_{k}^{2}(t)N_{k}^{2}(t).$
(4.5)

Let us define a quantity $\Delta\equiv\lim_{t\rightarrow\infty}\rho_{\rm
fluct}/\rho_{\Theta}$ as the ratio between the two energy densities at late
time limit. We refer this quantity as the “efficiency” of the fragmentation.
Then the transition from weak fragmentation to complete fragmentation should
occur approximately when $\Delta$ reaches to unity. This point will determine
the boundary between the regions $2$ and $3$.

We are interested in the late time limit of $\Delta$. In this limit
$\overline{\cos\Theta}=1$, and we can assume that all the modes which did grow
during the parametric resonance become non-relativistic so that the momentum
term in (4.5) can be ignored. The redshift factor $A_{k}^{2}(t)$ at late times
can be approximated by

$A_{k}^{2}(t)=\frac{\omega_{k}(t_{i})}{\omega_{k}(t)}\quantity(\frac{a_{i}}{a})^{3}\approx\frac{k/a_{i}}{m_{0}}\quantity(\frac{a_{i}}{a})^{3}=\frac{ka_{i}^{2}}{m_{0}a^{3}}.$
(4.6)

By plugging this result into (4.5) and using the expression for the initial
power spectrum, we find that the energy density in the fluctuations approach
an asymptotic value given by

$\rho_{\rm
fluct}\rightarrow\frac{m_{0}m_{\ast}f^{2}}{2}\quantity(\frac{1}{3})^{2}A_{s}\quantity(\frac{\dot{\Theta}_{i}}{H_{i}}\frac{a_{i}}{a_{\ast}})^{2}\quantity(\frac{a_{\ast}}{a})^{3}\int\differential{\kappa_{\ast}}\exp(\frac{m_{\ast}}{H_{\ast}}\mathcal{B}_{\kappa}^{\infty}),$
(4.7)

where $\mathcal{B}_{\kappa}^{\infty}$ is the asymptotic value of the growth
factor defined in (3.13). The late time limit of the energy density of the
homogeneous mode can be found by using the fact that for $\varrho\ll 1$ we
have

$\varrho\approx\frac{4}{\pi}\widetilde{I}=\frac{4}{\pi}\frac{m_{\ast}}{m}\quantity(\frac{a_{\ast}}{a})^{3}.$
(4.8)

With this result, the efficiency factor takes the form

$\Delta=\frac{\pi}{16}\quantity(\frac{1}{3})^{2}A_{s}\quantity(\frac{\dot{\Theta}_{i}}{H_{i}}\frac{a_{i}}{a_{\ast}})^{2}\int\differential{\kappa_{\ast}}\exp(\frac{m_{\ast}}{H_{\ast}}\mathcal{B}_{\kappa}^{\infty}).$
(4.9)

Note that at early time $\dot{\Theta}\propto a^{-3}$, therefore the pre-factor
in this expression does not depend on the choice of the initial time provided
that it is early enough so that $\dot{\Theta}\propto a^{-3}$, and any
fragmentation prior to this time can be neglected. By expressing the axion
velocity $\dot{\Theta}_{i}$ in terms of the yield parameter introduced in
Section 2.1 as $\dot{\Theta}_{i}=s(T_{i})Y/f^{2}$ where $s(T_{i})$ is the
entropy density at $T_{i}$, we can derive the following result:

$\frac{\dot{\Theta}_{i}}{H_{i}}\frac{a_{i}}{a_{\ast}}\approx\frac{4\pi}{3}\sqrt{\frac{g_{s}(T_{\ast})}{10}}\frac{YT_{\ast}M_{\rm
pl}}{f^{2}}.$ (4.10)

When arriving this, we have neglected the difference between the effective
degrees of freedom at $T_{i}$ and $T_{\ast}$, and also assume that the
effective degrees of freedoms in the entropy and the energy density are the
same. So our final expression for the efficiency factor is

$\Delta=\frac{\pi^{3}}{10}\quantity(\frac{1}{9})^{2}A_{s}g_{s}(T_{\ast})\quantity(\frac{YT_{\ast}M_{\rm
pl}}{f^{2}})^{2}\int\differential{\kappa_{\ast}}\exp(\frac{m_{\ast}}{H_{\ast}}\mathcal{B}_{\kappa}^{\infty}).$
(4.11)

With this factor, we can calculate the critical point at which the efficiency
becomes unity. We find that across the parameter space the transition happens
around $m_{\ast}/H_{\ast}\approx 42$ for $\gamma=8$, and around
$m_{\ast}/H_{\ast}\approx 38$ for constant axion mass ($\gamma=0$). This leads
to a critical axion decay constant given by

$\boxed{f_{\rm crit}^{2\text{ -- }3}\approx
10^{12}\,$\mathrm{GeV}$\quantity[\quantity(40\pi)^{-12-2\gamma}\quantity(\frac{g_{s}(T_{\ast})}{72})^{2+\gamma}\quantity(\frac{m_{0}}{$\mathrm{eV}$})^{-4-\gamma}\quantity(\frac{h^{2}\Omega_{\Theta,0}}{h^{2}\Omega_{\rm
DM}})^{8+2\gamma}]^{\frac{1}{16+3\gamma}}}$ (4.12)

For the QCD-Axion, this boundary corresponds to a decay constant of

$f_{\rm crit}^{2\text{ -- }3}\approx 1.5\times
10^{10}\;$\mathrm{GeV}$,\mbox{\quad QCD Axion.\quad}$ (4.13)

####  4.1.3 Boundary between regions 3 and 4:

In both regions the parametric resonance is efficient enough for complete
fragmentation. We want to determine whether the fragmentation ends before the
ALP field would have been trapped by the potential in the absence of
fragmentation. At early times when $\varrho\gg 1$, we can simplify the
fragmentation calculation to understand qualitative features of the
fragmentation. In this regime, we can also show that the procedure we have
outlined in Section 3 reproduces the results of Fonseca:2019ypl . We present
this calculation in Appendix C. Here we will give a quick summary.

In the $\varrho\gg 1$ regime, the exponential particle production happens when

$H(T)\lesssim\frac{\pi m^{4}(T)}{4\dot{\Theta}^{3}(T)}.$ (4.14)

The temperature at which this inequality is satisfied can be defined as the
onset of fragmentation provided that it is larger than the trapping
temperature $T_{\ast}$. We denote this temperature by $T_{\rm frag}$. Since
the back-reaction is not efficient in the beginning of fragmentation, the
onset of fragmentation does not depend on the initial conditions of the mode
functions.

The back-reaction becomes relevant when

$H(T)\lesssim\frac{\pi
m^{4}(T)}{4\dot{\Theta}^{3}(T)}\frac{1}{\ln\alpha^{-1}(T)},$ (4.15)

where for adiabatic initial conditions

$\alpha(T)=\quantity(\frac{4\pi}{9})^{2}\quantity(\frac{g_{s}(T_{\ast})}{80})^{2}\quantity(\frac{YT_{\ast}M_{\rm
pl}}{f^{2}})^{2}A_{s}\quantity(\frac{T}{T_{\ast}})^{2}.$ (4.16)

If the inequality (4.15) gets saturated before trapping, then the energy
density of the homogeneous mode gets diluted quite rapidly, and the
fragmentation is completed in a short amount of time. To estimate the
temperature at which the back-reaction becomes prominent, we can replace
$\dot{\Theta}(T)\approx s(T)Y/f^{2}$ and
$m(T)=m_{\ast}(T_{\ast}/T)^{\gamma/2}$ in (4.15) to get an equation in terms
of the temperature $T$ and other model parameters. Let $T_{\rm s}$ denote the
solution of this equation. For fragmentation to be completed before trapping,
this temperature $T_{\rm s}$ should be larger than the trapping temperature
$T_{\ast}$. Therefore, the Region 4 can be approximately defined by the
condition $T_{\rm s}\gtrsim T_{\ast}$. Across the parameter space, the
transition happens around $m_{\ast}/H_{\ast}\sim 9(5)\times 10^{2}$ for
$\gamma=8(0)$. The expression for the critical decay constant can be
approximated by

$\boxed{f_{\rm crit}^{3\text{ -- }4}\approx
10^{12}\,$\mathrm{GeV}$\quantity[\quantity(\mathcal{O}(1)\pi\times
10^{2})^{-12-2\gamma}\quantity(\frac{g_{s}(T_{\ast})}{72})^{2+\gamma}\quantity(\frac{m_{0}}{$\mathrm{eV}$})^{-4-\gamma}\quantity(\frac{h^{2}\Omega_{\Theta,0}}{h^{2}\Omega_{\rm
DM}})^{8+2\gamma}]^{\frac{1}{16+3\gamma}}}$ (4.17)

where the $\mathcal{O}(1)$ factor is $\gamma$-dependent. For the QCD axion,
the critical decay constant at this boundary is

$f_{\rm crit}^{1\text{ -- }2}\approx 7.3\times
10^{8}\;$\mathrm{GeV}$,\mbox{\quad QCD Axion.\quad}$ (4.18)

####  4.1.4 Breakdown of perturbativity for larger $m_{\ast}/H_{\ast}$

Finally we discuss what happens in the parameter space where
$m_{\ast}/H_{\ast}$ is even larger. For this we revisit the equation of motion
for the homogeneous mode (2.6), but we also include the back-reaction at
leading order. Then the equation of motion is modified to Fonseca:2019ypl

$\ddot{\Theta}+3H\dot{\Theta}+m^{2}\sin\Theta-\frac{1}{2}m^{2}\sin\Theta\expectationvalue{\quantity(\delta\theta)^{2}}=0.$
(4.19)

The last term is responsible for the back-reaction of fluctuations, and as a
result for the completion of the fragmentation. All of our analysis depends on
the assumption that this term is negligible before the onset of fragmentation.
Therefore we should check whether the adiabatic initial conditions for the
mode functions (3.33) do not violate this assumption at the beginning of
fragmentation.

By using (3.33), the variance can be calculated by

$\expectationvalue{\quantity(\delta\theta)^{2}}=\int\frac{\differential[3]{k}}{(2\pi)^{3}}\absolutevalue{\theta_{k}}^{2}\sim\quantity(\frac{1}{3})^{2}A_{s}\quantity(\frac{\dot{\Theta}}{H})^{2},$
(4.20)

where we have neglected the modes with $k/a_{\rm kin}H_{\rm kin}>1$, and
assumed $\ln(a/a_{\rm kin})\sim\mathcal{O}(1)$. At early times, the background
evolution is dominated by the first two terms in (4.19) since the mass is
negligible. However, if the variance becomes much larger than unity at these
times, it can affect the slow-roll evolution much earlier than the
fragmentation does. To avoid this case, we demand that

$3H(T)\dot{\Theta}(T)>\frac{m^{2}(T)}{2}\quantity(\frac{1}{3})^{2}A_{s}\quantity(\frac{\dot{\Theta}(T)}{H(T)})^{2}$
(4.21)

until the onset of fragmentation $T_{\rm frag}$ given by the solution of
(4.14). This puts a strong bound on the parameter space which can roughly be
approximated by

$\boxed{\frac{m_{\ast}}{H_{\ast}}\lesssim\mathcal{O}(1\text{---}10)\times
10^{3},}$ (4.22)

where the factor depends on $\gamma$. We want the stress that this result does
not exclude the parameter space where this bound is violated. It just implies
that the analysis we describe might not be reliable, and a dedicated study is
needed. In this paper, we concentrate on the region of the parameter space
where (4.22) is satisfied666We have also checked whether the variance term can
dominate the mass term in (4.19) before the onset of fragmentation. This
yields a weaker bound given by $m_{\ast}/H_{\ast}\lesssim\mathcal{O}(1)\times
10^{5}$..

###  4.2 ALP relic density with fragmentation

If the ALP field is completely fragmented, then all the energy density in the
homogeneous mode gets transferred into the fluctuations. This will have an
effect on the relic density today, since the redshift of the ALP energy
density is not necessarily the same with and without fragmentation. Naive
expectation is that the energy density will be diluted slightly since the
modes that are enhanced exponentially are mildly relativistic right after the
fragmentation. Therefore it is natural to ask whether this effect is
significant or not.

Figure 9: Left: Dilution factor (4.24) as a function of the axion decay
constant for an ALP with constant mass $m=10^{-15}\,$\mathrm{eV}$$. Right:
Amplification factor $N_{k}$ (3.13) as a function of the dimensionless
momentum $k/a_{\ast}m_{\ast}$. The solid lines show the spectrum at the end of
fragmentation, while the thin dashed lines together with the shaded regions
show the spectrum when the ALP is trapped. The line colors match with the
colored dots on the left plot.

At late times, the ALP energy density is given by

$\rho_{\theta}=\overline{\rho}_{\theta}+\frac{f^{2}}{2}\int\frac{\differential[3]{k}}{(2\pi)^{3}}\quantity[\dot{\theta}_{k}^{2}+\quantity(\frac{k^{2}}{a^{2}}+m^{2}_{0})\theta_{k}^{2}].$
(4.23)

where $\overline{\rho}_{\theta}$ is the energy density remaining in the
homogeneous mode which is negligible in the case of complete fragmentation.
Let $\rho_{\Theta}$ be the energy density without fragmentation. We define the
dilution factor as

$\mathcal{Z}\equiv\frac{\rho_{\theta}}{\rho_{\Theta}}.$ (4.24)

which will be relevant for the section on the gravitational-wave signal.
Unfortunately, it is not possible to calculate this factor precisely without a
proper lattice simulation. In the non-linear regime after the fragmentation,
the self-interactions between the enhanced momentum modes can modify the
momentum spectrum of the fluctuations which can cause $\mathcal{O}(1)$
modifications in the dilution factor. We leave the careful study of this
dilution factor with lattice simulations for future work.

Figure 10: Evolution of the energy density of the axion in the large
misalignment case (where the field is initially frozen and then redshifts as
matter once oscillations turn on) compared to the kinetic
misalignment/fragmentation case ($a^{-6}$ redshifting preceeding the matter
era, with almost no time for a radiation-like $a^{-4}$ behaviour in between).

At least we show in Fig. 9 our estimate for $\mathcal{Z}$. Here, we assume an
ALP with constant mass, fixed it to be $m=10^{-15}\,$\mathrm{eV}$$, and show
the behavior of the dilution factor as a function of the decay constant $f$ on
the left plot. For large decay constants, the fragmentation is inefficient,
therefore we have $\mathcal{Z}\approx 1$. As the fragmentation becomes more
efficient, i.e. for smaller decay constants, more and more energy is
transferred to the modes that are mildly relativistic after the fragmentation.
As a result, the dilution factor decreases confirming our expectations.
However, we observe the opposite trend at even lower decay constants where our
calculation predicts dilution factors larger than one, meaning that the relic
density is _enhanced_ by the fragmentation. This might be seen counter-
intuitive at first, but it is not. For small decay constants, most of the
fragmentation occurs before trapping as we can see from the right plot of Fig.
9. Recall that the energy density of the homogeneous mode scales as
$\rho_{\Theta}\propto a^{-6}$ in this regime. During fragmentation, the energy
density is transferred to the fluctuations that redshift slower compared to
$a^{-6}$. This yields to an enhancement in the relic density. Even though the
behavior of the dilution factor can be understood for large and small decay
constants, we cannot derive conclusions about the behavior in between with our
simplified calculation. The self-interactions between the excited modes can
modify the momentum distribution and therefore the dilution factor. This
effect cannot be captured by our approximation and require a lattice
calculation.

We also show the evolution of the energy density of the axion in figure 10.
This clearly shows that there is almost no time for a radiation-like equation
of state. The produced axions are not much relativistic when they are produced
and very quickly cool down so we transit quickly from a kination-like equation
of state to a matter-like equation of state for the axion.

###  4.3 Other initial conditions

Before concluding this section, we briefly comment on the modifications of our
results, if we had taken different initial conditions for the mode functions.
In particular, we consider the case where the mode functions are in the Bunch-
Davies vacuum initially. This implies that

$\evaluated{\delta\phi_{k}}_{i}=f\evaluated{\theta_{k}}_{i}=\frac{e^{-i(k/a_{i})t}}{a_{i}\sqrt{2k}}\quad\Rightarrow\quad
P(k)=\absolutevalue{\theta_{k}}^{2}=\frac{1}{2ka_{i}^{2}f^{2}}.$ (4.25)

This is the only modification we need to make. Then we can rederive all the
results we obtained in this section by using this power spectrum instead of
(3.40). We give a summary of modifications by skipping the details of the
calculation.

The efficiency factor $\Delta$ with which we have determined the boundary
between the regions 3 and 4 takes the form

$\Delta=\frac{1}{64\pi}\quantity(\frac{m_{\ast}}{f})^{2}\int\differential{\kappa}_{\ast}\kappa_{\ast}^{2}\exp(\frac{m_{\ast}}{H_{\ast}}\mathcal{B}_{\kappa}^{\infty}),\mbox{\quad
Bunch-Davies\quad}$ (4.26)

Note that the factor in front of the integral is much more suppressed compared
to factor with adiabatic initial conditions (4.11). Therefore larger
$m_{\ast}/H_{\ast}$ values are needed for complete fragmentation. We find that
$\Delta$ exceeds unity when $m_{\ast}/H_{\ast}\sim\mathcal{O}(1)\times 10^{2}$
where exact values show larger variations across the parameter space compared
to the adiabatic case. For the QCD axion, the critical decay constant is

$f_{\rm crit}^{2\text{ -- }3}\approx 2.1\times
10^{9}\;$\mathrm{GeV}$,\mbox{\quad QCD Axion, Bunch-Davies vacuum\quad},$
(4.27)

which is almost an order of magnitude larger than the critical value in the
adiabatic case (4.13).

In the $\varrho\gg 1$ regime, the fragmentation condition

$H(T)\lesssim\frac{\pi m^{4}(T)}{4\dot{\Theta}^{3}(T)}$ (4.28)

is not modified since it is independent of the initial conditions. The
condition for strong back-reaction does have a similar form

$H(T)\lesssim\frac{\pi
m^{4}(T)}{4\dot{\Theta}^{3}(T)}\frac{1}{\ln\alpha^{-1}_{\rm BD}(T)},$ (4.29)

except the expression for $\alpha(T)$ is modified to

$\alpha_{\rm BD}(T)=\frac{\dot{\Theta}^{2}(T)}{128\pi^{2}f^{2}}.$ (4.30)

##  5 Constraints on the ALP parameter space due to fragmentation

We need to impose a number of constraints to ensure that the fragmentation
process does not spoil existing cosmological observations. These are the
following777Note that axion fragmentation does not generate a domain wall
problem Morgante:2021bks .:

###  5.1 Structure formation

If all of the dark matter is made up of fragmented axions, they will need to
be sufficiently cold at matter-radiation equality in order to be consistent
with structure formation. This puts a constraint on the axion velocity $v_{\rm
eq}$ at matter-radiation equality $a_{\rm eq}$. In Section 3, we learned that
the momentum modes which grow most efficiently due to the parametric resonance
are those with $k\sim a_{\ast}m_{\ast}$. Then we can estimate the axion
velocity at matter-radiation equality as

$v_{\rm eq}\sim\frac{k/a_{\rm eq}}{m_{0}}\sim\quantity(\frac{a_{\ast}}{a_{\rm
eq}})\quantity(\frac{m_{\ast}}{m_{0}}).$ (5.1)

To get precise constraints from structure formation would require involved
numerical simulations, which is beyond the scope of this work. Instead, we
will use the bound $v_{\rm eq}<10^{-3}$ which is commonly considered in the
literature Berges:2019dgr . So we demand

$\quantity(\frac{a_{\ast}}{a_{\rm
eq}})\quantity(\frac{m_{\ast}}{m_{0}})\lesssim 10^{-3}.$ (5.2)

###  5.2 Big Bang Nucleosynthesis (BBN)

The presence of an additional energy component at the time of BBN will modify
the expansion rate and therefore be subject to constraints from observations
of primordial helium-4 and deuterium abundances. Conventionally, the energy
densities of new, dark and relativistic particle species are recast in terms
of a neutrino density to be constrained through the effective number of
neutrino species $N_{\rm eff}$. However, this recast is only possible when the
new energy density redshifts like radiation. Because our axions can either
behave as cold dark matter, radiation or be in a phase with a kination-like
scaling, $N_{\rm eff}$ does not provide a convenient language to cast our BBN
constraints in. Instead, we work directly in terms of energy densities, and
run a full numerical BBN calculation in the presence of such an additional
dark density component. Such a calculation can be done with the built-in
routines of the numerical code AlterBBN 2.2 Arbey:2018zfh .

The routine alter_standmod in AlterBBN 2.2 permits the addition of an
additional energy density component of the form

$\rho_{D}(T)=\zeta\rho_{\gamma}(T_{\rm BBN})\quantity(\frac{T}{T_{\rm
BBN}})^{n},$ (5.3)

where $\zeta=\rho_{D}(T_{\rm BBN})/\rho_{\gamma}(T_{\rm BBN})$ is the ratio of
the additional energy density to the photon density as measured at the
reference temperature $T_{\rm BBN}=1\,$\mathrm{MeV}$$ and $n$ is the scaling
exponent such that $\rho_{D}\propto a^{-n}$. Expressed in terms of the total
radiation density the parameter $\zeta$ is

$\zeta=\evaluated{\frac{\rho_{\theta}}{\rho_{\gamma}}}_{$\mathrm{MeV}$}=\frac{g_{\rho}(T_{\rm
BBN})}{2}\evaluated{\frac{\rho_{\theta}}{\rho_{\rm
rad}}}_{$\mathrm{MeV}$}\approx 5.29\evaluated{\frac{\rho_{\theta}}{\rho_{\rm
rad}}}_{$\mathrm{MeV}$}.$ (5.4)

If we assume that the ALP energy density scales cleanly with either $n=3$, $4$
or $6$ during BBN, we can then use this routine to calculate the produced
helium-4 and deuterium abundances. These abundances are then compared to the
most current PDG constraints Zyla:2020zbs , which at $1\sigma$ uncertainty are

$\displaystyle Y_{p}$ $\displaystyle=0.245\pm 0.003,$ (5.5)
$\displaystyle\evaluated{D/H}_{p}$ $\displaystyle=\quantity(2.547\pm
0.025)\times 10^{-5},$ (5.6)

where $Y_{p}$ and $D/H|_{p}$ are the helium-4 and deuterium fractions
respectively. Using these constraints with the AlterBBN routine alter_standmod
we obtain $2\sigma$ contstraints on ALP densities. Depending on the scaling of
the axion relic at $T\sim 1\,$\mathrm{MeV}$$, the constraints are the
following:

  * •

$\rho_{\theta}\propto a^{-3}$: Cold dark matter is not constrainted by BBN.

  * •

$\rho_{\theta}\propto a^{-4}$: Constrained by BBN if $\rho_{\theta}/\rho_{\rm
rad}|_{1\,$\mathrm{MeV}$}\gtrsim 3.2\times 10^{-2}$.

  * •

$\rho_{\theta}\propto a^{-6}$: Constrained by BBN if $\rho_{\theta}/\rho_{\rm
rad}|_{1\,$\mathrm{MeV}$}\gtrsim 1.9\times 10^{-1}$.

Only the last constraint is relevant for us. This corresponds to a lower bound
of

$T_{*}\gtrsim 20\mbox{ keV}.$ (5.7)

By using (4.10) and assuming that the field is still rolling at $T_{\rm BBN}$
the last condition above becomes equivalent to

$\frac{4\pi^{2}}{135}\quantity(\frac{g_{s}(T_{\ast})}{g_{s}(T_{\rm
BBN})})^{1/3}g_{s}(T_{\rm BBN})\quantity(\frac{YT_{\rm BBN}}{f})^{2}\lesssim
0.19.$ (5.8)

Expressing the yield in terms of the zero-temperature axion mass by using
(2.21) we find that the BBN constraint implies a bound on the zero-temperature
barrier height:

$\boxed{\Lambda_{b,0}\gtrsim 9\times
10^{-7}\,$\mathrm{GeV}$\times\quantity(\frac{h^{2}\Omega_{\Theta,0}}{h^{2}\Omega_{\rm
DM}})^{1/2}.}$ (5.9)

This result is independent of all other model parameters including the
temperature-dependence of the axion mass. We observe that this bound is always
stronger than constraint from structure formation.

Figure 11: Parameter space for ALP dark matter. Every point in the white
region can have the correct relic abundance to explain DM. We assume a
temperature-dependent axion mass with $\gamma=8$, see (2.29) for definition,
consistent with the QCD axion. The orange line separates the Standard and
Kinetic Misalignment regions. The region above the orange line can also be
reached without the initial kinetic energy, but by choosing the initial angle
very close to the top of the potential, i.e. Large Misalignment Mechanism.
Above the blue line, the fragmentation is efficient enough so that all the
energy density is transferred from the homogeneous mode to the fluctuations.
Above the green line, the fragmentation becomes efficient before the ALP gets
trapped by the potential; see Section 4.1 for the precise definitions. Above
the red line, the variance of the ALP angle becomes larger than unity before
the onset of fragmentation, so our calculation cannot be trusted. In the gray
region, the ALP field is rolling during BBN with a large enough kinetic energy
so that it spoils the BBN predictions, see Eq. (5.9). Thin lines correspond to
experimental projections. The bounds/projections on the axion-photon coupling
are translated into bounds on the axion decay constant by assuming a KSVZ-like
coupling given in (5.10). Orange constraints apply to any ALP while the green
ones assume the ALP is DM. Figure 12: Same as Figure 11 but now by assuming a
constant axion mass.

###  5.3 ALP parameter space with existing and future experiments

We now show the fragmentation regions of Section 4.1, and the model
constraints of Section 5 on the ALP parameter space together with a long list
of experiments constraints/projections of diverse nature in Figures 11 and 12.
Current constraints are included as filled regions while projections are
distinguished by having outlines only. In order to translate the
bounds/projections on the axion-photon coupling $g_{\theta\gamma\gamma}$ to
the bounds/projections on the axion decay constant $f$, we assumed a KSVZ-like
coupling:

$g_{\theta\gamma\gamma}^{\rm KSVZ}=\frac{\alpha_{\rm
EM}}{2\pi}\frac{1.92}{f}\approx\frac{2.23\times 10^{-3}}{f}.$ (5.10)

Figures 11 and 12 also include projections for constraints on the ALP-neutron
coupling, which are indicated with dashed lines. These assume a KSVZ-like
coupling of

$\displaystyle C_{\theta n}\approx-0.02,$ (5.11)

where the ALP-neutron coupling $C_{\theta n}$ is defined from

$\displaystyle\mathcal{L}\supset C_{\theta
n}\frac{\partial_{\mu}\theta}{2f_{a}}\bar{n}\gamma^{\mu}\gamma_{5}n.$ (5.12)

For other ALP couplings, the constraints and the projections of all the
experiments need to be adjusted with the exception of superradiance bounds
which depend on the axion mass and the axion decay constant directly. The
references for all the bounds/projections can be found in Appendix D and the
majority of the digitized bounds are sourced from the `AxionLimits` repository
AxionLimits .

The Kinetic Misalignment Mechanism (KMM) is not the only method to get ALP
dark matter with lower decay constants. Another possibility is the Large
Misalignment Mechanism (LMM) where one chooses an initial angle very close to
the top Zhang:2017dpp ; Arvanitaki:2019rax . The distinction between the
Standard and Large misalignment is arbitrary. In this work, we will define the
LMM region as the region above the orange line in Figures 11 and 12 without
the initial kinetic energy.

##  6 Gravitational waves

In the case of efficient parametric resonance, the ALP field at the end of
fragmentation can have significant anisotropic stress, and this stress can
source gravitational waves (GW). A precise study of this process is not the
topic of this paper, and will be studied in a future work. Here, we will
present a very rough estimate based on the method introduced in
Chatrchyan:2020pzh . A more detailed analysis can be found in Madge:2021abk .

We start by discussing the frequency of the gravitational waves. Today’s
frequency is related to the comoving momentum $k$ by

$\nu=\frac{1}{2\pi}\frac{k}{a_{0}}=\frac{1}{2\pi}\frac{k}{a_{\rm
emit}}\frac{a_{\rm emit}}{a_{0}},$ (6.1)

where $a_{\rm emit}$ is the scale factor at which the GW is emitted. By
assuming that emission happens around trapping $a_{\rm emit}\approx a_{\ast}$
we get

$\nu\approx\frac{m_{\ast}}{2\pi}\kappa_{\ast}\frac{a_{\ast}}{a_{0}}.$ (6.2)

We can calculate $a_{\ast}/a_{0}$ by solving (2.28). In this equation,
$\rho_{\Theta,0}$ is the ALP energy density today assuming no fragmentation.
The non-linear effects after the fragmentation can dilute the relic density by
a factor of $\mathcal{O}(1\text{--}10)$ Berges:2019dgr . Let $\mathcal{Z}$ be
this suppression factor so that $\rho_{\theta,0}=\mathcal{Z}\rho_{\Theta,0}$
is the correct ALP energy density with fragmentation. By requiring that
$\rho_{\theta,0}$ should match the dark matter density we get

$\frac{a_{\ast}}{a_{0}}=\quantity(\frac{3\pi}{8}\frac{\Omega_{\rm
DM}}{\mathcal{Z}}\frac{M_{\rm pl}^{2}H_{0}^{2}}{m_{0}m_{\ast}f^{2}})^{1/3}.$
(6.3)

We are interested to know the frequency at which the GW amplitude is peaked.
Since most of the fragmented modes have the momenta $\kappa_{\ast}\sim 1$, it
is likely that the GW spectrum will also be peaked at those momenta. Setting
$\kappa_{\ast}\sim 1$ in (6.2), and using (6.3) we obtain the peak frequency
in terms of model parameters as

$\nu_{\rm peak}\sim 8\times
10^{-11}\,$\mathrm{Hz}$\quantity(\frac{m_{\ast}}{m_{0}})^{2/3}\quantity(\frac{m_{0}}{10^{-16}\,$\mathrm{eV}$})^{1/3}\quantity(\frac{f}{10^{14}\,$\mathrm{GeV}$})^{-2/3}\mathcal{Z}^{-1/3}.$
(6.4)

Next we estimate the peak gravitational wave amplitude $\Omega_{\rm GW}(\nu)$.
It is defined as the fraction of energy density in gravitational waves per
logarithmic frequency:

$\frac{\rho_{\rm GW}}{\rho_{\rm crit}}=\int\differential{\ln\nu}\Omega_{\rm
GW}(\nu).$ (6.5)

After the emission, the energy density in gravitational waves dilutes as
radiation. Today’s GW amplitude is given in terms of the amplitude at emission
as

$\Omega_{\rm GW,0}=\frac{\rho_{\rm GW,0}}{\rho_{\rm crit,0}}=\frac{\rho_{\rm
GW,emit}}{\rho_{\rm crit,emit}}\frac{\rho_{\rm GW,0}}{\rho_{\rm
GW,emit}}\frac{\rho_{\rm crit,emit}}{\rho_{\rm crit,0}}=\Omega_{\rm
GW,emit}\quantity(\frac{a_{\rm emit}}{a_{0}})^{4}\quantity(\frac{H_{\rm
emit}}{H_{0}})^{2}.$ (6.6)

Let $k_{\rm peak}$ be the comoving momentum corresponding to the peak
frequency $\nu_{\rm peak}$. Ref. Chatrchyan:2020pzh gives the following
estimate for the peak amplitude at emission:

$\Omega_{\rm GW,emit}^{\rm peak}\sim\frac{64\pi^{2}}{3M_{\rm pl}^{4}H_{\rm
emit}^{2}}\frac{\rho_{\theta,\rm emit}^{2}}{\quantity(k_{\rm peak}/a_{\rm
emit})^{2}}\frac{\alpha^{2}}{\beta},$ (6.7)

where $\alpha\lesssim 1$ roughly measures the fraction of the energy stored in
the fluctuations, and $\beta\gtrsim 1$ is the typical logarithmic width of the
spectrum of fluctuations in momentum space. We set both of them to unity for
our estimates. Again we assume that the GW emission takes place at trapping.

Figure 13: Peak amplitude and frequency of gravitational waves induced by
axion fragmentation occurring in kinetic misalignment mechanism (KMM) versus
large misalignment mechanism (LMM) according to (Eq. 6.11). The lines assume
constant axion mass and should not be understood as GW spectra. The
predictions are compared to the sensitivity of future experiments. The
expression of $\nu_{\rm peak}$ in terms of the axion mass and decay constant
is given by Eq. 6.4, see contours in Fig. 14 and 15.

Then, the energy density of the ALP field at emission is

$\rho_{\theta,\rm emit}\approx 2m_{\ast}^{2}f^{2}.$ (6.8)

Also, the peak momentum becomes $k_{\rm
peak}=\kappa_{\ast}a_{\ast}m_{\ast}\sim a_{\ast}m_{\ast}$. Then (6.7) is
simplified to

$\Omega_{\rm GW,\ast}^{\rm
peak}\sim\frac{256\pi^{2}}{3}\quantity(\frac{m_{\ast}}{H_{\ast}})^{2}\quantity(\frac{f}{M_{\rm
pl}})^{4}.$ (6.9)

Evolving this amplitude until today by using (6.6) and (6.3) we obtain

$\Omega_{\rm GW,0}^{\rm peak}\sim 1.5\times
10^{-15}\quantity(\frac{m_{\ast}}{m_{0}})^{2/3}\quantity(\frac{m_{0}}{10^{-16}\,$\mathrm{eV}$})^{-2/3}\quantity(\frac{f}{10^{14}\,$\mathrm{GeV}$})^{4/3}\mathcal{Z}^{-4/3}.$
(6.10)

By combining this result with (6.4) we can obtain a simple relation between
the peak frequency and the peak amplitude:

$\Omega_{\rm GW,0}^{\rm peak}\sim
10^{-35}\quantity(\frac{m_{\ast}/m_{0}}{\quantity(\nu_{\rm
peak}/$\mathrm{Hz}$)\mathcal{Z}})^{2}.$ (6.11)

From this, we learn that ALP models with a constant mass have better prospects
for an observable gravitational-wave signal. Secondly, the models with a lower
peak frequency predict a larger GW amplitude. Finally, if there is an
additional dilution in the energy density due to the fragmentation, the
gravitational-wave amplitude is also enhanced. We show the contours of the
peak GW frequency calculated via (6.4), and the peak GW amplitude calculated
via (6.11) in the upper plots of Figure 14.

Figure 14: The peak frequency (Eq. 6.4) and peak amplitude (Eq. 6.11) contours
in the ALP parameter space, and the regions that can be probed by
gravitational waves (GW) induced by the fragmentation process by measuring the
$\mu$-distortions in the CMB Kite:2020uix . On the left plot we assume that
the fragmentation does not dilute the relic density, while the right plot
assumes a factor of $10$ dilution. Above the blue lines fragmentation is
complete, and it is likely that the efficient GW production only happens in
this region. The gray regions are excluded by the BBN constraints studied in
Section 5. GW prospects of kinetic fragmentation are thus severely constrained
by BBN.

A similar estimate has also been obtained in Madge:2021abk by a different
method, but their estimate is lower than ours by a factor of $\sim 20$. This
reference also presents a numerical calculation of the GW spectrum at the
linearized level, and confirms that their estimate predicts the peak amplitude
quite well. However, the GW amplitude can be enhanced significantly during the
non-linear phase which cannot be captured by the linearized analysis
Dufaux:2007pt . Therefore, one can interpret our estimate and the estimate of
Madge:2021abk as optimistic and conservative respectively.

The power-law-integrated-sensitivity for the SKA mission Janssen:2014dka can
reach up to $h^{2}\Omega_{\rm GW}\sim 3\times 10^{-16}$ at frequency $\nu\sim
2\times 10^{-9}\,$\mathrm{Hz}$$ assuming a signal-to-noise of $1$ and an
observation time of $20$ years Schmitz:2020syl . So even considering a
constant mass and assuming a suppression of factor of $10$, i.e.
$\mathcal{Z}=0.1$, our optimistic estimate tells us that the signal is barely
observable. Recently, Kite:2020uix did point out that the gravitational waves
with frequencies much smaller than the ones probed by SKA might be observable
by measuring the $\mu$-distortions in the CMB via the experiments such as
COBE/FIRAS 1994ApJ…420..439M ; Fixsen:1996nj , PIXIE 2011JCAP…07..025K ,
SuperPIXIE Kogut2019CMB , and Voyage 2050 Chluba:2019kpb . The COBE/FIRAS
experiment sets an upper limit on the $\mu$-distortions ($\mu<9\times
10^{-5}\,95$%CL), while the forecasted constraints for the experiments are
$\mu<3\times 10^{-8}$ for PIXIE, $\mu<7.7\times 10^{-9}$ for SuperPIXIE, and
$\mu<1.9\times 10^{-9}$ for Voyage 2050.

These upper bounds on the observable $\mu$-distortions can be translated into
lower bounds on the observable gravitational wave signal for a given frequency
Kite:2020uix . By using these bounds, we can estimate the potentially
observable regions in the ALP parameter space which we show in Figure 14. The
curve labeled “Voyage 2050 10x” assumes an upper limit of $\mu<1.9\times
10^{-10}$, see Kite:2020uix for details. We also show the critical lines
separating complete and incomplete fragmentation (blue line), as well as the
region which is excluded by the BBN constrains (gray region) discussed in
Section 5. Likely, the efficient GW production does happen only in the region
where the fragmentation is complete. We conclude that the BBN bound severely
constrains the parameter space which can be observable by GWs. The prospects
do improve if the fragmentation efficiently dilutes the relic density of ALPs,
since this will increase the energy budget for the GWs without overclosing the
universe. Precise estimations of the dilution factor, peak GW amplitude and
the GW frequency require a dedicated lattice analysis.

Figure 15: Analog of Figure 14 but for the case of large misalignment. Unlike
the case of kinetic misalignment there are no BBN constraints, however the
requirement that the ALPs are not too warm during the matter-radiation
equality gives a similar but weaker bound. In the purple shaded region on the
lower right, the initial angle is not close to the top, and the standard
misalignment is at play. The green line shows the parameter space where
$\absolutevalue{\pi-\Theta_{i}}\approx 10^{-10}$ above which the density
contrast of the ALP field becomes $\mathcal{O}(1)$, and it is expected that
complete fragmentation will happen above this line. It is likely that
efficient GW production will happen only in this region.

Since the BBN constraint is the main obstacle to get an observable
gravitational wave signal, one might wonder what happens in a model which is
immune to the BBN constraint. The Large Misalignment Mechanism is an example
of such a model since the energy density prior to the oscillations does not
redshifts as $a^{-6}$. Therefore, we did repeat the above analysis for the
Large Misalignment scenario, and obtained a result which is very close to
(6.11) except an $\mathcal{O}(1)$ factor which depends very mildly on the
value of $\absolutevalue{\pi-\Theta_{i}}$. We show the results in Figure 15.
Even though the BBN constraints are absent, there is still a similar but
weaker constraint that the ALPs should not be too warm at matter-radiation
equality in order to be consistent with the structure formation. The green
line shows approximately the critical initial angle
$\absolutevalue{\pi-\Theta_{i}}\approx 10^{-10}$ at which the density contrast
of the ALP field becomes $\mathcal{O}(1)$. Above this line, the initial angle
needs to be chosen closer to the top of the potential, and we expect that
complete fragmentation and the efficient GW production does happen in this
region. We see that the absence of the BBN constraint opens up a sizable
region which might be probed by future CMB surveys.

The ALP masses below $\sim 10^{-20}\,$\mathrm{eV}$$ are constrained by various
probes such as Lyman-$\alpha$ observations Rogers:2020ltq , and the galactic
rotation curves Bar:2018acw ; Bar:2019bqz ; Bar:2021kti . However, it is
possible that some of these contraints can be evaded as a consequence of the
fragmentation Leong:2018opi .

We close this section by stressing another much stronger signal of GW that can
arise in the Kinetic Misalignment Mechanism, and which is not related to
fragmentation. While in this work we have assumed a radiation-dominated
universe during the period when the equation of state of the axion scales as
$a^{-6}$, it is also possible that the axion energy density dominates
temporarily, inducing a kination era inside the radiation era. Such a kination
era enhances primordial signals of GW from inflation and cosmic strings
leading to striking peak features that are observable by upcoming
interferometers such as LISA and the Einstein Telescope Co:2021lkc ;
Gouttenoire:2021wzu ; Gouttenoire:2021jhk . The occurrence of a kination era
in specific UV completions and the precise parameter space region that lead to
observable GWs are worked out in Gouttenoire:2021jhk and will also be
discussed in Eroncel:2022a ; Eroncel:2022b .

##  7 Conclusion

In summary, the correct ALP relic abundance to explain dark matter can be
recovered naturally in the whole [axion mass $m$, decay constant $f$] plane.
The region which is usually discarded because it leads to insufficient dark
matter becomes open in the scenario of kinetic misalignment, in which
oscillations are delayed due to the initial velocity of the axion field. This
was pointed out in Co:2019jts ; Chang:2019tvx and the underlying framework was
extensively discussed in Co:2019wyp ; Co:2020jtv ; Co:2020xlh in some specific
UV completions. This is an extremely good news for a whole generation of
experiments as it provides them with a strong motivation for ALP dark matter
in the region of low values for the axion decay constant. It is in particular
exciting that experiments such as ALPSII which is about to start running at
DESY is sensitive to dark matter. It is also quite remarkable that IAXO can
probe the QCD axion as DM. What we have shown in this paper is that one can no
longer describe this regime in terms of the homogeneous zero-mode. In fact,
the axion field entirely fragments. We have provided a detailed analytical
derivation of the phenomenon.

One main model-independent consequence is a distinct prediction for dense
compact mini-clusters that is presented in detail in a companion paper
Eroncel:2022c . There is also a stochastic gravitational wave background
generated by axion fragmentation as we discuss in our Section 6. However, this
signal appears at extremely low frequency and is typically below the
sensitivity of Pulsar Timing Arrays. On the other hand, it would be visible in
CMB experiments such as Voyager.

Our results apply generically to any ALP, including the QCD axion. One key
question concerns the motivations for an initial axion velocity. To address
this, one should start with the complex scalar field from which the axion
originates as the angular direction. There is a whole class of models where
the axion receives a kick from the radial mode of the complex scalar field at
very early times in the cosmological evolution, see Co:2019jts ; Co:2019wyp ;
Co:2020jtv ; Co:2020xlh ; Gouttenoire:2021jhk . We will present the precise
determination of the viable parameter space of these models for ALP DM from
fragmentation in a companion paper Eroncel:2022a . Note that in a significant
region of parameter space, one can get a striking signal in gravitational
waves in these models from a kination era induced by the axion Co:2021lkc ;
Gouttenoire:2021wzu ; Gouttenoire:2021jhk . We do not discuss this in our
section 6 which only deals with the GW signal from the fragmentation effect,
ignoring the potential signal from a kination era, as in this paper we always
assume radiation domination at early times.

The main results of this paper, which are independent from any UV completion,
are summarised in the following expressions and figures:

  * •

The temperature $T_{\ast}$ when the field gets trapped by the barrier is given
in Eq. 2.30. Fig. 8 shows the contours of the trapping temperature.

  * •

The duration of the kination-like scaling underlying the kinetic misalignment
mechanism is constrained by (3.37).

  * •

Fig. 7 shows the contours of $m_{*}/3H_{*}$ at the time when the field gets
trapped, which is orders of magnitude larger compared to the unit value that
characterises the time when the field starts oscillating around its minimum in
the standard misalignment mechanism.

  * •

The amplification factor of a given mode depends crucially on the value of
$m_{*}/H_{*}$ as expressed in Eq. 3.13. It is also plotted in Fig. 9.

  * •

Fig. 11 and 12 show the contours of the different regimes of dark matter
production for respectively the constant axion mass and temperature-dependent
mass cases. The expressions for these region boundaries are given in (Eq. 4.2,
4.12, 4.17, 4.22). The BBN constraint, $T_{\ast}>20$ keV, translates into a
bound on the size of the zero-temperature barrier, see (5.9).

  * •

All phenomenological implications of fragmentation and in particular the
determination of the different regimes are controlled by the power spectrum
(3.33) of the mode functions (B.53).

  * •

Figure 9-left indicates that the impact of fragmentation on the relic
abundance prediction from kinetic misalignment is relatively weak. It
typically differs by a factor of order 1, we do not expect more than one order
of magnitude effect. This can be understood as the momentum of the produced
axions during fragmentation is of the same order as the axion mass (see Fig.
9-right), higher modes are not excited. A precise estimate of the ${\cal
O}(1)$ factor requires a lattice calculation which we will present in a future
work.

  * •

We estimated the GW signal in Fig. 13 and the parameter space that can be
probed this way in Fig. 14. We showed that there are better prospects for GW
from large misalignment in Fig. 15, something which had not been mentioned in
Arvanitaki:2019rax .

There are many interesting questions that deserve further investigation.
Primarily, the precise predictions for the QCD axion DM definitely deserve
better scrutiny, in particular using lattice calculations. Another compelling
aspect is the correlation with GW signatures that can be derived once a
specific UV completion is defined. Finally, the possibility to relate the halo
spectrum properties to the parameters of the axion potential is an exciting
opportunity in view of the upcoming observational prospects.

## Acknowledgements

This work is supported by the Deutsche Forschungsgemeinschaft under Germany
Excellence Strategy - EXC 2121 “Quantum Universe” - 390833306.

##  Appendix A Detailed discussion of the parametric resonance

In this Appendix we discuss in detail the solutions of (2.7). If we neglect
the expansion of the universe then (2.7) becomes

$\ddot{\theta}_{k}+\quantity(\frac{k^{2}}{a^{2}}+m^{2}\cos\Theta)\theta_{k}=0,$
(A.1)

where both the scale factor $a$ and the axion mass $m$ are constant. We define
a dimensionless time $t_{m}\equiv mt$, and dimensionless momentum variable
$\kappa\equiv k/am$ so that (A.1) takes a simpler form:

$\theta_{\kappa}^{\prime\prime}+\quantity(\kappa^{2}+\cos\Theta)\theta_{\kappa}=0,$
(A.2)

where primes denote derivatives with respect to $t_{m}$. The solutions of this
equation have been studied in Greene:1998pb in the case of oscillations after
trapping, i.e. $\varrho<1$. We review the method of getting the solutions, at
the same time we generalize it to obtain the solutions for the rolling axion,
$\varrho>1$.

We start by defining a new time variable $z$ by

$z(t_{m})\equiv\cos\Theta(t_{m}),$ (A.3)

where $\Theta$ is the solution of the homogeneous mode. Without expansion the
energy density $\varrho$ becomes a constant of motion. By taking the
derivative of $z(t_{m})$ by using the fact that $\varrho$ is conserved we
obtain

$\derivative{z}{t_{m}}=\pm\sqrt{2(1-z^{2})(2\varrho-1+z)},$ (A.4)

where the sign depends of the sign of $\sin\Theta$ and $\Theta^{\prime}$.
Plugging this result into (A.2) gives

$2(1-z^{2})(2\varrho-1+z)\theta_{\kappa}^{\prime\prime}+\quantity[1+2(1-2\varrho)z-3z^{2}]\theta_{\kappa}^{\prime}+\quantity(\kappa^{2}+z)\theta_{\kappa}=0,$
(A.5)

where primes are now denoting derivatives with respect to $z$. Note that $z$
is not single-valued for the whole oscillation. Thus this equation can
describe the solutions only for quarter of the period for $\varrho<1$, and
half of the period for $\varrho>1$. However we can still use this equation to
get the solutions in each patch, and then glue them together to get the full
solution.

Let $\theta_{\pm}(z)$ denote the two linearly independent solutions to (A.5),
where we have omitted the momentum subscript for cleaner notation. It turns
out that the linear combinations $\theta_{+}^{2}$, $\theta_{-}^{2}$, and
$\theta_{+}\theta_{-}$ obey a third order equation:

$2(z^{2}-1)(2\varrho-1+z)M^{\prime\prime\prime}+\quantity[9z^{2}-6(1-2\varrho)z-3]M^{\prime\prime}+2(z-1+2\varrho-2\kappa^{2})M^{\prime}-2M=0.$
(A.6)

This equation has a polynomial solution given by

$M(z)=z-1+2\varrho-2\kappa^{2}.$ (A.7)

The original equation we are trying to solve (A.2) is a _Hill differential
equation_ so according to the _Floquet’s theorem_ , the solutions must be of
the form

$\theta_{\pm}(t_{m};\kappa)=\psi_{\pm}(t_{m};\kappa)e^{\pm\mu_{\kappa}t_{m}}.$
(A.8)

Therefore the polynomial solution (A.7) should correspond to
$\theta_{+}\theta_{-}$ which implies

$\theta_{+}(z)\theta_{-}(z)=\mathcal{N}^{2}\quantity(z-1+2\varrho-2\kappa^{2}),$
(A.9)

where $\mathcal{N}$ is a normalization constant. To obtain the individual
solutions we note that the Wronskian of the system
$W\equiv\theta_{+}^{\prime}\theta_{-}-\theta_{+}\theta_{-}^{\prime}$ obeys

$W^{\prime}(z)=-\frac{1+2(1-2\varrho)z-3z^{2}}{2(1-z^{2})(2\varrho-1+z)}W(z)=-\derivative{z}\ln\sqrt{(1-z^{2})(2\varrho-1+z)}W(z).$
(A.10)

This can easily be solved as

$W(z)=\theta_{+}^{\prime}(z)\theta_{-}(z)-\theta_{+}(z)\theta_{-}(z)=\frac{c_{\kappa}\,\mathcal{N}^{2}}{(1-z^{2})(2\varrho-1+z)},$
(A.11)

where $c_{\kappa}$ is an integration constant coming from the integration of
(A.10) which we will determine shortly. By combining (A.8) and (A.11) we
obtain the following differential equations for $\theta_{\pm}$:

$2M(z)\derivative{\ln\theta_{\pm}}{z}=M^{\prime}(z)\pm\frac{c_{\kappa}}{\sqrt{(1-z^{2})(2\varrho-1+z)}}.$
(A.12)

The solutions are888Note that while $(1-z^{2})(2\varrho-1+z)$ is always non-
negative, $M(z)$ does cross zero at $z=1-2\varrho+2\kappa^{2}$. If this point
lies in the range of the integral, then the integral is understood as its
Cauchy principal value.

$\theta_{\pm}(z)=\mathcal{N}\sqrt{\absolutevalue{M(z)}}\exp(\pm\frac{c_{\kappa}}{2}\int^{z}\frac{\differential{z^{\prime}}}{M(z^{\prime})\sqrt{(1-z^{\prime
2})(2\varrho-1+z^{\prime})}}).$ (A.13)

The integration constant $c_{\kappa}$ can be obtained by plugging this
solution into (A.5). The result is

$c_{\kappa}^{2}=8\kappa^{2}(\varrho-\kappa^{2})(1-\varrho+\kappa^{2}).$ (A.14)

These coefficients determine the instability bands in the parametric
resonance. If $c_{\kappa}^{2}>0$, the Floquet exponent will be real, and
parametric resonance happens. Otherwise they will be imaginary, and mode
functions will only have oscillatory solutions. Therefore by using (A.14) we
can directly find the modes which are inside the instability bands. They are
given by

$\displaystyle\varrho-1<\kappa^{2}<\varrho,$ $\displaystyle\mbox{\quad
for\quad}\varrho>1,$ (A.15) $\displaystyle 0<\kappa^{2}<\varrho,$
$\displaystyle\mbox{\quad for\quad}\varrho<1.$ (A.16)

To find the value of the Floquet exponent we need to do a little bit more
work. We work out the cases before and after trapping separately.

##### Before trapping:

Without loss of generality we can assume that the homogeneous mode travels
from $\Theta=0$ to $\Theta=2\pi$ during one period. In the first half of the
period, $\Theta$ moves from $0$ to $\pi$, and $z$ decreases from $z=1$ to
$z=-1$. In this patch the exponentially growing solution is $\theta_{+}$. By
choosing the normalization factor $\mathcal{N}$ such that the solution is
unity initially, its value after half a oscillation is given by

$\theta_{\kappa}^{(1/2)}=\sqrt{\frac{\absolutevalue{M(-1)}}{\absolutevalue{M(1)}}}\exp(\frac{c_{\kappa}}{2}\int_{1}^{-1}\frac{\differential{z^{\prime}}}{M(z^{\prime})\sqrt{(1-z^{\prime
2})(2\varrho-1+z^{\prime})}})$ (A.17)

In the second half of the period, $\Theta$ moves from $\pi$ to $2\pi$ while
$z$ increases from $z=-1$ to $z=1$. Now the exponentially growing solution is
$\theta_{-}$. In order to glue the solutions the normalization $\mathcal{N}$
should be chosen such that the full solution is continuous. Then we obtain the
solution after a full oscillation as

$\begin{split}\theta_{\kappa}^{(1)}&=\theta_{\kappa}^{(1/2)}\sqrt{\frac{\absolutevalue{M(1)}}{\absolutevalue{M(-1)}}}\exp(-\frac{c_{\kappa}}{2}\int_{-1}^{1}\frac{\differential{z^{\prime}}}{M(z^{\prime})\sqrt{(1-z^{\prime
2})(2\varrho-1+z^{\prime})}})\\\
&=\exp(c_{\kappa}\int_{1}^{-1}\frac{\differential{z^{\prime}}}{M(z^{\prime})\sqrt{(1-z^{\prime
2})(2\varrho-1+z^{\prime})}}).\end{split}$ (A.18)

The Floquet exponent $\mu_{\kappa}$ can be obtained by matching this solution
to (A.8) and using the fact that $\psi_{\pm}$’s are periodic functions. We
then obtain

$\mu_{\kappa}^{(\varrho>1)}=\frac{c_{\kappa}}{\mathcal{T}_{>}}\int_{1}^{-1}\frac{\differential{z^{\prime}}}{M(z^{\prime})\sqrt{(1-z^{\prime
2})(2\varrho-1+z^{\prime})}},$ (A.19)

where $\mathcal{T}_{>}$ is the period of oscillation of the homogeneous mode
in physical time $t$ when $\varrho>1$. This can be derived from the
conservation of energy as

$\mathcal{T}_{>}=\frac{2}{m\sqrt{\varrho}}\textrm{K}(1/\sqrt{\varrho}).$
(A.20)

The integral in (A.19) can be expressed in a more useful form by changing the
integral path using the Cauchy’s integral theorem:

$-\int_{-\infty}^{1-2\varrho}\frac{\differential{z^{\prime}}}{M(z^{\prime})\sqrt{(1-z^{\prime
2})(2\varrho-1+z^{\prime})}}.$ (A.21)

By defining $\vartheta$ such that $z=-1/\sin^{2}\vartheta$ we obtain our final
result.

$\mu_{\kappa}^{(\varrho>1)}=\frac{c_{\kappa}}{\mathcal{T}_{>}}\int_{0}^{\arcsin(1/\sqrt{2\varrho-1})}\frac{\differential{\vartheta}}{1+\quantity(1-2\varrho+2\kappa^{2})\sin^{2}\vartheta}\frac{2\sin^{2}\vartheta}{\sqrt{\quantity(1+\sin^{2}\vartheta)\quantity[1+(1-2\varrho)\sin^{2}\vartheta]}}.$
(A.22)

##### After trapping:

Here we start the oscillation at the minimum $\Theta=0$. In the first half of
the oscillation $\Theta$ travels from $0$ to $2\arcsin(\sqrt{\varrho})$ while
$z$ decreases from $z=1$ to $z=1-2\varrho$. The growing solution is the
$\theta_{+}$ solution. Again normalizing the mode functions to unity initial
amplitude we find the solution after the first half of the oscillation as

$\theta_{\kappa}^{(1/2)}=\sqrt{\frac{\absolutevalue{M(1-2\varrho)}}{\absolutevalue{M(1)}}}\exp(\frac{c_{\kappa}}{2}\int_{1}^{1-2\varrho}\frac{\differential{z^{\prime}}}{M(z^{\prime})\sqrt{(1-z^{\prime
2})(2\varrho-1+z^{\prime})}}).$ (A.23)

In the second half of the oscillation $\Theta$ travels from
$2\arcsin(\sqrt{\varrho})$ back to $0$ while $z$ increases from $z=1-2\varrho$
to $1$. Now the growing solution is $\theta_{-}$. So after a full period the
solution is

$\theta_{\kappa}=\exp(c_{\kappa}\int_{1}^{1-2\varrho}\frac{\differential{z^{\prime}}}{M(z^{\prime})\sqrt{(1-z^{\prime
2})(2\varrho-1+z^{\prime})}}).$ (A.24)

The Floquet exponent can directly be read from this result as

$\mu_{\kappa}^{(\varrho<1)}=\frac{c_{\kappa}}{\mathcal{T}_{<}}\int_{1}^{1-2\varrho}\frac{\differential{z^{\prime}}}{M(z^{\prime})\sqrt{(1-z^{\prime
2})(2\varrho-1+z^{\prime})}},$ (A.25)

where the period $\mathcal{T}_{<}$ in this case is

$\mathcal{T}_{<}=\frac{2}{m}\textrm{K}(\sqrt{\varrho}).$ (A.26)

Applying the integral transformations that we did in deriving (A.22) yields to
the result

$\mu_{\kappa}^{(\varrho<1)}=\frac{c_{\kappa}}{\mathcal{T}_{<}}\int_{0}^{\pi/2}\frac{\differential{\vartheta}}{1+\quantity(1-2\varrho+2\kappa^{2})\sin^{2}\vartheta}\frac{2\sin^{2}\vartheta}{\sqrt{\quantity(1+\sin^{2}\vartheta)\quantity[1+(1-2\varrho)\sin^{2}\vartheta]}}.$
(A.27)

##### Final result:

Our final result for the Floquet exponents can be summarized as

$\mu_{k}=\sqrt{8\kappa^{2}\quantity(\varrho-\kappa)\quantity(1-\varrho+\kappa^{2})}\times\begin{cases}\mathcal{T}_{>}^{-1}(\varrho)\,\mathcal{I}\quantity(\arcsin\quantity(1/\sqrt{2\varrho-1})),&\varrho>1\\\
\mathcal{T}_{<}^{-1}(\varrho)\,\mathcal{I}\quantity(\pi/2),&\varrho<1\end{cases},$
(A.28)

where

$\mathcal{I}(\varphi)=\int_{0}^{\varphi}\frac{\differential{\vartheta^{\prime}}}{1+\quantity(1-2\varrho+2\kappa^{2})\sin^{2}\vartheta^{\prime}}\frac{2\sin^{2}\vartheta^{\prime}}{\sqrt{\quantity(1+\sin^{2}\vartheta^{\prime})\quantity[1+(1-2\varrho)\sin^{2}\vartheta^{\prime}]}}.$
(A.29)

We show a plot of the Floquet exponents together with the instability bands in
Figure 16. The boundaries of the instability bands are shown in white lines.
We can observe that the parametric resonance is most efficient around trapping
$\varrho\approx 1$, and for the modes $\kappa\sim 0.5$.

Figure 16: Floquet exponents as a function of the dimensionless energy density
$\varrho$, and the dimensionless momentum $\kappa$, using the analytical
result ((A.28)). The white lines denote the boundaries of the instability
bands as given in ((A.15)) and ((A.16)). We see that the parametric resonance
is most efficient around trapping $\varrho\approx 1$, and for the modes
$\kappa\sim 0.5$.

##  Appendix B Calculation of the adiabatic initial conditions

In this Appendix our primary goal is to derive the initial conditions for the
mode functions given by (3.33). In Appendix B.1 we briefly review the
cosmological perturbation theory for completeness, and at the same time
introduce our notation and conventions. We adopt the notation of
maggiore2008gravitational . The derivation of the initial conditions is
presented in Appendix B.2.

###  B.1 Brief review of cosmological perturbation theory

In order to study the perturbations, it is more convenient to use conformal
time which we shall denote by $\eta$. We employ the _Newtonian (conformal)
gauge_ where the perturbed metric takes the form

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

The perturbed stress-energy tensor can be put into the following form:

$\displaystyle T^{0}_{\;0}$ $\displaystyle=-\overline{\rho}-\delta\rho,$ (B.2)
$\displaystyle T^{i}_{\;0}$
$\displaystyle=-\quantity[\overline{\rho}+\overline{p}]v^{i},$ (B.3)
$\displaystyle T^{0}_{\;i}$
$\displaystyle=\quantity[\overline{\rho}+\overline{p}]v_{i},$ (B.4)
$\displaystyle T^{i}_{\;j}$ $\displaystyle=\quantity[\overline{p}+\delta
p]\delta^{i}_{j}+\quantity(\partial^{i}\partial_{j}-\frac{1}{3}\delta^{i}_{j}\nabla^{2})\Pi,$
(B.5)

where $\rho$ and $p$ are the energy density and the pressure respectively, and
quantities with an overline denote the background values. In these equations,
$v^{i}$ is the _peculiar velocity field_ , and $\Pi$ is the anisotropic
stress. The peculiar velocity field can be decomposed into a transverse and a
longitudinal part:

$v^{i}(\eta,\mathbf{x})=v^{T,i}(\eta,\mathbf{x})+\partial^{i}v(\eta,\mathbf{x}),\quad\partial_{i}v^{T,i}(\eta,\mathbf{x})=0.$
(B.6)

Here, $v$ is the _velocity potential_. Another useful quantity is the
_velocity divergence_ defined by

$\vartheta\equiv\partial_{i}v^{i}=\nabla^{2}v.$ (B.7)

It is more usual to express the perturbations in terms of _density contrast_
$\delta$, and the _sound speed_ $c_{s}$ defined respectively by

$\delta\equiv\frac{\delta\rho}{\overline{\rho}}\mbox{\quad
and\quad}c_{s}^{2}\equiv\frac{\delta p}{\delta\rho}.$ (B.8)

One should keep in mind that both $\delta$ and $c_{s}^{2}$ are gauge-dependent
quantities due to the fact that $\delta\rho$ and $\delta p$ are gauge-
dependent. To overcome this, one can work with gauge-independent variables. In
Newtonian gauge, the gauge-invariant generalizations of $\delta\rho$ and
$\delta p$ are given by

$(\delta\rho)^{\rm
GI}\equiv\delta\rho-3\mathcal{H}(\overline{\rho}+\overline{p})v\mbox{\quad
and\quad}(\delta p)^{\rm GI}\equiv\delta
p-3\mathcal{H}(\overline{\rho}+\overline{p})c_{s,(a)}^{2}v,$ (B.9)

where $\mathcal{H}=aH=a^{\prime}/a$ is the conformal Hubble parameter, and
$c_{s,(a)}^{2}$ is the _adiabatic speed of sound_ :

$c_{s,(a)}^{2}\equiv\evaluated{\frac{\delta
p}{\delta\rho}}_{s}=w+\frac{\overline{\rho}\dot{w}}{\dot{\overline{\rho}}},$
(B.10)

where $w=\overline{p}/\overline{\rho}$ is the equation of state. Then the
gauge-invariant versions of the density contrast and the sound speed are given
by

$\delta^{\rm GI}=\delta-3\mathcal{H}(1+w)v\mbox{\quad
and\quad}(c_{s}^{2})^{\rm GI}=\frac{(\delta p)^{\rm GI}}{(\delta\rho)^{\rm
GI}}.$ (B.11)

For modes deep inside the horizon $\delta$ and $c_{s}^{2}$ are approximately
equal to their gauge-invariant friends, so one can use them interchangebly.
However one should be careful about the super-horizon modes.

##### Einstein equations:

The linearized Einstein equations $\delta G_{\mu\nu}=8\pi G\delta T_{\mu\nu}$
can be used to relate metric perturbations to the perturbations in the stress-
energy tensor. This way, we obtain four equations:

$\displaystyle\nabla^{2}\Phi-3\mathcal{H}\quantity(\Phi^{\prime}-\mathcal{H}\Psi)$
$\displaystyle=-4\pi Ga^{2}\delta\rho,$ (B.12)
$\displaystyle\Phi^{\prime}-\mathcal{H}\Psi$ $\displaystyle=4\pi
Ga^{2}\quantity(\overline{\rho}+\overline{p})v,$ (B.13)
$\displaystyle\Phi+\Psi$ $\displaystyle=-8\pi Ga^{2}\Pi,$ (B.14)
$\displaystyle\Phi^{\prime\prime}+2\mathcal{H}\Phi^{\prime}-\mathcal{H}\Psi^{\prime}-(2\mathcal{H}^{\prime}+\mathcal{H}^{2})\Psi$
$\displaystyle=-4\pi Ga^{2}\quantity(\delta p+\frac{2}{3}\nabla^{2}\Pi).$
(B.15)

In the rest of this section we will assume that for all fluids, the
anisotropic stress does vanish $\Pi=0$999The anisotropic stress vanishes for
scalar fields at linear order in perturbation theory. Photons develop
anisotropic stress only during the matter era when they are sub-dominant. Only
free-streaming neutrinos have non-negligible anisotropic stress, but we do
neglect them in this work., so (B.14) implies $\Psi=-\Phi$.

##### Perturbations of a scalar field:

Consider an ALP field with the Lagrangian:

$\mathcal{L}=-\frac{f^{2}}{2}g^{\mu\nu}\partial_{\mu}\theta\partial_{\nu}\theta-V(\theta,T).$
(B.16)

where the potential is in general depends on the temperature $T$ of the
universe. The corresponding stress-energy tensor is given by

$f^{-2}T_{\mu\nu}=\partial_{\mu}\theta\partial_{\nu}\theta-
g_{\mu\nu}\quantity(\frac{1}{2}g^{\rho\sigma}\partial_{\rho}\theta\partial_{\sigma}\theta+\frac{1}{f^{2}}V(\theta,T)).$
(B.17)

We now write the ALP field $\theta(\eta,\mathbf{x})$ as the sum of the
homogenous component $\Theta(\eta)$ and fluctuations
$\delta\theta(\eta,\mathbf{x})$ which are assumed to be small. We also denote
the temperature fluctuations101010The contribution of the temperature
perturbations to the ALP density constrast has recently been pointed out in
Sikivie:2021trt ; Kitajima:2021inh . by $T=\overline{T}+\delta T$, where
$\overline{T}$ represents the background value. By expanding the stress-energy
tensor up to linear order in $\delta\theta$, $\Phi$ and $\delta T$, and
matching with the expressions (B.2)—(B.5), we obtain the zeroth order
quantities

$\displaystyle\overline{\rho}_{\theta}\equiv\rho_{\Theta}$
$\displaystyle=\frac{f^{2}}{2}\quantity(\frac{\Theta^{\prime}}{a})^{2}+V(\Theta,\overline{T})$
(B.18) $\displaystyle\overline{p}_{\theta}\equiv p_{\Theta}$
$\displaystyle=\frac{f^{2}}{2}\quantity(\frac{\Theta^{\prime}}{a})^{2}-V(\Theta,\overline{T}),$
(B.19)

and the first order quantities

$\displaystyle\delta\rho_{\theta}$
$\displaystyle=f^{2}\quantity[\frac{\Theta^{\prime}\delta\theta^{\prime}}{a^{2}}+\Phi\frac{\Theta^{\prime
2}}{a^{2}}+\frac{1}{f^{2}}\quantity(\evaluated{\partialderivative{V}{\theta}}_{\Theta,\overline{T}}\delta\theta+\evaluated{\partialderivative{V}{T}}_{\Theta,\overline{T}}\delta
T)]$ (B.20) $\displaystyle\delta p_{\theta}$
$\displaystyle=f^{2}\quantity[\frac{\Theta^{\prime}\delta\theta^{\prime}}{a^{2}}+\Phi\frac{\Theta^{\prime
2}}{a^{2}}-\frac{1}{f^{2}}\quantity(\evaluated{\partialderivative{V}{\theta}}_{\Theta,\overline{T}}\delta\theta+\evaluated{\partialderivative{V}{T}}_{\Theta,\overline{T}}\delta
T)]$ (B.21) $\displaystyle v_{\theta}$
$\displaystyle=-\frac{f^{2}\Theta^{\prime}\delta\theta}{a^{2}\quantity(\rho_{\Theta}+p_{\Theta})}=-\frac{f^{2}\Theta^{\prime}\delta\theta}{a^{2}\rho_{\Theta}(1+w_{\Theta})},$
(B.22)

where $w_{\Theta}$ is the equation of state of the background ALP field.
Perturbing the Klein-Gordon equation up to linear order in fluctuations gives
at zeroth order the evolution of the background field

$\Theta^{\prime\prime}+2\mathcal{H}\Theta^{\prime}+\frac{a^{2}}{f^{2}}\evaluated{\partialderivative{V}{\theta}}_{\Theta,\overline{T}}=0,$
(B.23)

where $\mathcal{H}=aH$ is the conformal Hubble factor. At the first order, we
obtain the evolution of the fluctuations:

$\delta\theta^{\prime\prime}+2\mathcal{H}\theta^{\prime}-\nabla^{2}\delta\theta+\frac{a^{2}}{f^{2}}\evaluated{\partialderivative[2]{V}{\theta}}_{\Theta,\overline{T}}\delta\theta=2\Phi\frac{a^{2}}{f^{2}}\evaluated{\partialderivative{V}{\theta}}_{\Theta,\overline{T}}-4\Phi^{\prime}\Theta^{\prime}-\frac{a^{2}}{f^{2}}\evaluated{\partialderivative{V}{\theta}{T}}_{\Theta,\overline{T}}\delta
T.$ (B.24)

In Fourier space this becomes

$\theta_{k}^{\prime\prime}+2\mathcal{H}\theta_{k}^{\prime}+k^{2}\theta_{k}+\frac{a^{2}}{f^{2}}\evaluated{\partialderivative[2]{V}{\theta}}_{\Theta,\overline{T}}\theta_{k}=2\Phi_{k}\frac{a^{2}}{f^{2}}\evaluated{\partialderivative{V}{\theta}}_{\Theta,\overline{T}}-4\Phi_{k}^{\prime}\Theta^{\prime}-\frac{a^{2}}{f^{2}}\evaluated{\partialderivative{V}{\theta}{T}}_{\Theta,\overline{T}}\delta
T.$ (B.25)

##### Curvature perturbations in radiation era:

The Fourier modes of the curvature perturbations $\Phi_{\mathbf{k}}$ have
exact solutions valid during the radiation era:

$\Phi_{\mathbf{k}}(\eta_{k})=3\Phi_{\mathbf{k}}(0)\quantity(\frac{\sin\eta_{k}-\eta_{k}\cos\eta_{k}}{\eta_{k}^{3}}),$
(B.26)

where $\eta_{k}\equiv k\eta/\sqrt{3}$. The initial conditions
$\Phi_{\mathbf{k}}(0)$ are imprinted by the inflation, and they are related to
the _comoving curvature perturbation_ $\mathcal{R}_{\mathbf{k}}(0)$ via

$\Phi_{\mathbf{k}}(0)=\frac{2}{3}\mathcal{R}_{\mathbf{k}}(0).$ (B.27)

Since the comoving curvature perturbations are generated by the quantum
fluctuations of the inflaton, they are stochastic variables. This means that
we can only calculate and measure the _power spectrum_ of them:

$P_{\mathcal{R}}^{\rm
in}(k)=\expectationvalue{\absolutevalue{\mathcal{R}_{\mathbf{k}}(0)}^{2}}.$
(B.28)

More specifically the experiments measure the _dimensionless power spectrum_ :

$\mathcal{P}_{\mathcal{R}}^{\rm
in}(k)=\frac{k^{3}}{2\pi^{2}}P_{\mathcal{R}}^{\rm in}(k)\equiv
A_{s}\quantity(\frac{k}{k_{\star}})^{n_{s}-1}.$ (B.29)

The Planck 2018 measurements (TT,TE,EE+lowE+lensing 68%) Aghanim:2018eyx are
consistent with

$A_{s}=2.1\times 10^{-9},\quad n_{s}=0.9649,\quad
k_{\star}=0.05\,\rm{Mpc}^{-1}.$ (B.30)

In this work we neglect the spectral tilt by setting $n_{s}=1$.

###  B.2 Adiabatic initial conditions in Kinetic Misalignment

In this section, we derive the initial conditions for the axion perturbations
that arise purely from adiabatic fluctuations. In other words, we neglect all
other types of isocurvature fluctuations to obtain a result which is as model-
independent as possible.

At early times, the axion mass can be neglected. Then the homogeneous mode
equations of motion (B.23) becomes

$\Theta^{\prime\prime}+2\mathcal{H}\Theta^{\prime}=0\quad\Rightarrow\quad\Theta^{\prime}\propto
a^{-2},$ (B.31)

where $\Theta$ is the axion homogeneous mode, and $\prime$ denotes derivative
with respect to conformal time. The mode function equation of motion (B.24)
takes the form

$\theta_{k}^{\prime\prime}+2\mathcal{H}\theta_{k}^{\prime}+k^{2}\theta_{k}=-4\Phi_{k}^{\prime}\Theta^{\prime},$
(B.32)

where $\Phi_{k}$’s are given by (B.26). In radiation domination
$\mathcal{H}=\eta^{-1}$. It will be more convenient to write this equation of
motion in terms of $\eta_{k}$. Doing so we obtain

$\derivative[2]{\theta_{k}}{\eta_{k}}+\frac{2}{\eta_{k}}\derivative{\theta_{k}}{\eta_{k}}+3\phi_{k}=-4\derivative{\Phi_{k}}{\eta_{k}}\derivative{\Theta}{\eta_{k}}\equiv\mathcal{F}(\eta_{k}).$
(B.33)

The right-hand-side of this equation is a forcing term which will induce
fluctuations even if they are zero initially. This equation has a homogeneous
solution

$\theta_{k}^{\rm
hom}(\eta_{k})=\frac{1}{\eta_{k}}\quantity[c_{1}\cos(\sqrt{3}\eta_{k})+c_{2}\sin(\sqrt{3}\eta_{k})],$
(B.34)

and a particular solution

$\theta_{k}^{\rm
par}=\theta_{s}(\eta_{k})\int^{\eta_{k}}\differential{\eta_{k}^{\prime}}\frac{\theta_{c}(\eta_{k}^{\prime})\mathcal{F}(\eta_{k}^{\prime})}{W\quantity[\theta_{c}(\eta_{k}^{\prime}),\theta_{s}(\eta_{k}^{\prime})]}-\theta_{c}(\eta_{k})\int^{\eta_{k}}\differential{\eta_{k}^{\prime}}\frac{\theta_{s}(\eta_{k}^{\prime})\mathcal{F}(\eta_{k}^{\prime})}{W\quantity[\theta_{c}(\eta_{k}^{\prime}),\theta_{s}(\eta_{k}^{\prime})]}.$
(B.35)

Here $\theta_{c}=\cos(\sqrt{3}\eta_{k})/\eta_{k}$ and
$\theta_{s}=\sin(\sqrt{3}\eta_{k})/\eta_{k}$ are the homogeneous solutions,
and $W$ is their Wronskian which is $\sqrt{3}/\eta_{k}^{2}$.

Let $\eta_{i}=\mathcal{H}_{i}^{-1}$ denote an arbitrary time during which the
axion homogeneous mode scales as kination, i.e. $\rho_{\Theta}\propto a^{-6}$.
We also define

$\eta_{k,i}\equiv\frac{k\eta_{i}}{\sqrt{3}}$ (B.36)

which tells whether a given mode is super- or sub-horizon at $\eta=\eta_{i}$.
Then the forcing term $\mathcal{F}(\eta_{k})$ is given by

$\mathcal{F}(\eta_{k})=-12\Phi_{k}(0)\frac{\dot{\Theta}_{i}}{H_{i}}\frac{\eta_{k,i}}{\eta_{k}^{2}}\quantity(\frac{\sin(\eta_{k})}{\eta_{k}^{2}}+\frac{3\cos(\eta_{k})}{\eta_{k}^{3}}-\frac{3\sin(\eta_{k})}{\eta_{k}^{4}}).$
(B.37)

So the full solution becomes

$\theta_{k}(\eta_{k})=\frac{1}{\eta_{k}}\quantity{\cos(\sqrt{3}\eta_{k})\quantity[c_{1}-\mathcal{I}_{s}(\eta_{k})+\mathcal{I}_{s}(\eta_{k,i})]+\sin\quantity(\sqrt{3}\eta_{k})\quantity[c_{2}+\mathcal{I}_{c}(\eta_{k})-\mathcal{I}_{c}(\eta_{k,i})]},$
(B.38)

where $\mathcal{I}_{c,s}$ are indefinite integrals

$\mathcal{I}_{c,s}(\eta_{k})\equiv\int^{\eta_{k}}\differential{\eta_{k}^{\prime}}\frac{\theta_{c,s}(\eta_{k}^{\prime})\mathcal{F}(\eta_{k}^{\prime})}{W\quantity[\theta_{c}(\eta_{k}^{\prime}),\theta_{s}(\eta_{k}^{\prime})]},$
(B.39)

where $\theta_{c}(\eta_{k})=\cos(\sqrt{3}\eta_{k})/\eta_{k}$, and
$\theta_{s}(\eta_{k})=\sin(\sqrt{3}\eta_{k})/\eta_{k}$ are the homogeneous
solutions of (B.33), and $W=\sqrt{3}/\eta_{k}^{2}$ is their Wronskian. Even
though these integrals can be performed analytically, their full expressions
are given by trigonometric integrals, and are too complicated to be
informative. At large $\eta_{k}$, both integrals decay as $\eta_{k}^{-3}$.
This means that in the deep sub-horizon limit, we can neglect the
$\mathcal{I}_{c,s}(\eta_{k})$ terms in (B.38), and keep only the
$\mathcal{I}_{c,s}(\eta_{k,i})$ terms.

In order to get the full solution one also needs to specify the coefficients
$c_{1}$ and $c_{2}$. These are determined by the solution of the mode
functions at $\eta=\eta_{i}$. To calculate them, we will assume that the mods
we are interested are super-horizon at $\eta_{i}$, which implies that
$\eta_{k,i}\ll 1$. The assumption is justified towards the end of Section 3.3.
Then it is possible to determine both $c_{1}$ and $c_{2}$ by using the
adiabacity conditions. The adiabatic perturbations are defined such that the
local state of matter at some spacetime point $(\eta,\mathbf{x})$ of the
perturbed universe is the same as the background universe at some slightly
different time. For a universe filled with multiple fluids, the adiabatic
perturbations are induced by a common, local shift in time of all background
quantities Baumann:2018muz . Therefore at early times, the perturbations of
all species can be described by

$\delta\rho_{i}(\eta,\mathbf{x})=\overline{\rho}_{i}(\eta+\delta\eta(\mathbf{x}))-\overline{\rho}(\eta)\approx\overline{\rho}_{i}^{\prime}\delta\eta(\mathbf{x}).$
(B.40)

Using the continuity equation
$\overline{\rho}_{i}^{\prime}=-3\mathcal{H}(1+w_{i})\overline{\rho}_{i}$, this
equation implies

$\frac{\delta_{i}}{1+w_{i}}=\frac{\delta_{j}}{1+w_{j}}$ (B.41)

for all species $i$ and $j$. Let us consider the adiabatic perturbations in
the rotating axion $\delta_{\theta}$, and in the radiation $\delta_{r}$. At
early times, the axion equation of state is $w_{\theta}=1$, while for
radiation it is $w_{r}=1/3$. Then, two perturbations should be related to each
other by

$\delta_{\theta}=\frac{3}{2}\delta_{r}.$ (B.42)

Since all perturbations are comparable, the total density perturbation
$\delta\rho\equiv\sum_{i}\overline{\rho}_{i}\delta_{i}$ is dominated by the
species that carry the dominant energy density. Since the universe is
radiation dominated at early times, we can assume
$\delta\rho\approx\overline{\rho}_{r}\delta_{r}$. Using (B.12) and the
Friedmann equation $3\mathcal{H}^{2}\approx 8\pi Ga^{2}\overline{\rho}_{r}$,
one can show that $\delta_{r}\approx 2\Phi$ in the super-horizon limit. This
implies that the adiabatic initial conditions for the rotating axions should
have the following super-horizon limit:

$\delta_{\theta}\approx 3\Phi.$ (B.43)

We now present an alternative derivation of this result, which will also give
us the correct initial conditions for the axion mode functions. Since
adiabatic initial conditions arise due to the fluctuations in the total energy
density, they should vanish in a time slicing where there is no
perturbatiation in the energy density $\evaluated{\delta\rho}_{\rm
unif}=0$111111Note that this is true only for the super-horizon adiabatic
modes. Once they enter the horizon they don’t vanish even in the uniform-
density gauge. Classification of perturbations as adiabatic or isocurvature
are done by specifying their super-horizon behavior.. This is known as the
_slicing of uniform energy density_. From a generic slicing, one can go this
slicing via a time displacement given by Riotto:2002yw

$\eta\rightarrow\eta+\xi,\quad\xi=\frac{\delta\rho}{\overline{\rho}^{\prime}}.$
(B.44)

Under this time displacement, the perturbations $\delta S$ of a scalar
quantity $S$ transforms as $\delta S\rightarrow\delta
S-\overline{S}^{\prime}\xi$. Since both the axion angle $\theta$ and the axion
energy density $\rho_{\theta}$ are scalar quantities, we can write the
following relations between the values of these quantities in the uniform
density gauge, and the Newtonian gauge:

$\displaystyle 0$ $\displaystyle=\evaluated{\delta\theta}_{\rm
unif}=\evaluated{\delta\theta}_{\rm new}-\Theta^{\prime}\xi$ (B.45)
$\displaystyle 0$ $\displaystyle=\evaluated{\delta\rho_{\theta}}_{\rm
unif}=\evaluated{\delta\rho_{\theta}}_{\rm new}-\rho_{\Theta}^{\prime}\xi.$
(B.46)

Since the dominant energy density is in radiation, we can write

$\xi=\frac{\delta\rho}{\overline{\rho}^{\prime}}\approx\frac{\delta\rho_{r}}{\overline{\rho}_{r}^{\prime}}=-\frac{\delta\rho_{r}}{3\mathcal{H}\overline{\rho}_{r}(1+w_{r})}=-\frac{\delta_{r}}{4\mathcal{H}}\approx-\frac{\Phi}{2\mathcal{H}}.$
(B.47)

By plugging this result into (B.45) and (B.46), and taking the Fourier
transform we find the following super-horizon initial conditions for the axion
mode functions in the Newtonian gauge:

$\displaystyle\lim_{\eta_{k}\rightarrow 0}\theta_{k}(\eta_{k})$
$\displaystyle\approx-\frac{1}{2}\Phi_{k}(0)\frac{\Theta^{\prime}}{\mathcal{H}}=-\frac{1}{2}\Phi_{k}(0)\frac{\dot{\Theta}}{H},$
(B.48) $\displaystyle\lim_{\eta_{k}\rightarrow 0}\delta_{\theta,k}(\eta_{k})$
$\displaystyle=2\quantity[\frac{\theta_{k}^{\prime}(\eta)}{\Theta^{\prime}(\eta)}+\Phi_{k}(\eta_{k})]\approx
3\Phi_{k}(0)$ (B.49)

We see that (B.49) agrees with (B.43). Now we can take the full solution
(B.38), evaluate at $\eta_{k}=\eta_{k,i}$, take the $\eta_{k,i}\rightarrow 0$
limit, and find the coefficients $c_{1}$ and $c_{2}$ by matching with the